Remote sensing of tropospheric turbulence using GPS radio occultation

Radio occultation (RO) measurements are sensitive to the small-scale irregularities in the atmosphere. In this study, we present a new technique to estimate tropospheric turbulence strength (namely, scintillation index) by analyzing RO amplitude fluctuations in impact parameter domain. GPS RO observations from the COSMIC (Constellation Observing System for Meteorology, Ionosphere, and Climate) satellites enabled us to calculate global maps of scintillation measures, revealing the seasonal, latitudinal, and longitudinal characteristics of the turbulent troposphere. Such information are both difficult and expensive to obtain especially over the oceans. To verify our approach, simulation experiments using the multiple phase screen (MPS) method were conducted. The results show that scintillation indices inferred from the MPS simulations are in good agreement with scintillation measures estimated from COSMIC observations.


Introduction
Atmospheric turbulence associated with fluctuation of temperature, humidity, and water vapor are prevalent in the tropospheric 10 region. Irregularities in the turbulence cause the index of refraction of the tropospheric medium to fluctuate. Electromagnetic signals transmitted from Global Navigation Satellite System (GNSS) satellites (for example), carrying communication and navigation information, propagate through the turbulent troposphere. The spatial changes of the index of refraction introduce irregular fluctuations in the intensity and phase of the traversing electromagnetic signals by causing scintillation (Wheelon, 2004). Signal scintillation can affect the performance of satellite communication and navigation systems such as Global Posi- 15 tioning System (GPS). On the other hand, scintillation characteristics inferred from global GPS signal observations (employed here) are valuable resources to study atmospheric turbulence properties.
In this paper, we have employed data analysis and model simulation to investigate and quantify the effects of tropospheric turbulence on L-band signals propagating from a GPS satellite to a low earth orbit (LEO) satellite such as COSMIC (Constellation Observing System for Meteorology, Ionosphere, and Climate) (Anthes et al, 2008). We estimated global distribution 20 of the turbulence strength using a scintillation parameter (scintillation index) from COSMIC Radio Occultation (RO) measurements. To understand and quantify the observed scintillation, we have simulated the effect of tropospheric turbulence on L-band signals propagating through multiple phase screens (MPS). In the MPS model runs, the phase screens are assumed to be various realizations of random perturbations of index of refraction profiles through which electromagnetic waves are propagating (Knepp, 1983).

Global Scintillation Maps Inferred from COSMIC RO Measurements
This section presents investigation of the effect of tropospheric turbulence on L-band propagation utilizing RO observations from a GPS to a COSMIC satellites radio links. We estimated the impact of turbulence strength on L-band signals in terms 5 of scintillation index. The COSMIC satellites provide a significant number of RO profiles (up to about 2000 profiles per day) observed by the six micro-satellites covering the entire globe. Utilizing the RO profiles, we were able to estimate the global distribution of the effect of scintillation on GPS signals. The technique provides valuable scintillation data especially over the oceans where ground-based measurements are both difficult and expensive to perform. RO data were first used to determine the intensity and location of turbulent regions (Cornman et al, 2012). Our study differs in that we aim to study the global 10 climatology of tropospheric turbulence. In addition, we suggest that a different observable than raw amplitude and phase measurements is more effective in estimating scintillation index.

COSMIC signal amplitude measurements
The basic observations of RO soundings on COSMIC are time series of amplitude (the 1-second voltage signal-to-noise ratio) and phase of L-band signals transmitted by a GPS satellite (e.g., Hajj et al (2002)). As the radio signal propagates through 15 the troposphere, the amplitude of the signal (raw amplitude) suffers from the effects of defocusing, multipath propagation, and diffraction effects. Therefore, the raw amplitude fluctuation does not truly represent the effect of turbulence. To suppress amplitude fluctuations due to these effects, the canonical transform (CT) has been applied on the received complex signal to transform it from the time domain to the impact parameter domain (Gorbunov, 2002). In RO retrieval processing, the CT phase is considered the important quantity since it is used to retrieve the bending angle profile and subsequently the refractivity 20 profile, which are the primary retrieval variables. In our study, however, the focus is on the CT amplitude. We assume that CT amplitude fluctuations are dominated by small-scale irregularities caused by turbulent processes.
We note that CT operates under the assumption of a spherically symmetric atmosphere where each ray is uniquely identified by its impact parameter. Thus the presence of mesoscale or large-scale horizontal inhomogeneity can result in fluctuations in the CT amplitudes; however, these tend to occur at a larger spatial scales than the turbulent effects being considered so that 25 their contribution to the scintillation estimates is expected to be small.

Scintillation Effects on COSMIC RO Signals
The scintillation index σ I can be viewed as a measure of the effect of tropospheric turbulence on L-band propagation having signal intensity I. σ I can be described as a normalized standard deviation of intensity I fluctuation: where ... stands for expected value and I is the CT intensity, which is the square of the CT amplitude. The expected value is computed using data in 120 m intervals. In the global map of scintillation estimates presented in Figures 1, the scintillation index is an altitudinal average of σ I between altitudes 1 to 4 km. Figure 2 shows an average of σ I between 4 and 8 km. In this study, we restrict our results to above 1 km due to possible data quality issues affecting the near-surface retrievals (Ao et al, 2012). The global scintillation maps 7. The turbulence strength difference between summer and winter seasons over the Antarctica is small compared to the turbulence strength difference between summer and winter seasons over the Arctic. Similar results based on radiosonde, satellite, and atmospheric reanalyses observation were reported (Serreze et al, 2012;Bromwich et al, 2007). Compared to other regions, the Arctic and Antarctic have the least observations available. The RO inferred map presented in Figures   15 1 and 2, therefore, helps to build and enrich the polar region database.

Scintillation Index Estimates
8. For all seasons, σ I changes significantly with season in the North-East Pacific compared to South-East Pacific.

The Index of Refraction Structure Parameter
The index of refraction structure parameter C 2 n is valuable for investigating propagation of electromagnetic waves in the atmosphere. Amplitude and phase of the waves propagating through the atmosphere get degraded as C 2 n is intensified (Andreas 20 et al, 1989). Unfortunately, measuring C 2 n is expensive and difficult to perform, especially over the oceans. Due to that C 2 n profiles are available from only a few locations over the globe ((Andreas et al, 1989;Wheelon, 2004) and the references). RO offers the possibility of filling these gaps. Here, we present C 2 n profiles inverted from CT amplitude scintillation incurred on COSMIC RO observations and the Rytov variance (equation 2).
Assuming that the wavelength of propagation is small compared to the scale of index of refraction fluctuations, σ I can be 25 expressed in terms of turbulence properties (C 2 n ) of the medium by the Rytov variance. The Rytov variance assumes wave propagation in weak fluctuation regimes (σ 2 I < 1) caused by turbulence processes having characteristics of Kolmogorov spectra. Since the amplitude of fluctuations of the RO signals considered are characterized by σ 2 I < 1, the Rytov variance is valid for the present analysis. Note that the Rytov variance was used to describe characteristics of weak scintillation of radio and microwave propagation in the atmosphere (Clifford and Strohbehn, 1970;Woo and Ishimaru, 1973;Ishimaru, 1978;Martini 30 et al, 2006;Blaunstein and Christodoulou, 2007). For plane waves, σ 2 I is proportional to C 2 n (Andrews et al, 1999): κ • is wavenumber of the electromagnetic wave, λ is wavelength, and L is the propagation path length between a transmitter and a receiver through the turbulent medium.
The C 2 n profile can be calculated from σ 2 I profile for each occultation by making use of equation 2 with λ = 0.2 m and κ o = 31.4 m −1 . Figure 3A, 3B, and 3C show the distribution of σ 2 I averaged over 1 to 8 km in the tropics, mid-latitude, and high-latitude troposphere, respectively. The histograms reveal that about 95% of the σ 2 I values are less than 1 in Figure 3A. The 5 distribution of σ 2 I shows values less than 1 more 95% of the cases.. Since measured σ 2 I < 1 in Figures 1 and 2 predominantly, the application of the Rytov variance is justified for estimating C 2 n profile. To compute C 2 n profiles from σ 2 I profiles, we still have to know the value of an effective signal propagation length L in the lower troposphere corresponding to RO geometry. For that purpose, L has been estimated using a multi phase screen (MPS) model simulation technique described in Section 3. We use the MPS model to calculate the scintillation index σ 2 model I by varying the value of C 2 n (the MPS model runs use index of 10 refraction as input derived from Kolmogorov spectra, equations 3 and 4). We used the following procedures in the MPS model to estimate L: The relationship between the scintillation index σ 2 model I and C 2 n is linear (equation 2). Using linear regression, we determine the slope in the σ 2 model I and C 2 n relationship. The propagation path length L is estimated from the slope, and have an average value of L = 650 km in the lower troposphere. We recognize that L is not exactly a constant and would vary with altitudes and differ from sounding to sounding; however, we expect using a single average L to estimate C 2 n for all occultations 15 to be a good approximation to first order. in Figures 4A, 4B and 4D (Central Pacific, New Mexico, and Arctic (North of Alaska)), the C 2 n profiles show clear seasonal behavior (summer-winter contrast). C 2 n is larger in the July compared to January (New Mexico and Arctic (North of Alaska), Figure 4B & 4D). It is larger in January compared to July (South East Pacific, Figure 4C). On the other hand, C 2 n is only a little larger in July compared to January (Central Pacific, Figure 4A). The mean C 2 n value (0.25 × 10 −14 m −2/3 ), reported in (Wheelon (2004) and the references therein), agrees with an approximate mean summer values of C 2 n displayed in Figures 4B   25 and 4D.

Multiple Phase Screen (MPS) Model
The merit of the phase screen technique to find solution of electromagnetic wave propagation through turbulent medium has been described in the literature (e.g., (Knepp, 1983;Knepp and Lickisch, 2009)). In the context of this paper, the MPS model 30 is applied to simulate and quantify the effect of tropospheric turbulence (and the resulting index of refraction variations) on Lband propagation. The MPS technique employed uses the geometry of L-band propagation in a radio link from a GPS satellite (incident wave) to a low-Earth orbit satellite (receiver, observation screen) (Sokolovskiy, 2001).
The MPS technique involves finding the solution of the parabolic wave equation (Knepp, 1983;Knepp and Lickisch, 2009).
The solution comprises application of numerical and Fourier transform techniques repeatedly for waves propagating through MPSs (Knepp, 1983;Knepp and Lickisch, 2009) idealized by parallel planes (phase changing screens). The medium of signal propagation is divided into a number of thin layers, each modeled as phase changing screens encompassed by a free space.
Accordingly, we seek solution of the parabolic equation by performing successive iteration of the propagation calculations 5 from one screen to the next, and ultimately to the observation screen.

Index of Refraction Profiles: Model Input
The main input for the MPS model run is an index of refraction profile of the lower troposphere. The phase screens are constructed as random perturbations of an exponential background refractivity profile (first term in equation 3): where the index of refraction profile n(y) is calculated from the refractivity profile using N (y) = 10 6 [n(y) − 1], N (0) = N • = 320, and the scale height H s = 7 km. The resulting background refractivity is shown in Figure 5A. In the refractivity perturbation models (Ñ (y)), tropospheric turbulence is taken to be nonzero only up to 8 km.Ñ (y) is an inverse Fourier transform of the Kolmogorov spectra (equation 4) modulated by a Gaussian random number with zero mean and standard 15 deviation one (Knepp, 1983). In order to make the refractivity perturbations realistic, the refractivity profiles on the phase changing screens are specified by the characteristics of Kolmogorov spectra, Φ(κ) = 0.033C 2 n κ − 11 3 , κ = 2π y (4) Figure 5B shows example refractivity profiles calculated using equation 4 (an input for the MPS model runs).

Input Parameters for MPS model runs 20
The input parameters for MPS model runs were: Number of phase screens = 4000; vertical spacing = 1 m; spacing between phase screens = 1 km; a Gaussian random number initializes each phase screen differently; on each phase screen, each altitude is initialized differently.

Scintillation Index Profiles: MPS model
The MPS model runs were performed for two cases where: (a) C 2 n = 1.0 × 10 −13 m −2/3 and (b) C 2 n = 5.0 × 10 −14 m −2/3 . For 25 each case, 50 random realizations of the MPS runs were used to calculate CT amplitude profiles. The average CT intensity profiles (average of the 50 CT intensity profiles) were then used to calculate σ I profiles every 50 m in the lower troposphere (shown in Figure 5C).
For comparison and validation purposes, Figure 5D  with the COSMIC RO inferred σ I . Figures 5C and 5D show a very good agreement between the MPS model σ I -COSMIC RO σ I over the Central Pacific (January and July 2008) and New Mexico (July 2008). Figure 5E presents C 2 n estimates derived from σ 2 I profiles (MPS model, Figure 5C). The Rytov variance has been used to invert for the C 2 n profiles. Figure 5E reproduces the input C 2 n values fairly well. The C 2 n profiles shown in Figure 5E have similarities with the C 2 n profiles shown in Figure 4 validating C 2 n profiles derived from COSMIC RO data. In Figure 5E, the blue vertical line shows C 2 n = 5.0 × 10 −14 m −2/3 and 5 the brown vertical line shows C 2 n = 1.0 × 10 −14 m −2/3 . The decrease of C 2 n starting ∼6 km is believed to be due to the finite vertical extent of the turbulence model in the MPS simulation.

Conclusions
We have used (i) radio occultation observations on board the COSMIC satellites, and (ii) multiple phase screen model calculations to investigate and quantify the effect of tropospheric turbulence on L-band propagation. Instead of regular amplitude 10 and phase data, we have used the CT amplitudes in estimating the scintillation indices from each occultation. This has the advantage of removing signal fluctuations due to atmospheric multipath and diffraction from sharp vertical layers.
Global maps of scintillation measures across different seasons were obtained from one year of COSMIC RO data. The resulting global scintillation maps reveal very strong seasonal dependence, with the northern hemisphere summer exhibiting relatively large turbulent activities compared to the southern hemisphere summer. Irrespective of seasons, the tropical regions 15 are generally characterized by a relatively large scintillation index. The maps also show clear ocean-continent contrast of the scintillation estimates in the Southern East pacific, South America, and South Atlantic regions. The scintillation estimates appear to be positively correlated with water vapor, precipitation, and convection. We have also presented C 2 n profiles estimated from COSMIC RO data using the Rytov variance method for weak scintillation. This represents the first ever observational estimates of global tropospheric turbulence strength. While certain simplifications have been used in this initial study, the 20 results are encouraging, and future work can be done to refine the results and cross-validate with other observations. We have also performed numerical simulations of radio propagation through a random phase changing screen (in which refractivity profiles were specified by the Kolmogorov spectra). Scintillation index profiles inferred from the MPS technique are in a reasonable agreement with scintillation index profiles inferred from COMSIC RO data and provide confidence to our