A Statistical Analysis of Sporadic-E Characteristics Associated with GNSS Radio Occultation Phase and Amplitude Scintillations

: Statistical GNSS-RO measurements of phase and amplitude scintillation are analyzed at the mid-latitudes in the local summer for a 100 km altitude. These conditions are known to contain frequent sporadic-E, and the S 4 - σ φ trends provide insight into the statistical distributions of the sporadic-E parameters. Joint two-dimensional S 4 - σ φ histograms are presented, showing roughly linear trends until the S 4 saturates near 0.8. To interpret the measurements and understand the sporadic-E contributions, 10,000 simulations of RO signals perturbed by sporadic-E layers are performed using length, intensity, and vertical thickness distributions from previous studies, with the assumption that the sporadic-E layer acts as a Gaussian lens. Many of the key trends observed in the measurement histograms are present in the simulations, providing a key for understanding the complex mapping between layer characteristics and impacts on the GNSS-RO signals. Ad-ditionally, the inclusion of Kolmogorov turbulence and a diffusion-limited threshold on the lens strength/(vertical thickness) 2 ratio helps to make the layers more physically realistic and improves agreement with the observations.


Introduction
With modern day society's strong reliance on GNSS for both navigation and timing, scintillation caused by ionospheric irregularities can induce loss of lock and increase navigation errors [1][2][3], thereby seriously impacting our technological infrastructure [4].While ground-based GNSS receivers are most severely perturbed by equatorial plasma bubbles (EPBs) [5], sporadic-E (E s ) is the leading cause of spaced-based GNSS cycle slips from an E-region ionospheric irregularity [6].Sporadic-E manifests as enhanced ion layers caused by metallic ion convergence through neutral wind shears [7,8].The sporadic-E layers are typically referred to as "clouds" because of their non-uniform, patchy nature.The E s impact on radio wave propagation, in terms of over-the-horizon radar and global positioning system (GPS) timing for low Earth orbit (LEO) satellites, relies not only on the presence of sporadic-E but also on the specific characteristics.
Many methods for measuring or characterizing E s have been proposed and implemented in previous studies, including incoherent scatter radars (ISRs) [9], ionosondes [10][11][12][13][14][15], and global navigation satellite system (GNSS) radio occultation (RO) [16,17].While ISR measurements provide detailed information on the sporadic-E layers, the limited number of global sites and measurement cadence inhibit statistical analyses of E s layers.Ionosondes provide direct measurements of sporadic-E intensity and height, but layer characteristics such as length or vertical thickness are difficult to obtain, given the spatial separation between ground sites and the vertical orientation of the measurements.
GNSS-RO has shown promise for monitoring [17][18][19] and characterizing [20][21][22] sporadic-E layers due to the relatively thin vertical structures that significantly perturb GNSS signals as they propagate horizontally across the ionosphere to a LEO satellite.This geometry allows for characterization of sporadic-E using GNSS-RO, while the relatively long horizontal lengths require exceptionally strong layers to observe in ground-based GNSS sensors [23].While the large number of global GNSS-RO measurements with high spatiotemporal resolutions (see [24] for more detail) provides a unique method for characterizing global E s characteristics and patterns, the amplitude and phase measurements are often difficult to interpret because of the integrated nature of the observations [25,26].
To help with the interpretation of GNSS-RO observations of sporadic-E, multiple-phase screen simulations of the signal perturbation have guided the mapping of sporadic-E characteristics to amplitude perturbations [17,20,27].As shown in the simulations, the magnitude of amplitude fluctuations in terms of a scintillation index such as S 4 are strong functions of the sporadic-E layer's vertical thickness, length, intensity (foEs), and orientation.While the vertical thickness can be estimated from GNSS-RO observations [20], the length and intensity are both unknown from GNSS-RO measurements alone, which makes the intensity estimates uncertain [25].However, simulations can help in our understanding of the complex mapping between the layer properties and observed GNSS-RO signals.
Here, we provide a statistical view of GNSS-RO scintillation observations in terms of the phase standard deviation σ φ and amplitude scintillation index S 4 over a month-long period over the mid-latitudes during the local summer at an E-layer altitude of 100 km.Then, to help interpret the results, we simulate GNSS-RO signals after traversing sporadic-E layers in an attempt to reproduce the measured data.Ultimately, we are able to reproduce many of the unique trends from the observations using standard E s distributions for length, intensity, and vertical thickness with the assumption that the E s layer acts as a Gaussian lens.Finally, inclusion of turbulence to the idealized Gaussian lens helps to improve the agreement between the simulation and measurements.

S 4 and σ φ Calculations
In this study, the S 4 index is calculated from the Constellation Observing System for Meteorology, Ionosphere, and Climate (COSMIC-1) 50-Hz GNSS-RO atmPhs amplitude data, namely the so-called signal-to-noise ratio (SNR), using the conventional definition where I = (SNR) 2 represents the signal power and the SNR is reported as a ratio in receiver voltage units (V/V) [28].denotes the running average using a truncation window defined by N cut from a time series, and thus the S 4 index is an N cut -dependent variable.In this study, we used three different truncation values (N cut = 51, 121, and 201) for the S 4 calculations, which corresponded to 2.2, 5, and 9 km in tangent height for the typical rate of rising and setting ROs, respectively.These three truncation windows were used to analyze various scale sizes in the scintillation, which helped to both interpret observations and compare them against simulations.
For the phase scintillation index σ φ , we computed the standard deviation of the detrended phase perturbation φ': where φ' = φ − φ is derived from the difference between the raw atmPhs excess phase measurement and its running mean.Because the φ varies exponentially with the height, φ needs to be run iteratively twice for each profile in order to obtain an unbiased φ'.
The definition for is same as in the S 4 calculation and determined by the truncation parameter N cut .
Sporadic-E (E s ) produces strong scintillations in the amplitude and phase of GNSS-RO signals, but the amplitudes of S 4 and σ φ are not necessarily proportional to each other.While an early study [17] compared the SNR and phase perturbations (φ'), here, we focus on the relationship between S 4 and σ φ exclusively for GNSS-RO L 1 measurements.Figure 1 shows two cases where S 4 appears to be capped at 1.5, while σ φ can reach up to ∼0.4 m.
The first case (top panel in Figure 1) was a typical and popular E s event where the S 4 and σ φ oscillations were confined mostly to altitudes of 80-110 km.This narrow height distribution suggests that the E s layers were generally thin and horizontally stratified.The RO limb sounding was particularly sensitive to these thin layers, causing strong scintillations when the RO tangent heights were parallel to the layers.However, the scintillation amplitude was significantly reduced at heights lower than the E s layer when RO passed through the layers at a slanted angle.A detailed discussion of the geometric smearing effect on the E s amplitude by RO can be found in [17].φ and σ φ represent the L 1 phase measurement and its scintillation, respectively.LST refers to Local Solar Time, Lat is the latitude, Lng is the longitude, and the Straight Line Ht is the straight line tangent altitude (SLTA).
The second case in Figure 1, although not occurring very frequently, was also an E s phenomenon.This event had an extended impact on the scintillation measurements below a ∼115 km tangent height where an E s occurred.The extended E s can happen when its horizontal extent is large (i.e., 1000 km or greater) and associated with spatial variations of many different scales [17].Strong gravity waves coupled with planetary wave activities, which can create embedded structures in the E s layer, are capable of producing this type of scintillation.Observational evidence (e.g., [29]) and a model illustration [3] all supported such extended impacts from E s .Alternatively, the extended scintillation could be due to equatorial F-region irregularities, which are known to produce the C-type scintillation [30] displayed here.

Joint S 4 -σ φ Histograms
To analyze the GNSS-RO scintillation statistics associated with E s , we focused on the mid-latitude measurements (−45 • latitude) over a month-long period during the local summer (January 2008) at an E-region altitude of 100 km.These parameters were selected to maximize the likelihood of sporadic-E causing scintillation, which is known to occur mostly in the mid-latitudes during the local summer [15,[17][18][19]] at altitudes near 100 km [20].While other ionospheric irregularities such as EPBs and polar cap electron density gradients are also known to cause significant satellite-based GNSS scintillation [6], our focus on mid-latitude measurements reduced the likelihood of these other ionospheric irregularities contributing to the observed scintillation (see [3] for more details).However, it did not completely remove the influence of equatorial and polar F-region irregularities.Furthermore, it is known that these ionospheric irregularities have a strong dependence on the local time and other geomagnetic and wind conditions [31][32][33][34][35], but here, we included all measurements taken during the month of interest to increase the total number of data points.By analyzing the scintillation at a single altitude instead of locating the maximum S 4 or σ φ for each occultation, we sampled various regions of the diffraction patterns produced by the E s layers with various altitudes centered around 100 km.This sampling was exactly what we attempted to simulate for comparison.
The S 4 saturation from E s was further quantified with a joint two-dimensional (2D) S 4 -σ φ histogram, which clearly showed limited growth by S 4 at large σ φ values.As shown in Figure 2a, where N cut = 51 was used to extract only thin E s layers, the S 4 saturation started to occur at σ φ ∼ 0.03 m.Although most of the S 4 and σ φ measurements were linearly correlated at small σ φ values, a significant fraction of E s had a very large σ φ value that was not linearly proportional to S 4 .The S 4 values seemed to have a mean of 0.8 as σ φ increased.
At the larger truncations (i.e., N cut = 121 and 201), where scintillations from thicker E s layers were included, the S 4 saturation was still evident, but it had a gradual transition from unsaturated to saturated.In these cases, the S 4 values tended to approach the same mean at 0.8, and the statistics showed very few samples with S 4 > 1, especially for N cut = 201.The reduction in S 4 for the larger windows was caused by the narrow E s vertical thicknesses relative to the window lengths, which reduced the S 4 values when scintillation did not occur over the entire window length.For σ φ , we saw an increase with increasing window lengths, which was due to larger-scale phase perturbations that were not captured using thin windows.More details on the phase and amplitude perturbations will be provided in Section 3.1.From this, there are benefits to using three different truncation windows for analyzing various intensity and phase perturbations.

Gaussian Lens Simulations
To model the scintillation observed in the RO data (Figure 2), the E s characteristics from previous studies were randomly sampled with the assumption that the layer acted as a perturbing Gaussian lens.An analytic solution for the electric field pattern following a Gaussian lens follows where E is the complex electric field, i = √ −1, X = x/r 0 is the dimensionless transverse coordinate, and Z = z/kr 2 0 is the dimensionless coordinate in the propagation direction [36].Here, r 0 is the 1/e thickness for the lens and k = 2π/λ is the wavenumber.The Gaussian lens in the transverse plane follows where φ is the phase contribution from the lens and φ 0 is the peak lens strength [36].For our simulations, we used an upper limit of p = 100 to calculate the complex electric field at a distance of z = 3000 km after the lens, which is characteristic of RO geometry and the distance traveled to the LEO satellite after reaching the tangent point.A transverse resolution of ∆x = 4λ L 1 was used for the simulations, where λ L 1 was the L 1 wavelength.
While the assumption of a Gaussian lens profile for a sporadic-E layer is certainly an oversimplification, it provides a method to quickly simulate a variety of E s lengths, thicknesses, and strengths to analyze the general trends.The complex, turbulent structures observed by incoherent scatter radars [37,38] will provide additional variability and will be explored in Section 4.
The sporadic-E parameters were derived from previous studies for the length, thickness, tilt, and intensity.The length follows a lognormal distribution with a mode of 170 km and σ = 0.7 [39].We also reduced the lengths by 65%, following the geometry argument in [27] for randomly slicing through sporadic-E layers with horizontal lengths of ∼10 times the horizontal thickness.
The vertical thicknesses follow a lognormal distribution with a mode of 1.5 km and σ = 0.4 [20].This vertical thickness corresponds to the thickness where the signal intensity reaches 95% of the background intensity within the typical u-shape profile observed for E s .As observed in Figures 3 and 4, the Gaussian lens is generally thinner than the ushaped intensity dip, especially for the moderate and stronger lenses.This increase in thickness for the intensity profiles is caused by the negative lens, which diverges the signal away from the lens at large distances.While the relationship between the 95% intensity definition of thickness from [20] and the Gaussian lens thickness defined for simulations is rather complex and requires more research to fully understand, here, we define the vertical thickness as the altitude range between the points where the phase reaches 20% of the lens peak.This 20% value was selected instead of the full width at 1/e to reduce the lens thicknesses and match the observed S 4 -σ φ trends in Figure 2. Additional research is required to properly map the thicknesses derived in [20] to Gaussian lens thicknesses for simulation or observation.However, this simple mapping provides a starting point for simulation of GNSS-RO statistical observations.
Sporadic-E intensities in the form of foEs were taken from [22], assuming a Gaussian profile with a mean of ∼3 MHz and a standard deviation of ∼1 MHz (estimated from Figures 8 and 9 in [22]).The background E-layer density should be taken into account to properly capture the lens strength magnitudes through conversion of the foEs to an foµEs, as outlined in [40], where foµEs refers to the metal plasma critical frequency of the E s layer.However, we did not explicitly convert the foEs to an foµEs here, as we found that the foEs intensity distribution from [22] provided a decent match to the observed S 4 -σ φ scintillation trends.Furthermore, as foµEs calculations require knowledge of the local E-layer densities at the time and location of measurement, it is difficult to calculate an foµEs distribution directly from an foEs distribution without additional information.We hope that the statistical approach outlined here helps provide insight into foµEs distributions by providing a method for adjusting the distribution to improve agreement between simulations and observations.The E s parameters were converted to effective Gaussian lenses using where n p is the index of refraction at the center of the E s layer, l is the horizontal length of the layer, f L1 is the GPS L 1 frequency of 1575.42MHz, and k is the L 1 wavenumber.With n p < 1, φ 0 < 0, which creates a negative lens.Here, we ignored the background ionosphere to simplify the simulations, and the E s tilts were not included to focus on the moderate-to-strong scintillation trends.As shown in [20], layer tilts significantly reduced the amplitude scintillation at the LEO measurement plane, which introduced additional data points near the origin of our S 4 -σ φ histograms.As these points did not significantly contribute to the moderate-to-strong scintillation trends observed in Figure 2, they were ignored for this analysis.

Simulation Examples
Examples of the Gaussian lens simulations are displayed in Figures 3 and 4 for the phase and intensity of the GPS L 1 signal, respectively, as measured at the LEO measurement plane after propagating through a sporadic-E layer with a vertical thickness of 1.5 km.Three varying lens strengths are displayed to show the different scintillation responses, using a 2.2 km window (N cut = 51) to calculate S 4 and σ φ .The lens strengths were −1 rad (−3 cm), −5 rad (−15 cm), and −10 rad (−30 cm), showing weak, moderate, and strong lensing, respectively.For perspective, the −5 rad lens strength may correspond to an E s layer with a 3.4 MHz foEs and a length of 65 km.As observed in Figure 3, the measured signal phase closely tracked the Gaussian lens for the weak and moderate lenses, but significant phase scintillation reduced the measured phase for the strong lens.Interestingly, while the scintillation on the edges of the lens was most severe for the strong lens, the largest σ φ was observed for the moderate lens scenario with a peak of ∼6 cm.This counter-intuitive result was due to the large phase gradient from the lens that was mostly captured by the moderate lens, but it was less severe for the strong lens because of the phase scintillation observed at the edges, which reduced the peak phase measured in the lens.These results indicate that thin, strong lenses can produce relatively low peak σ φ values, but the phase scintillation at the edges of the lenses will be enhanced.The signal intensities at the LEO measurement plane for the same lenses are displayed in Figure 4. Here, the characteristic U-shaped profile described in [20] is present for the moderate and strong lenses, where the low intensity region at the bottom of the U shape was caused by the negative lensing effect of the E s layers.Additionally, the ringing around the u-shape fade is due to diffraction [20].However, the weak lens did not produce a strong diffraction pattern, and only the beginning of a U-shaped pattern was present.S 4 showed elevated peaks of 1.3 for the moderate lens and 1.5 for the strong lens, with a relatively low peak of approximately 0.3 for the weak lens.In the center of the moderate and strong profiles, the S 4 magnitude was elevated by the decreased intensity mean at the center of the U-shaped profile.For the S 4 techniques calculated using a longer running average of the background intensity, such as those provided by the COSMIC Data Analysis and Archive Center (CDAAC) [28], the S 4 calculated over one-second intervals (50 points for 50 Hz data) would be nearly equal to the standard deviation of the normalized intensity.Following the results from [27], the peak S 4 plateaued for the moderate-to-strong lenses, and the roughly linear peak S 4 relationship with the lens strength was mostly observed for the weak-to-moderate lenses.Additionally, the vertical lens thickness played a large role in the phase and amplitude scintillation, with thinner lenses generally causing more scintillation in both phase and amplitude.However, as discussed previously, σ φ also depends on the large phase perturbation induced directly by the E s layer and not strictly the scintillation from diffraction.These examples help to interpret the statistical results outlined below.

Simulated Statistics
In an attempt to reproduce the observations displayed in Figure 2, the E s length, vertical thickness, and foEs distributions were randomly sampled 10,000 times, and the resulting electric field profiles in the LEO measurement plane (z = 3000 km) were stored for each of the 10,000 simulations.Histograms of the sampling along with the distributions are displayed in Figure 5a for the length, Figure 5b for the vertical thickness, and Figure 5c for the foEs.[39], the vertical thicknesses from the lognormal distribution outlined in [20], and the foEs magnitudes from the Gaussian distribution outlined in [22], along with histograms of the 10,000 random samples used in the simulations.
Joint S 4 -σ φ histograms of the simulated S 4 and σ φ at the LEO measurement plane are displayed in Figure 6.Each grid point in the measurement plane was used as a single data point to create the histogram.Extremely low values (below 0.001) were removed for both S 4 and σ φ to improve readability.The majority of points occurred in the region of S 4 < 0.3 and σ φ < 0.02 m, corresponding to weak intensity and phase fluctuations, respectively.This region corresponds to the weak lens simulations in Figures 3 and 4.
Outside of this weak lens region, there were three regions of interest.The first most noticeable region was caused by significant phase scintillation with elevated S 4 and low σ φ .When strong phase scintillation occurred, the peak phase perturbation magnitude was reduced, which in turn reduced the σ φ peak, as displayed in Figure 3.In this region, both phase and amplitude scintillation are important, and techniques such as those outlined in [17,19,22] will work well, while techniques relying on phase perturbation alone will struggle with properly characterizing the layer.With the exception of the 2.2 km window, the simulations showed significant counts in this region while the observations (Figure 2) did not.This indicates that the simulations were oversampling the strong-thin regions, and a method to reduce this sampling is presented in Section 5.The second region of interest was produced by the moderate lenses, resulting in a roughly linear increase in σ φ with respect to S 4 .For the 2.2 km window, this moderate lens zone was much weaker than the strong scintillation zone, but the strong scintillation linear trend from the origin to S 4 ≈ 0.7, σ φ ≈ 0.02 m roughly matched the observations in Figure 2. The moderate lens region was more obvious for the 5.0 km and 9.0 km windows, with a ∆σ φ /∆S 4 slope of 0.05/0.6m for both windows.These slopes were in line with the observed slopes from the origin to σ φ ≈ 0.05 m, although the observations showed a slightly elevated slope of 0.05/0.5 m for the 9.0 km window.
However, the simulated moderate lens region extended throughout S 4 ≈ 1.2, which was not observed in the GNSS-RO measurements.Instead, Figure 2 shows the S 4 plateau between 0.6 and 0.8.This difference was likely due to the assumption of a Gaussian lens for the E s layers, which was certainly an oversimplification, as shown by the incoherent scatter radar measurements (e.g., [37,38]).The realistic E s layers were patchy and billowy, which was caused by the Kelvin-Helmholtz instability [41], reducing the ideal strong lensing and also introducing additional perturbations and scintillation caused by the billows.The Gaussian lens's focusing diffraction pattern is capable of producing very large peak intensity scintillation, while the billowy layers from turbulence are not capable of producing the same focusing.This important point will be explored further in Section 4.
Finally, the third zone of interest was characterized by elevated σ φ and low S 4 (Figure 6a) due to thicker layers with relatively weak lens strengths.These thicker, weaker layers caused phase perturbations without causing significant diffraction in the signal intensity.The GPS-RO techniques relying on phase perturbation (e.g., [18,25]) will work well in this region, while the techniques relying on the signal intensity (e.g., [19,22]) will have difficulty characterizing these thicker, weaker layers.
While the simulations reported here were strictly for GNSS-RO geometries, it must be noted that many previous studies have compared and validated RO techniques for monitoring sporadic-E with ground-based ionosondes [19,21,25,26].The foEs distribution [22] was derived from RO observations with an empirical fit to ionosonde measurements, and the length distribution was taken from a topside sounder (Explorer 20) [39].There is certainly room for improvement in GNSS-RO techniques for monitoring sporadic-E, but the general agreement with ground-based sensors provides confidence in the results obtained only from GNSS-RO observations.

Simulations Including Turbulence
Following [42], the spectral density function Φ δn (k) for random variations in the index of refraction δn( r) can be described by where δn 2 is the refractive index variance and Q(k) is the Shkarofsky [43] spectral density function: Here, k S corresponds to the inner scale (smallest structure), k L corresponds to the outer scale (largest structure), K ν is the modified Bessel function, and ν is the power law index parameter.Additionally, we have where C s is the turbulent strength parameter.
To include a turbulence spectrum into our simulations, we used the following conditions: k S = 2π/100 rad/m, k L = 2π/∆x rad/m where ∆x is the vertical thickness of the layer, ν = 4/3 for the Kolmogorov turbulence, and C s = 10 −8 .A Kolmogorov power law was selected based on the results in [44], and the value of C s was adjusted to improve the simulations agreement with the observations displayed in Figure 2 while also producing vertical refractivity profiles with variations similar to those observed for the horizontally integrated E s profiles from the ISR measurements as displayed in [37].This turbulence was added to the initial Gaussian lens by generating complex frequency samples of the Shkarofsky spectrum using a normal distribution with a standard deviation of one and a mean of zero for both the real and imaginary components.The Fourier transform of these spectral random realizations provided random fluctuations in space that were added to the original Gaussian lens to provide a turbulent sporadic-E layer: where n turb (x, z) is the two-dimensional index of refraction with the inclusion of turbulence, n 0 (x, z) is the original refractive index profile, and δn(x) is the random index of refraction variation in the vertical (transverse) direction.Here, we used the same vertical δn(x) for all z in the horizontal (propagation) direction to focus on vertical lens irregularities.
Example phase and amplitude profiles for turbulent lenses are displayed in Figures 7 and 8.The same lens strengths of −1, −5, and −10 radians used for the Gaussian lenses in Figures 3 and 4 are analyzed here.However, the vertical thickness was increased to 2.5 km to clearly see the turbulence.The phase profiles showed similar behavior to the Gaussian lenses, with the exception of reduced phase scintillation for the −10 radian lens due to the thicker vertical thickness that reduced the overall vertical gradient in the index of refraction.A slight asymmetry was observed for the −10 radian lens, with a stronger asymmetry observed for the −5 rad lens.While the phase profiles were similar to the Gaussian lenses, the amplitude profiles were significantly perturbed (Figure 8).This was most obvious for the −5 and −10 radian lenses, which showed large asymmetries in the idealized U-shaped profiles.This asymmetry was a result of the asymmetric lens caused by the inclusion of turbulence.The asymmetric U-shaped profile presented here is commonly measured in GNSS-RO observations (e.g., Figure 1), resulting from the small scale structures observed in E s layers.Additionally, the strong S 4 caused by lens focusing was reduced in the turbulent lenses because of the lens imperfections.This reduced focusing ultimately resulted in the S 4 plateau observed in Figure 1.It must be noted here that many GNSS-RO observations of sporadic-E show multiple layers [45].These multi-layered structures result in either multiple U shapes or a blend caused by the interference of the two diffraction patterns.Examples of these multi-layered profiles can be found in previous studies (e.g., Figure 3 of [20]) and throughout the CDAAC catalog (e.g., atmPhs_C001.2010.029.01.35.G14_2010.2640_ncand atmPhs_C002.2008.001.00.28.G04_2013.3520_nc).The simplified modeling approach taken here did not account for these multi-layered structures, as we were focused on the scintillation produced by single layers.While the multi-layered events with sufficient altitude spread would not change the scintillation statistics of the individual layers, the interfering diffraction patterns were not accounted for in this study.After the turbulence was added to the Gaussian lens, an analytic solution to the electric field (i.e., Equation (3) for the Gaussian lens) no longer existed, and the electric field at the GNSS-RO measurement plane had to be calculated using a numerical technique.Here, we used the multiple-phase screen approach due to the long propagation distance of 3000 km after the perturbing layer.This multiple-phase screen approach is commonly used for RO propagation through sporadic-E layers [17,20,27], and more information on the technique can be found in [42,46].Specifically, we used a transverse (vertical) grid spacing of 76 cm (four L 1 wavelengths) for a total transverse length of 50 km.In the propagation (horizontal) direction, we took 1000 equal steps through the E s layer followed by ∼20 km steps after the layer for a total propagation distance of 3000 km.
Using the same distributions implemented for Figure 6 with the addition of turbulence provided the results displayed in Figure 9.The linear slopes were similar to the slopes observed for the Gaussian lens, but plateauing of S 4 ≈ 1.0 better matched the observations (Figure 2).This was due to a reduction in the Gaussian lenses focusing with the lens "imperfections" added by the turbulence.Additionally, the 9.0 km window ∆σ φ /∆S 4 slope from the origin was increased to 0.05/0.5 m, which matched the observed slope better than the Gaussian lens simulations.The similarity of these simulations to the observations supports previous sporadic-E measurements that showed the turbulent nature of these layers [37,38].Including a relatively small degree of turbulence in our simulations provided more realistic sporadic-E layers and improved our agreement with the GNSS-RO observations.While the inclusion of turbulence helped to reduce the high σ φ and S 4 region caused by focusing, the strong phase scintillation region with strong S 4 and low σ φ remained more enhanced than the observations suggest.Similar to the Gaussian lens results, the 2.2 km window showed reasonable agreement with the observed ∆σ φ /∆S 4 slopes from the origin, but the 5.0 and 9.0 windows did not show this enhancement.This implies that we were oversampling the strong, thin layers with large vertical gradients in the index of refraction, and some physical limitation on the coupled E s characteristics may have existed that was unaccounted for in our uncoupled random sampling.

Diffusion-Limited
The region of elevated S 4 and low σ φ observed in our simulations (Figures 6 and 9) was not present in the observations displayed in Figure 2, indicating that the strong, thin layers were oversampled from the length, strength, and vertical distributions.To limit the number of strong, thin lenses in an attempt to improve agreement with the observations, we removed the randomly sampled layers with φ 0 /r 2 0 > 13.5 rad/km 2 .This threshold value was obtained by averaging the ratios obtained for the φ 0 values, where various vertical thicknesses in the range of 0.5 ≤ r 0 ≤ 1.5 km began to show significant phase scintillation.Here, we define significant phase scintillation as a peak measured phase difference of greater than 30% from the applied phase, illustrated by the black and red lines in Figure 3. Physically, this limit on φ 0 /r 2 0 may be related to a diffusion limited layer as D i d 2 n i dx 2 ∝ φ 0 /r 2 0 at the center of the layer, where D i is the diffusion coefficient for metal ions, n i is the metal ion density, and x is in the vertical direction.Ion convergence caused by wind shears is counteracted by this vertical diffusion such that the final layer characteristics will be dependent on both processes [8].Layers that become too strong and thin will experience elevated diffusion that will ultimately place a limit on the allowed φ 0 /r 2 0 ratio.Using this diffusion limit on the randomly sampled distributions resulted in the remaining parameters displayed in Figure 10.Out of the 10,000 initial samples, 4383 layers with φ 0 /r 2 0 > 13.5 rad/km 2 were removed.Compared with Figure 5, the layers with elevated values for length or foEs (φ 0 is proportional to the length-foEs 2 product) and low r 0 were removed from the distributions.The simulated S 4 -σ φ histograms for these diffusion-limited layers are displayed in Figure 11.Compared with the original distributions used to create Figure 6, the strong phase scintillation region with elevated S 4 and low σ φ was greatly reduced.Similarly, the S 4 -σ φ histograms for the turbulent lenses were also improved by the removal of layers with φ 0 /r 2 0 > 13.5 rad/km 2 (Figure 12).While the strong phase scintillation region was not completely removed in the simulations, inclusion of the diffusion-limited threshold on the φ 0 /r 2 0 ratio significantly improved agreement with the observations.The threshold enforced here was selected to show a proof of concept, but placing physical thresholds on the randomly sampled distributions should be enforced to ensure that unrealistic sporadic-E layers are not produced.Determination of these thresholds requires insight into the coupling between each of the parameters, which we treated as independent up to this point.Alternatively, the reduction of longer, thinner, and stronger layers may indicate that the E s distributions used here require adjustment to match the observations.For the length distribution, the authors of [23] found that strong layers have a median length of 100 km, which is a reduction from the 170 km mode used here from [39].Additionally, taking the background E-layer into account through a conversion of foEs to foµEs as outlined in [40] would help to reduce the effective lens strengths, since φ 0 ∝ length-foEs 2 .Certainly, additional research is warranted to help improve the E s distributions and to determine the coupling characteristics.

Conclusions
A statistical view of GNSS-RO phase and amplitude scintillation at an E-region altitude in the mid-latitudes during the local summer indicates that sporadic-E is likely the primary cause of the unique S 4 -σ φ trends.The GNSS-RO observations showed a nearly linear trend in joint two-dimensional S 4 -σ φ histograms until the S 4 saturated near 0.8.To understand the trends, 10,000 RO simulations through sporadic-E layers were performed using length, intensity, and vertical thickness distributions from previous research with the assumption that the E s layer acted as a Gaussian lens.This simplified view of the layers provides insight into different scintillation statistics for weak, moderate, and strong lenses.
Using the Gaussian lens assumption, a roughly linear ∆σ φ /∆S 4 slope from the origin to σ φ ≈ 0.05 m was observed in both the observations and simulations.In the simulations, this trend was produced strictly by E s layers, which implies that the observed GNSS-RO scintillation was also likely caused by E s .While the Gaussian layer simulations showed agreement with the primary trends, a region of elevated σ φ and S 4 was present in the simulations but not the observations.This elevated S 4 -σ φ was caused by focusing from these idealized lenses.In reality, the E s layers were turbulent and billowy.Therefore, we included Kolmogorov turbulence by adding random realizations of a Shkarofsky spectral density function to perturb the idealized Gaussian lenses.The addition of these random turbulence fluctuations made the E s layers more physically realistic, which in turn improved agreement with the observations by removing the large number of points in the elevated σ φ and S 4 region.Additionally, a plateau of S 4 ≈ 1.0 was produced by these turbulent layers, which better matched the observed trends.
Simulations using the original distributions without constraints on the coupled layer parameters showed large counts with significant phase scintillation (elevated S 4 and low σ φ ) produced by strong, thin lenses.This strong phase scintillation region was not present in the observations, which indicated an oversampling of these strong, thin layers.To reduce this oversampling, a threshold was placed on the φ 0 /r 2 0 ratio which was proportional to the ion diffusion rates.After placing a limit on this ratio, counts in the elevated S 4 and low σ φ region were significantly reduced, resulting in better agreement with the observations and implying that there may be a diffusion limit on E s layers.The improvement was best illustrated by the 2.2 km window simulations, likely due to the narrow window being more sensitive to small-scale phase fluctuations.
This research provides a new method for interpreting GNSS-RO results in bulk through large numbers of simulations guided by predefined probability distribution functions.For future work, the distributions can be adjusted to find the best fits to the statistical observations and improve our understanding of the sporadic-E layer characteristics.However, it must be noted that certain ambiguities arise with respect to the length × foEs 2 product contributing to the strength of the lens, and the diffraction pattern is also highly dependent on the lens thickness, which requires physically realistic bounds and starting points placed on the distributions using previous E s measurements.Improving the turbulence estimates and implementation can also help improve agreement and provide insight into the complex layer structure.Finally, it may be possible to apply this approach to equatorial and polar observations to enhance characterization of the ionospheric disturbances in these regions.While the scale sizes, general characteristics, and challenges of equatorial and polar observations are markedly different than those for mid-latitude E s , a statistical scintillation approach taking advantage of RO geometry provides a complementary dataset to the traditional ground-based GNSS measurements.

Figure 2 .
Figure 2. Joint S 4 -σ φ histograms: (a), (b), and (c) for truncation N cut = 51, 121 and 201, respectively, over all times during January 2008 between −40 • and −50 • latitude and 100 km altitude where E s is maximized.The symbol line in each panel highlights the S 4 peak value at each σ φ bin.The S 4 and σ φ bin sizes are 0.02 and 2 mm, respectively, and the white bins indicate no data.All 2D histograms are normalized by their peak values, and a total of ∼1.1 million samples were used to create the histograms.

Figure 3 .
Figure 3.Example simulations for Gaussian lenses with vertical thickness of 1.5 km and lens strengths of −1, −5, and −10 rad.The Gaussian lens strength is displayed along with the simulated signal phase and σ φ at the LEO measurement plane.Note that the different lenses were shifted by 10 km in altitude for readability, and a window of 2.2 km was used for the σ φ calculations.

Figure 4 .
Figure 4. Simulated signal intensity and S 4 at the LEO measurement plane for the Gaussian lenses displayed in Figure 3.Note the elevated S 4 near the center caused by the reduced mean at the bottom of the U-shaped profile.The lenses were shifted by 10 km in altitude for readability, and a window of 2.2 km was used for S 4 calculations.

Figure 5 .
Figure5.Probability distribution functions (PDFs) of the E s lengths from the lognormal distribution outlined in[39], the vertical thicknesses from the lognormal distribution outlined in[20], and the foEs magnitudes from the Gaussian distribution outlined in[22], along with histograms of the 10,000 random samples used in the simulations.

Figure 6 .
Figure 6.Simulated S 4 -σ φ histograms, assuming the E s layer is a Gaussian lens with 10,000 randomly sampled parameters using distributions obtained from previous E s observations.The histograms are normalized by their peak values.White bins indicate no data, and the symbol line in each panel highlights the S 4 peak value at each σ φ bin.

Figure 7 .
Figure 7. Same as Figure 3 but with the inclusion of turbulence.Additionally, a 2.5 km vertical thickness was used to clearly show the turbulent layers.

Figure 8 .
Figure 8. Same as Figure 4 but with the inclusion of turbulence.Additionally, a 2.5 km vertical thickness was used to clearly show the turbulent layers.

Figure 9 .
Figure 9. Simulated S 4 -σ φ histograms using 10,000 randomly sampled E s parameters with the addition of turbulence to the Gaussian layer.