© D. Vasylyev et al., Published by EDP Sciences 2024

1 Introduction

The trans-ionospheric electromagnetic wave experiences random modulation of its amplitude and phase as it propagates through the inhomogeneous ionosphere. The strength of such modulation depends on several factors such as time of the day and season, period of high or low solar activity, geomagnetic conditions, geographical location, propagation link geometry, etc. As a result, the morphology of scintillation events is complex and varied (Aarons, 1982).

Recent morphological studies have shown that there is a certain correlation between the occurrence of strong spatial gradients in the ionospheric electron density and intense phase scintillations. The evidence for such a correlation has been reported mainly in polar and high latitude regions, e.g., by studying scintillation at the boundaries of the ionospheric trough (Vo & Foster, 2001), storm enhanced density (SED) regions (Doherty et al., 2004; Coster & Skone, 2009; Sun et al., 2013), the polar tongue of ionization (TOI) (Kinrade et al., 2012; van der Meeren et al., 2014), at the edges of polar cap patches (Jin et al., 2015; Lamarche et al., 2022), or in cusp flow channels (Spicher et al., 2020). The onset of intense scintillation associated with the occurrence of strong total electron content (TEC) gradients in equatorial regions, e.g., due to scattering at the edges of equatorial plasma bubbles (EPBs), has also been reported (Vijayakumar and Pasricha, 1997; Kintner et al., 2007; Olwendo et al., 2012; McNamara et al., 2013; Sato et al., 2021). The signatures of gradient-associated scintillation can also be observed at mid-latitudes (Vadakke Veettil et al., 2017; Fallows et al., 2020; Lee et al., 2021).

Such gradient-associated scintillation is important to account for, e.g., relative GNSS positioning, which is affected not only by the absolute TEC of the ionosphere, but often even more by horizontal gradients in the electron distribution (Wanninger, 1993; Strangeways, 2000). The ionospheric electron density gradients also result in the azimuth shifts on SAR interferograms that cause the target images to be displaced from the true azimuth location (Zhu et al., 2016). In some cases, the anomalous ionospheric spatial gradients caused by geomagnetic storms may escape the detection and lead to integrity failures in Ground-Based Augmentation Systems of airports (Luo et al., 2004; Stankov et al., 2009). The dual-frequency positioning techniques can be useful for mitigation of ionospheric temporal gradients; however, the traces of large spatial gradients may still remain (Konno et al., 2006) and pose a threat to various positioning, navigation, timing, and reconnaissance services.

Scintillation indices have long been used to characterize the ionospheric perturbation of radio waves (Yeh & Liu, 1982). Additionally to this, the values of these indices have been used to diagnose the level of the ionospheric disturbance. Because of the mentioned correlation of scintillation indices with the presence of strong electron density gradients, the alternative ionospheric disturbance indices based on gradients have been proposed to fit certain research problems (Borries et al., 2020). The most notable examples of such gradient-based disturbance indices are: the Rate of change of TEC (ROT) (Wanninger, 1993; Aarons et al., 1996), the spatial gradient of TEC (Carrano & Groves, 2007; Cesaroni et al., 2015), the Rate Of TEC Index (ROTI) (Pi et al., 1997), the Along Arc TEC Rate (AATR) (Juan et al., 2018), the Disturbance Ionosphere Index (DIX) (Jakowski et al. (2012), the DIXSG index (Wilken et al., 2018), the Rate Of change of electron Density Index (RODI) (Jin et al., 2019; Wood et al., 2022), just to name a few. On one hand, the advantage of some of these indices, such as the ROTI, is that their values can be obtained without the need for specialized and costly scintillation detectors. The data acquisition rate for the gradient indices is also lower than for scintillation indices, which is an advantage when the data storage capacities are limited. On the other hand, a sufficient linear correlation between the scintillation indices and the ROTI has been reported by Basu et al. (1999), Yang & Liu (2016), Luo et al. (2018), Acharya & Majumdar (2019), De Franceschi et al. (2019), which means that the ROTI can be used for effective scintillation monitoring. Other indices, such as the RODI, are constructed from in-situ measurements, correlate with scintillation levels observed on the ground, and thus can provide good insight into scintillation-related ionospheric irregularities (Olwendo et al., 2019; Kotova et al., 2023).

One can ask whether the scintillation theory can explain the correlation of the scintillation indices and the ionospheric gradients. The standard theoretical approach relates the random phase contributions, gained by a radio wave propagating through a randomly inhomogeneous ionosphere, to the fluctuations of the TEC along the propagation path. A ground-based receiver operating in the Fresnel diffraction zone records just these original phase variations. At the same time, amplitude scintillation would additionally experience the so-called Fresnel filtering effect, i.e., ionospheric irregularities of scales larger than the radius of the first Fresnel diffraction zone would not contribute to the recorded amplitude fluctuations (see e.g., Bhattacharyya et al., 2000). Both the phase and amplitude scintillation indices are then related to the variance of the fluctuations of the TEC along the wave propagation path.

The relationship between the scintillation indices and the gradient-based ionospheric disturbance indices is not obvious. For example, such a relationship would suggest that the random phase increment gained by propagation through the ionosphere would depend on the spatial or temporal gradients of the fluctuating part of the TEC¹. In fact, some early theoretical investigations recognized that the phase fluctuation obtained by transmission through a random medium should also depend on such gradients (Chernov, 1960; Salpeter, 1966; Tatarskii, 1971). Such a dependence was attributed to random refraction at medium inhomogeneities, which introduces the excess ray path lengthening, i.e., the additional phase delay. However, under the assumption that the gradients of the refractive index are small, such ray path lengthening was neglected in the mentioned works. To the best of our knowledge, the gradient-dependent correction term is preserved in the complex phase method (Gherm et al., 2005) as well as in the recently proposed method by Hamza et al. (2023). An interesting result is also presented in Bhattacharyya et al. (2000), where the authors showed the relationship between the amplitude scintillation index and the second-order derivative of the standard deviation of variations in TEC (DROTI index) using the transport-of-intensity equation.

In this paper, the relationship between refractive scintillation and ionospheric disturbance indices is established within the framework of geometrical optics. Retaining the term of the second order of smallness of the medium refractive index fluctuations, we express the random phase increment, gained by the signal wave while propagating in the randomly inhomogeneous ionosphere, in terms of the fluctuating part of the TEC and gradients of the refractive index fluctuations. We show that for very strong gradients, the second-order corrections become significant and the correlation between the scintillation index and the gradient-based disturbance indices can be established.

In our investigations, we focus primarily on refractive scintillation due to scattering on medium-scale irregular structures at high latitudes, such as those found at the edges of polar TOI, SED, or at the boundary of the auroral oval. More specifically, we consider scattering from ionospheric irregularities with scales larger than the typical Fresnel scale for the specified signal frequency, e.g., larger than 300 m for L-band GNSS signals. The associated scintillation events are related to the phenomenon of “phase-without-amplitude” scintillation (Fremouw et al., 1978; Li et al., 2010; Forte et al., 2017; Hong et al., 2020; Nishimura et al., 2023), where intense phase fluctuations are observed, but little increase in amplitude fluctuations². This effect can be explained by the inhibition of the large to small-scale energy cascade in the process of formation of ionospheric irregularities caused by high recombination rates and/or strong conductivity along the field lines (Forte et al., 2017). The large-scale structures are then able to distort the wavefront of the signal wave while being larger than the radius of the first Fresnel zone, they do not contribute to amplitude fluctuations. The primary wavefront distortion, i.e., its tilt, is associated with random refraction at interfaces of media characterized by different refractive indices. The tilt grows with the growth of the refractive index gradient, the effect that is accounted for within the formalism presented in this article. We also show that the proposed method can be effectively applied to the description of amplitude scintillation caused by scattering on equatorial plasma bubbles in the low-latitude region. In this case, under the condition of well-developed and highly structured irregularities at the borders of plasma bubbles, both phase and amplitude scintillation have refractive contributions (Vats et al., 1981), and the “phase-without-amplitude” scintillation effect is absent.

For the simulation of refractive scintillation we extend the conventional method of random phase screen method to include the dependence on the refractive index gradients. This extension is partly similar to the phase gradient method presented by Schmidt and Miller (2019), which was primarily developed in order to study wavefront gradient sensors in adaptive optics applications. These authors used the gradient of turbulent screens to obtain the expansion coefficients of the detected wavefront in circular orthogonal modes. In this work, the originally generated phase screen together with its gradient values are used for radio wave propagation simulations within the aforementioned second-order geometric optics approximation. One parameter of the current approach, namely the signal ray deflection vector, is related to the novel data product NeGIX (spatial electron density gradient index) of the Swarm satellite mission. This allows us to integrate the in-situ measurements into a phase screen model to compare the results of numerical simulations with detected scintillation data.

The organization of this paper is as follows. Section 2 summarizes the conventional method of random phase screen simulation of transionospheric wave propagation. Section 3 presents the theoretical and algorithmic extension of this method to phase gradient screens. Section 4 applies the developed formalism to the simulation of phase scintillation at high latitudes and amplitude scintillation at low latitudes. Finally, conclusions and outlook are presented in Section 5. The main theoretical background is based on the method of Tatarskii (1971) and is summarized in Appendix A and Appendix B for completeness and referencing convenience. In Appendix C we show that the refractive tilt error originates primarily from the scattering on topside ionospheric irregularities. Appendix D gives a short overview of the electron density gradient index NeGIX derived from the Swarm electron density data.

2 Conventional random phase screen method

One of the popular models for simulating wave propagation through a random medium is based on the concept of a random phase screen. This phase screen is placed along the wave propagation path and randomly modulates the wavefront of the signal wave. The random modulation process is governed by the fundamental statistical properties of the random medium. The electromagnetic wave transmitted through the random phase screen forms a complex diffraction pattern at the observation plane located at some distance from the screen, while the wavefront arriving at the observation plane is randomly corrugated. The degree of distortion of the received signal amplitude and phase characteristics from the expected undisturbed ones is characterized by amplitude and phase scintillation indices, respectively.

2.1 Theoretical formulation

For a simple theoretical formulation of the problem we choose the coordinate system with the center located at the phase screen and with the z-axis aligned along the propagation direction. The receiver is located at a distance z = L from the phase screen. The phase screen imposes a phase perturbation

$$\delta \phi \left({\mathit{r}}_{\perp }\right)=k\underset{0}{\overset{s}{\int }}\mathrm{}\delta n({\mathit{r}}_{\perp },z^{\text{'}})\mathrm{d}{z}^{\text{'}},$$(1)

where k = |k| is the wave number corresponding to the wave vector k and δn is the random variation of the medium refractive index with respect to its mean. The integration is performed along the whole propagation path within the layer of the randomly inhomogeneous medium of slant thickness s. In the case of the ionospheric random media equation (1) can be rewritten in terms of the fluctuations $\delta N_{\mathrm{TEC}}={\int }_0^s\mathrm{}\delta N_e\mathrm{d}z$ of the TEC with respect to its mean value as

$$\delta \phi \left({\mathit{r}}_{\perp }\right)=-\lambda r_e\delta N_{\mathrm{TEC}}\left({\mathit{r}}_{\perp }\right),$$(2)

where δN_e is the random component of the ionospheric electron density, λ = k/2π is the signal wavelength and r_e is the classical electron radius (= 2.8 × 10⁻¹⁵ m). Equation (2) is derived by substituting the cold plasma approximation of the Appleton-Hartree dispersion relation for plasma and neglecting the terms of order $\delta N_e^2$ and higher (Rino, 2011; Vasylyev et al., 2022).

If the signal wave is spatially bounded in the transverse to the propagation direction, the complex field amplitude at the receiver location can be written in the paraxial approximation of the Fresnel diffraction formula (Mercier, 1962)

$$u({\mathit{r}}_{\perp },L)=\frac{\mathrm{ik}}{2\pi L}\int \mathrm{}u({\mathit{r}}_{\perp }^{\text{'}},0)\exp \left[-i\frac{k}{2L}|{\mathit{r}}_{\perp }-{\mathit{r}}_{\perp }^{\text{'}}{|}^2\right]{\mathrm{d}}^2{\mathit{r}}_{\perp }^{\text{'}},$$(3)

$$u({\mathit{r}}_{\perp },0)=u_0\left({\mathit{r}}_{\perp }\right)e^{-\mathrm{i\delta }\phi \left({\mathit{r}}_{\perp }\right)},$$(4)

where u₀ is the amplitude of the wave hitting the random phase screen. In the Fourier domain equation (3) can be written as

$$\tilde{u}({\mathit{\kappa }}_{\perp },L)=\tilde{u}({\mathit{\kappa }}_{\perp },0)\exp \left[\mathrm{i\lambda }L {\kappa }_{\perp }^2/4\pi \right],$$(5)

where $\tilde{u}$ is the Fourier image of u, and the Fourier transform is performed only with respect to the transverse coordinates. For weak scintillation conditions, the phase departure from the expected value at the receiver can be determined within the model of a single phase screen from equation (3) as (Yeh & Liu, 1982)

$$S({\mathit{r}}_{\perp },L)=\arg \left[u\right({\mathit{r}}_{\perp },L\left)\right]-\arg [\langle u({\mathit{r}}_{\perp },L)\rangle ]=\frac{k}{2\pi L}\int \mathrm{}\delta \phi \left({\mathit{r}}_{\perp }^{\text{'}}\right)\sin \left(\frac{k}{2L}|{\mathit{r}}_{\perp }-{\mathit{r}}_{\perp }^{\text{'}}{|}^2\right){\mathrm{d}}^2{\mathit{r}}_{\perp }^{\text{'}},$$(6)

where the averaging is done over possible realizations of the random process δφ.

It can be seen from equation (5) that for spatial frequencies ${\kappa }_{\perp }<<1/\sqrt{\lambda L}$ the exponential multiplier can be omitted. This means that if u(r_⊥, 0) has variations on a scale much larger than the Fresnel scale size $\sqrt{\lambda L}$, the field distribution at the receiver site would be directly related to the TEC variations along (slightly different) propagation paths ending at the receiver. At the same time, within the framework of the phase screen model, the variations at such large scales would not contribute to amplitude fluctuations and would be associated only with phase fluctuations. The resulting phase scintillation is due to random refraction. For shorter length scales, i.e., for ${\kappa }_{\perp }\gtrsim 1/\sqrt{\lambda L}$, the exponential term in equation (5) would generally contribute to the modulation of both amplitude and phase. The resulting field u(r_⊥, L) would exhibit both amplitude and phase scintillations, which are primarily associated with diffractive effects.

The levels of amplitude and phase disturbance of the transmitted signal are conventionally given by the corresponding scintillation indices:

$$S_4=\sqrt{\langle |u{|}^4\rangle /\langle |u{|}^2{\rangle }^2-1},$$(7)

and

$${\sigma }_S=\sqrt{\langle S^2\rangle -\langle S{\rangle }^2}.$$(8)

The statistical moments that enter these definitions are related to the statistical moments for the random variable δφ or, in the view of equation (1), with the statistical moments of the stochastic component of the medium refractive index δn. In the case of weak scintillation, the expressions (3) and (6) can be used for the determination of scintillation indices.

2.2 Phase screen generation

The phase and amplitude expressions given by equations (1), and (3) are random functions of a random process δn. For the characterization of the wave propagation through the irregular medium, one is usually interested in statistically averaged quantities, such as the scintillation indices, rather than in particular realizations of the amplitude or phase corresponding to a random realization of the random process δn. In this regard, it is convenient to perform the multiple numerical simulations by using the phase screen method and use the resulting set of amplitude and phase values to construct the required statistics.

In the process of generating random phase screens for such simulations, one should match the statistical properties of the screen to those of the random medium. The latter are given in terms of the autocorrelation function of the medium’s refractive index fluctuations, defined in the spatial domain r, or, equivalently, in terms of the corresponding power spectral density (PSD), defined in the Fourier domain of spatial frequencies κ. Several experiments suggest that the spectral characteristics of ionospheric irregularities obey a power-law dependence in the spectral range κ₀ ≤ κ ≤ κ_m. Here, the lower frequency limit κ₀ = 2π/L₀ is related to the largest irregularity scale L₀ called the outer scale, while the frequency κ_m is inversely proportional to the smallest inner irregularity scale l₀. The power law dependence of the spectrum means that the random media can be considered as a composition of irregular structures with sizes varying from the inner up to the outer scale, such that the strength, with which each particular scale contributes to the refractive index fluctuations, depends on the scale size and is governed by the power law dependence. The most versatile analytical model of the three-dimensional PSD of the refractive index fluctuations has been proposed by Shkarofsky (1968) and reads in the case of an isotropic medium as:³

$${\mathrm{\Phi }}_{\delta n}\left(\kappa \right)=\frac{\langle \delta n^2\rangle }{(2\pi {)}^{3/2}}{\left(\frac{{\kappa }_0}{{\kappa }_m}\right)}^{\frac{p-1}{2}}{\kappa }_m^{-3}{t}^{-\frac{p+2}{2}}\frac{K_{\frac{p+2}{2}}\left(t\right)}{K_{\frac{p-1}{2}}\left(\frac{{\kappa }_0}{{\kappa }_m}\right)}, t=\frac{1}{{\kappa }_m}\sqrt{{\kappa }^2+{\kappa }_0^2},$$(9)

i.e., the spectrum depends on the magnitude of the spatial frequency vector $\kappa =\sqrt{{\kappa }_x^2+{\kappa }_y^2+{\kappa }_z^2}$ but not on its direction. Here 〈δn²〉 is the strength parameter of the refractive index fluctuations, p is the (one-dimensional) spectral index, and K_ν(x) is the modified Bessel function of the second kind.

The two-dimensional power spectrum for a thin slab of irregular medium of thickness s in a plane orthogonal to the direction of propagation, i.e., the z-axis, is given by

$$\begin{array}{c}{F}_{\delta \phi }\left({\kappa }_{\perp }\right)=2\pi k^{2}s{\mathrm{\Phi }}_{\delta n}({\kappa }_x,{\kappa }_y,{\kappa }_z=0)\\ =\frac{\langle \delta {\phi }^2\rangle }{\sqrt{2\pi }}{\left(\frac{{\kappa }_0}{{\kappa }_m}\right)}^{\frac{p-1}{2}}{s}^{-1}{\kappa }_m^{-3}{t}_{\perp }^{-\frac{p+2}{2}}\frac{K_{\frac{p+2}{2}}\left(t_{\perp }\right)}{K_{\frac{p-1}{2}}\left(\frac{{\kappa }_0}{{\kappa }_m}\right)}, t_{\perp }=\frac{1}{{\kappa }_m}\sqrt{{\kappa }_{\perp }^2+{\kappa }_0^2}.\end{array}$$(10)

This relation can be obtained using the formula (1). For ease of reference, we write in the second line of equation (10) the explicit formula for the phase spectrum for the Shkarofsky model (9). Here $\sqrt{\langle \delta {\phi }^2\rangle }=\mathrm{ks}\sqrt{\langle \delta n^2\rangle }$ is the variance of the phase screen realizations and ${\kappa }_{\perp }=\sqrt{{\kappa }_x^2+{\kappa }_y^2}$ is the magnitude of the spatial frequency component transverse to the wave propagation direction.

The spectrum (10) can now be used to generate phase screens that have the required statistical properties of the random medium. For this purpose, one uses the well-known properties of a random process, whose realization can be generated by filtering the Gaussian white noise process with a linear filter, followed by an inverse Fourier transform (Rice, 1944). The filter function should be chosen as the square root of the corresponding PSD. Using this method, a phase screen can be generated as:

$$\delta \phi \left({\mathit{r}}_{\perp }\right)=\mathrm{Re}\left[\int \mathrm{}h\left({\mathit{\kappa }}_{\perp }\right)\sqrt{F_{\delta \phi }\left({\kappa }_{\perp }\right)}{e}^{i{\mathit{\kappa }}_{\perp }·{\mathit{r}}_{\perp }}{\mathrm{d}}^2{\mathit{\kappa }}_{\perp }\right],$$(11)

where h(κ_⊥) is a zero-mean Hermitian complex Gaussian random variable with unit-variance, also called as the circular process (Goodman, 2015; Ishimaru, 2017). The expression (11) serves as the basis for several phase screen simulation techniques such as the fast Fourier transformation (FFT) method (Welsh, 1997; Schmidt, 2010), the sparse spectrum approach (Charnotskii, 2013), or the Paulson-Wu-Davies technique (Paulson et al., 2019).

Once the phase screen is generated, the corresponding random realization of the field amplitude is obtained using equation (5) followed by the inverse Fourier transform. When the propagation is simulated by placing multiple phase screens along the path, i.e., within the multiple phase screen approach, the split-step FFT method is usually used (Hardin et al., 1973; Knepp, 1983; Liu et al., 2016; Vasylyev et al., 2022). The propagation path between the transmitter and the receiver is divided into parts whose length is determined by the intersections of the phase screens with the path. The split-step method assumes that the propagation along each segment of the path is equivalent to the propagation in the vacuum. However, each time the signal ray crosses a random phase screen along the propagation path, the initial condition used to simulate propagation in the subsequent (vacuum) path segment, is redefined according to equation (4). The results of such simulations can be further used for obtaining such statistical quantities as the scintillation indices of equations (7) and (8).

3 Phase gradient screens

The phase screen model outlined in the previous section is capable of reproducing refractive scintillation when the length scales of the phase screen variations are much larger than the Fresnel scale size. The functional dependence in equation (1) inherently assumes that the refraction of the signal ray within the phase screen layer can be discarded. If the variations in the refractive index of the medium are characterized by strong gradients, one should also consider the refractive scattering within the medium. For example, the inhomogeneities formed close to the surface of an equatorial plasma bubble have a strong refractive index gradient as being developed on the interface between the depleted electron density region and the ambient ionosphere. Consequently, the electromagnetic wave propagating through the region of the ionosphere occupied by the plasma bubble would experience refractive deflection by both entering and exiting the bubble. To account for these refractive effects caused by strong gradients, the conventional random phase screen should be modified to include the dependence not only on the electron density fluctuations but also on the gradients of such fluctuations.

In this section, we show that this modification consists of adding an additional layer to the conventional phase screen. We will call this layer the phase gradient screen. The generation of the phase gradient screen is based on the gradient field of the initial conventional phase screen. This gradient field is scalar multiplied by the vector field of the refractive deflections of the signal rays from their original directions. We show that the vector field of ray deflections is related to electron density gradients in the topside ionospheric layers.

3.1 Theoretical formulation

Equation (1) implicitly assumes that the fluctuating part of the refractive index satisfies the inequality δn ≪ 1 and that the terms of order δn² and higher can be dropped. It turns out that to account for the refractive effects caused by scattering on refractive index gradients, one must also retain the terms of order δn². To do this, we use the methods of ray (geometric) optics (Tatarskii, 1971; Wheelon, 2004). In the Appendix A some details on the derivation of this correction term are given by solving the eikonal equation up to the δn² order of smallness. The corresponding contribution to the random phase increment is obtained from equation (A23) in Appendix A and can be written in terms of the fluctuating part of the TEC δN_TEC as

$$\delta \phi =\left\{-\lambda r_e\delta N_{\mathrm{TEC}}+k\underset{0}{\overset{s}{\int }}\mathrm{}{\nabla }_{\perp }\delta n\left(s^{\text{'}}{\mathcal{l}}_0\right)·\delta {\mathit{r}}_1\left(s^{\text{'}}\right) \mathrm{d}{s}^{\text{'}}\right\}+\lambda r_e{\nabla }_{\perp }\delta N_{\mathrm{TEC}}·\delta {\mathit{r}}_1\left(s\right),$$(12)

where ℓ₀ is the unit vector pointing in the initial direction of propagation, i.e., ℓ₀ ≡ k/k at the point z = 0. The value of δN_TEC is determined along the straight line given by sℓ₀. The transverse component of the nabla operator is obtained as ∇_⊥ = ∇ − ℓ₀ (ℓ₀ · ∇), while the gradients in equation (12) are evaluated according to the formula (A16). Here and in the following we omit the argument r_⊥ in expressions like equation (12) for the sake of simplicity, but keep this dependency in mind.

The vector

$$\delta {\mathit{r}}_1\left(s\right)=\underset{0}{\overset{s}{\int }}\mathrm{}(s-s^{\text{'}}){\nabla }_{\perp }\delta n\left(s^{\text{'}}{\mathcal{l}}_0\right) \mathrm{d}{s}^{\mathrm{\prime }}=-\frac{r_e{\lambda }^2}{2\pi }\underset{0}{\overset{s}{\int }}\mathrm{}(s-s^{\text{'}}){\nabla }_{\perp }\delta N_e\left(s^{\text{'}}{\mathcal{l}}_0\right) \mathrm{d}{\mathrm{s}}^{\mathrm{\prime }}$$(13)

is the component of the ray curve that is transverse to the original direction ℓ₀, see Figure 1b. Here δN_e (sℓ₀) denotes the fluctuation of the electron density from the mean background value at the point (r_⊥, s), where r_⊥ is transverse to the initial propagation direction ℓ₀. The magnitude of the vector δr₁ indicates the deflection distance of the ray traveling the path length s in the medium from the original direction ℓ₀.

Figure 1 Scheme of typical effects experienced by the electromagnetic wave propagating through a slab of randomly inhomogeneous medium of thickness s. The conventional phase screen approach simulates the random corrugation of the incoming wavefront (a). The phase derivative screen adds the random refraction of the incoming signal ray (b). In the case of a sufficiently strong refractive index gradient ∇δn, the ray bends within the layer and appears to be displaced at distance δr1 from the original direction of propagation ℓ0. The dashed vector shows that, in contrast to case (a), a different ray emerges in the initial direction of ray propagation. The fluctuating part of the gained phase along the initial direction is then given by equation (14).

In Appendix B we show that the mean value of the phase, cf. equation (12), vanishes for quite ionospheric conditions when the power-law spectral model, cf. equation (9), is applicable. This implies that the variables ∇_⊥δn and δr₁ can be considered as independent random variables with zero mean. This condition may be violated under disturbed space weather conditions, but we exclude such situations from our further consideration.

Equation (12) shows that if the second-order terms in the smallness of δn are retained, the phase gained by propagating through the layer of the randomly inhomogeneous ionosphere is not only proportional to the TEC calculated along the straight ray path but also depends on the transverse component of the TEC gradient. Thus, if the refractive scattering within the layer of random medium is significant to be accounted for, the resulting scintillation of the electromagnetic wave will correlate with the TEC gradients. In particular, one can expect the correlation of the phase and amplitude scintillation indices not only with the standard deviation of the TEC, but also with the standard deviation of the TEC gradients. Additionally, one can also notice that the deflection vector, given by equation (13), is expressed in terms of the gradient of the electron density fluctuations. This dependence explains the correlation between the scintillation indices and the ionospheric indices based on electron density gradients.

3.2 Generation of phase gradient screens

To include the refractive scattering on electron density gradients in the phase screen generation algorithm, we note that equation (12) is derived from the expansion of the type, cf. with equation (A20),

$$\delta \phi \left(s{\mathcal{l}}_0\right)=\delta \phi \left(\right(s{\mathcal{l}}_0+\delta {\mathit{r}}_1)-\delta {\mathit{r}}_1)\approx \delta \phi (s{\mathcal{l}}_0+\delta {\mathit{r}}_1)-\delta {\mathit{r}}_1·{\nabla }_{\perp }\delta \phi (s{\mathcal{l}}_0+\delta {\mathit{r}}_1).$$(14)

Here, the first term on the right-hand side corresponds to the terms in the curly bracket of equation (12), while the second term corresponds to the term proportional to the gradient of the TEC variation in equation (12). The expression (14) gives the random phase increment gained by the electromagnetic wave propagating in the random medium of thickness s when the random refraction within the medium layer cannot be neglected, i.e., when δr₁ ≠ 0. Due to the accumulated internal refraction the initial signal ray changes its original propagation direction from ℓ₀ to ℓ₀ + δr₁/s. The phase increment at the point along the initial direction is expressed as the phase increment along the changed direction corrected by the term proportional to its gradient. Since the displacement vector δr₁ in the considered approximation is transverse to the original direction ℓ₀, the gradient operator in equation (14) is also transverse to this direction and would be further denoted as ∇_⊥δφ.

The expression (14) suggests that the original ordinary phase screen generation method of Section 2.2 can be extended as follows. It is crucial to assume that the random variables ∇_⊥δn and δr₁ are considered to be statistically independent, such that the terms δr₁ and ∇_⊥δφ in (14) can be generated independently from each other. The validity of such an assumption is discussed in detail in Appendix B.

Now consider the phase screen plane transverse to the initial propagation direction of the electromagnetic wave ℓ₀. The coordinate system is chosen in a similar way as in Section 2.1 so that the z-axis is chosen along ℓ₀ and any point within the plane is given by the vector r_⊥ = (x y)^T. The ordinary phase screen is then generated by using the procedure of Section 2.2. Then the corresponding phase gradient screen is obtained by taking the gradient as

$${\nabla }_{\perp }\delta \phi \left({\mathit{r}}_{\perp }\right)=\mathrm{Re}\left[i\int \mathrm{}{\mathit{\kappa }}_{\perp } h\left({\mathit{\kappa }}_{\perp }\right)\sqrt{F_{\delta \phi }\left({\kappa }_{\perp }\right)}{e}^{i{\mathit{\kappa }}_{\perp }·{\mathit{r}}_{\perp }}{\mathrm{d}}^2{\mathit{\kappa }}_{\perp }\right].$$(15)

This procedure can be thought of as taking the gradient from the image representing the phase screen realization (Sect. 2.2).

The random displacement vector δr₁ can be generated independently of ∇_⊥δφ by assuming their statistical independence. If the transverse component of the refractive index gradient is |∇_⊥δn| ≤ M along the entire propagation path sℓ₀ in medium, then clearly $\left|\delta {\mathit{r}}_1\right|\le M{\int }_0^s\mathrm{}(s-s^{\text{'}})\mathrm{d}{s}^{\text{'}}=s^{2}M/2$. Thus, the magnitude of the displacement vector is bounded by the greatest absolute value of the refractive index gradient along the propagation path. This suggests to approximate the random displacement (13) as

$$\delta {\mathit{r}}_1\approx \frac{s^2}{2}{\nabla }_{\perp }\delta n\left(s_{\mathrm{ref}}\right),$$(16)

where s_ref is some reference distance along the path sℓ₀, where ∇_⊥δn reaches its maximum value. It is reasonable further to make the following approximation:

$$\delta {\mathit{r}}_1\approx \frac{s^2}{2}{\nabla }_{\perp }\delta n\left(0\right)=-\frac{r_e{s}^2{\lambda }^2}{4\pi }{\nabla }_{\perp }\delta N_e({\mathit{r}}_{\perp },z=0).$$(17)

In fact, the refractive index gradients at the boundary of the random medium, s_ref = 0, introduce the ray deflections that grow as the rays propagate through the medium. The gradients within the medium may be higher than at the boundaries, but their contribution to the cumulative ray displacement may still not be as large as that produced by the gradient at the surface of the layer. Strictly speaking, this condition can of course be violated for very steep gradients occurring within the bulk random medium, so equation (16) is more restrictive than the assumed condition (17). Some additional arguments of the importance of the topside ionospheric gradients in producing the major contribution to the signal ray deflections and phase scintillations associated with the tilt errors are given in Appendix C.

The approximations (16) or (17) seem to contradict our assumption that ∇_⊥δn and δr₁ are statistically independent. Within the random phase generation algorithm of this section, this is not an issue, since in this case it is assumed that δr₁ and ∇_⊥δφ are statistically independent. Since the phase is a cumulative random variable, i.e., it is an integral over δn, its correlation with any particular value δn(s_ref) is weak, especially if the propagation path s is long, as it is in the case of trans-ionospheric propagation.

In the present paper, we focus our attention on the method of phase gradient screen generation, where the simulation is supported by in-situ measurements of the electron density gradients by Swarm satellites. The proposed algorithm for the generation of random phase gradient screens is summarized in Table 1. The generated phase gradient screen can then be used in electromagnetic wave propagation simulations, for example by using the split-step method, as it is done for conventional phase screen simulations. The method based on a purely theoretical model of the electron density spectral characteristics would be considered elsewhere.

4 Simulation examples

In this section, we give some illustrative examples of scintillation simulation using the phase gradient screen method. To emphasize the difference between the conventional random phase screen approach and the proposed one, we apply both methods to the toy example of signal propagation through an irregular layer of ionosphere with the embedded region of depleted electron density. The remaining examples validate the phase gradient method by using scintillation data and other auxiliary information for high and low-latitude regions.

4.1 Toy example

Let us illustrate the capabilities of the phase gradient screen algorithm to simulate refractive scintillation by considering the following toy example, see Figure 2a. We assume that the region of depleted electron density is immersed in the ionospheric layer. For simplicity, we neglect the curvature of the Earth and the ionospheric layer and assume that the depletion region has a cylindrical shape with the axis aligned along the vertical direction. The height of the cylinder is assumed to be less than the thickness of the ionospheric layer. We suppose that both the ionospheric layer and the depleted region are randomly inhomogeneous media.

Figure 2 Schematic geometric arrangement for the toy example of cylindrical depletion region embedded in the slab of the ambient ionospheric medium (a). Three typical scenarios of signal propagation through the considered ionospheric structure are marked with i, ii, iii. Namely, in case i, the refractive scattering of the signal occurs at the top and bottom sides of the depleted structure; in case ii, the additional scattering occurs at the vertical wall of this structure, while in case iii, no such scattering occurs since the depleted region is absent. An example of the vector fields used to construct the phase gradient screen is shown in (b) and (c). The screen itself is placed perpendicular to the initial propagation direction of the signal wave and has a circular cross-section with the cylindrical depletion region. Figure (b) shows the contribution to the displacement vector field δr1 that comes from the scattering on the lateral surface of the cylinder, cf. the ray path ii in scheme (a). Figure (c) shows the generated phase screen with random amplitude δφ and the corresponding gradient field ∇⊥δφ shown as black arrows.

The scintillation of radio signals transmitted through the model ionospheric medium is simulated using the ordinary and gradient phase screen approaches described in Sections 2.2 and 3.2, respectively. In both cases the phase screens are placed transverse to the vertical direction and their cross-section with the depletion region would be a circular disk. The electron density perturbations are smaller in high ambient density regions than in the depleted regions (Huang et al., 2014). Therefore, we assume that the strength parameter of the refractive index fluctuations, 〈δn²〉, is larger inside the circular region within the screen compared to the outer ambient electron density region. For definiteness, we set $\sqrt{\langle \delta n^2\rangle }$ for the depleted region to be 5/4 of the value for the ambient ionosphere.

The generation of the ordinary phase screen follows the procedure described in Sections 2.2. An example of the generated phase screen and the resulting phase scintillation pattern at the observation plane located 350 km behind the screen are shown in Figures 3a and 3c. When the ordinary phase screen is generated we can extend it to the phase gradient screen according to equation (14). Firstly, we construct the gradient field ∇_⊥δφ using the values of δφ generated by the conventional procedure, see Figure 2c. Secondly, we calculate the vector field for the deflection parameter δr₁ according to equation (17). Here the gradient of the refractive index fluctuations, ∇_⊥δn, is evaluated by using the distribution of the fluctuation strength parameter 〈δn²〉 according to the model described above. In the considered model with sharp boundaries between the regions of the depleted and ambient electron densities, the deflection vector is non-zero at the boundaries of the circular region, see Figure 2b. The realization of the phase gradient screen according to equation (14) is shown in Figure 3b.

Figure 3 Comparison of the simulation results of the split-step algorithm using the conventional random phase screen approach (a), (c), (e) and the phase derivative method (b), (d), (f). Figures (a) and (b) show the random phase screens generated by the respective methods. Multiple simulations of radio wave propagation through these screens yield the corresponding phase scintillation indices (c) and (d). Here, the red dotted line shows the possible satellite ground track and the figures (e), (f) show the corresponding scintillation index records along the track. The markers i, ii, iii along the track in plot (d) refer to the propagation scenarios shown schematically in Figure 2 (a).

Using the multiple simulations of radio signal propagation through such a screen and collecting statistics on the wave phase fluctuations recorded at the observation plane we calculate the phase scintillation index according to the definition (8). The simulated phase scintillation indices are shown in Figures 3c and 3d for the ordinary and the phase gradient phase screen approaches. In both cases the phase scintillation is higher in the circular region as a consequence of the higher level of refractive index fluctuations in the region of depleted electron density.

The simulation with phase gradient screens shows the additional increase of the scintillation index at the circular edge of this structure, the effect attributed to the refractive scattering on steep electron density gradients. The jump-like increase of the phase scintillation index in situations, when the signal wave crosses the edge of a large-scale ionospheric structure, has been reported in References (Ledvina et al., 2002; Li et al., 2010; Prikryl et al., 2011; van der Meeren et al., 2014; Nishimura et al., 2021; Spogli et al., 2023). According to the phase gradient screen approach, such an increase in scintillation is related to the refractive scattering at the interfaces between large-scale structures embedded in the ionosphere. The physical interpretation of refractive scattering at such interfaces is as follows: the region at or near the interface is the region where ionospheric irregularities are likely to develop, e.g., due to Rayleigh-Taylor or Kelvin-Helmholtz instabilities. The irregularities with scales larger than the Fresnel scale contribute to the refractive scattering of the propagating signal wave and increase the associated refractive scintillation levels recorded by the receiving station.

The example considered in this section is rather artificial and is used to illustrate the main feature of the phase gradient screen approach. The method yields enhanced scintillation levels due to scattering on ionospheric structures characterized by strong electron density gradients. In the given example such a gradient appears on the border of the cylindrical depletion region. This border is regular in the sense it is given in the analytical form. In realistic situations, the boundary of the large-scale ionospheric formation, e.g., SED, EPB, etc., would have an irregular structure. At the interface of such a structure with the ambient ionosphere, various ionospheric irregularities of different scales form and evolve, some of which, characterized by large gradients, contribute to refractive scintillation. In the following sections, we consider some examples of such structures.

4.2 Scintillation over the European sector under disturbed conditions

In this section we consider the European region during a severe geomagnetic storm (a level 4 of 5 on NOAA’s space weather G-scale) on April 23, 2023. The storm was caused by a coronal mass ejection that erupted from the Sun on April 21 and reached its peak on April 23 at about 19:00 UTC. The Kp-index reached a value of 8 at 18:00 UTC on April 23, and this value remained above 5 until 12:00 UTC on April 24.

Some minor phase scintillation activity was observed in the night sector over northern Europe during the period April 23–24. The colored squares in Figures 4 and 5 show empirical values of the phase scintillation index obtained by some GNSS receiving stations located in Europe. The geographic coordinates of each square correspond to the coordinates of the ionospheric pierce points of the corresponding GNSS satellite link to the ground station.

Figure 4 Comparison of simulated phase scintillation indices (colored circles along the orbit of the Swarm satellite) with the empirical values (colored squares) for L1 radio signals. Simulations are performed using the conventional random phase screen technique (a) and by using the phase gradient screens (b). The empirical values are obtained from several GNSS receivers located in Europe and their geographical locations are given by the position of their ionospheric piercing points. The black arrows correspond to the NeGIX vectors (arbitrary scale), cf. Appendix D for the NeGIX definition. For reference, the contours for the magnitude of the vertical TEC gradient (in units of mm/km on GPS L1) are also shown.

Figure 5 Phase scintillation indices are simulated (circles) within the phase gradient screen approach and compared to empirical values (squares) like Figure 4b, but for the post-storm period for two successive passes of the Swarm satellites.

We use the phase gradient screens to simulate the phase scintillation over this region. The procedure of screen generation follows the procedure outlined in Section 3.2 and equation (14) defines the values of the random phase over a generated phase screen. While the phase δφ and its gradient ∇_⊥δφ are simulated using the power-law spectral model according to equation (10), the deflection parameter δr₁ in equation (14) is determined from the empirical value for the electron density gradient according to the approximation (17). By using this approximation it is assumed that the steep electron density gradient at the topside of the ionospheric layer, which is relevant for scintillation formation, contributes to the deflection of the propagating signal ray on distance |δr₁| at the bottom of this layer.

If we place the phase screen at the bottom of the scintillation-producing layer, say at the bottom of the ionospheric F-layer, and relate the thickness of the layer to the outer scale of ionospheric electron density fluctuations, we can estimate the height at which electron density gradients are of interest for the modeling. Various studies give different values for the length of the outer scale, ranging from 50 km (Vats et al., 1981) up to 200 km (Rino, 1979), with the preference to use values within these limits (Patel et al., 2011; Aol et al., 2020). Taking these considerations into account, it seems that the in-situ measurements from the Swarm satellites are optimal for the needs of scintillation modeling with the phase gradient screens.

The Swarm constellation consists of the Alpha, Bravo, and Charlie (A, B, and C) satellites, which were launched into near-polar orbits on November 22, 2013. A detailed description of the Swarm satellite mission can be found in (Friis-Christensen et al., 2006). Swarm A and C fly side by side (1.4° separation in longitude at the equator) at an altitude of 462 km (initial altitude) and at 87.35° inclination angle. Swarm B is decoupled from satellites A and C and flies in a higher orbit of 511 km (initial altitude), and due to these characteristics the data from this satellite are not suitable for estimating gradient values relevant for scintillation simulation purposes.

Appendix D summarizes the definition of the electron density gradient index (NeGIX) obtained from in-situ measurements by Swarm A and C satellites. The mean value of the deflection vector for the propagating signal ray at the bottom of the scintillation-producing ionospheric layer is related to the NeGIX index as

$$\langle \delta {\mathit{r}}_1\rangle \approx -\frac{r_e{s}^2{\lambda }^2}{4\pi }\langle \nabla \mathrm{Ne}\rangle .$$(18)

This relation follows from equation (17) with the NeGIX value 〈∇N_e〉 defined in equation (D3). For the NeGIX value 〈∇N_e〉 = 250 el cm⁻³ km⁻¹ and assuming that the thickness of the layer is 200 km, e.g., for the low-latitude ionosphere, the deflection distance for the L1 signal is 8 cm. Even for extremely high gradient values, the deflection distance reaches the values of half of a meter. The smallness of this parameter allows us to neglect it in the argument of the random phases δφ on the right-hand side of equation (14). This means that we can proceed with the generation of phase gradient screens by assimilating the in-situ gradient values according to Table 1, where we now substitute the mean value for the deflection vector according to equation (18).

The range of NeGIX values is related to the relative position of Swarm satellites A and C and to the data acquisition rate. These two factors limit the extent of ionospheric scales to the range of 30–150 km. The outer scale of ionospheric scintillation (Patel et al., 2011) is comparable to the scales within this range, and thus the NeGIX is the proper in-situ data for phase scintillation modeling. However, the empirical data for scintillation index estimation are processed in such a way that any trend components and multipath contributions in the signal records are minimized or eliminated. This is accomplished by setting the low-frequency cutoff equal to 0.1 Hz, which is achieved by the appropriate filtering procedure. The effect of introducing the low-frequency cutoff is equivalent to ignoring scintillation-associated scales larger than about one kilometer. The values of the spatial gradients associated with such kilometer-scale structures are higher than the NeGIX index derived for scales of two orders of magnitude larger. It is reasonable to assume that the values of gradients for kilometer-scale inhomogeneities are two orders of magnitude larger than the NeGIX values. This assumption is well supported by the electron density gradients of kilometer-scale structures reported in Loucks et al. (2017). In the following, we take this scaling factor into account when simulating the scintillation levels. This allows us to compare the results of the simulations with the empirical scintillation indices.

Figures 4 and 5 demonstrate the scintillation modeling capabilities of the phase gradient method. These figures show the empirical values of the phase scintillation indices for the L1-band signals from GNSS satellites (colored squares) that have been recorded by the Javad receivers installed in Kiruna, Sweden (67.84° N, 20.41° E); Ramfjordbotn, Norway (69.58° N, 19.22° E); Teneriffe, Spain (28.48° N, −16.32° E); Toulouse, France (43.56° N, 1.47° E); Tromsø, Norway (69.68° N, 18.98° E). Since the propagation geometry for phase screen simulations assumes vertical propagation, only data from satellites with elevation angles greater than 30° were used for validation purposes. For mapping of scintillation values the ionospheric pierce point (ipp) heights of 250 km have been used for high latitude stations (latitude >60° N) and the ipp-heights of 400 km for stations in mid and low latitudes. The ground tracks of the midpoints between Swarm A and C satellites are filled with colored circles showing the phase scintillation indices simulated along these tracks.

For reference, we also show the NeGIX vectors (black arrows), arbitrarily scaled for ease of presentation. The contour lines also show the values of the TEC gradients in units of mm/km on GPS L1. The TEC gradients are estimated from the ground-based TEC measurements assimilated into the NTCM model (Jakowski et al., 2011). Due to the coarse spatial resolution of the NTCM model, the use of TEC gradients in the scintillation modeling with phase gradient screens is not possible, but these contours show the variability of ionospheric conditions under geomagnetic disturbances. Note that the TEC gradient contours are missing from the upper left corners of the plots. This is an artifact related to the unreliable TEC values caused by the calibration problems for this region of the disturbed ionosphere. The boundary of this region correlates with the equatorward boundary of the auroral oval. The increase of the NeGIX magnitude and the remarkable variability of the gradient directions in the high latitude region are also related to the plasma structuring at the boundary of the auroral oval and in the regions of auroral emissions.

In Figure 4 we compare the simulation capabilities of the ordinary (Fig. 4a) and the phase gradient screen (Fig. 4b) methods under quite geomagnetic activity (Kp-index < 2). The empirical phase scintillation index shows elevated values in the northern and southwestern parts of Europe. The NeGIX gradient values show an increase in the region of southern and northern Europe, with significantly higher values and large variability in gradient directions over the northern sector of Europe. The ordinary phase screen simulation technique shows a uniform distribution of scintillation index values along the Swarm track, cf. Figure 4a. In addition, the phase gradient screen approach shows some variability in the northern part of the track due to the gradient component proportional to NeGIX. When compared with the empirical scintillation index, it can be seen that the phase gradient screen approach shows a better resemblance to the distribution of the empirical scintillation data. The simulation parameters used to model the phase scintillation indices along the track are summarized in Table 2.

Figure 5 shows the Swarm passages during the storm’s recovery phase along with the recorded phase scintillation indices. Similar to Figure 4, the elevated values of scintillation indices and the electron density gradients are associated with the edge of the auroral oval. Under disturbed conditions the oval boundaries are shifted more towards the equator compared to Figure 4. One can also see that the scintillation levels are smaller in the post-storm recovery phase compared to the quiet conditions shown in Figure 4. The increased scintillation activity can also be seen in the middle and lower parts of the map. These events seem to be intermittent in nature and are associated with small-scale patches (scale size < 30 km) that are not visible on the NeGIX traces. In this respect, the electron density gradients along the Swarm satellite tracks calculated from high-rate (16 Hz) electron density recordings would provide valuable information. However, for the times in this work (shown in Figs. 4–6) both Swarm A and C satellites were in a configuration did not allow reliable estimating of the plasma density from faceplate current measurements, and the 16 Hz density product could not be used in this study.

Figure 6 Comparison of the simulated amplitude scintillation indices (circles) with the empirical data (squares) for the low latitude region over Africa, the Atlantic Ocean, and South America. As a reference to the possible plasma bubble shapes, the background image shows the atmospheric emission radiance obtained from the GOLD mission data. The green vectors represent the gradient values of the mask extracted from the GOLD images and the black vectors correspond to the Swarm NeGIX index.

4.3 Low latitudes

The phase gradient screen approach is applicable to situations where steep electron density gradients cause random refraction of the propagating radio signal. Amplitude scintillation is primarily related to the diffractive effects. In this respect, the irregularities with sizes equal to the Fresnel scale $\sqrt{\lambda L}$ contribute the most to the perturbation of the signal amplitude. For example, in the case considered above of scintillation over the European sector under undisturbed conditions, cf. Figure 5, the amplitude scintillation does not show much variability in the high latitude sector, with the S₄ value remaining below the weak scintillation level of 0.2. The values of the index S₄ do not correlate with the gradient index NeGIX, which shows the effect of “phase without amplitude scintillation” (Mushini et al., 2014; van der Meeren et al., 2014; Priyadarshi et al., 2018). This is because the ionosphere in this sector is structured on scales larger than the Fresnel scale, and the refractive scattering has no effect on amplitude scintillation.

The situation can change in the case when the variance of the phase fluctuation reaches large values (Crain et al., 1979; Booker & MajidiAhi, 1982; Forte, 2008) such that the value of its root mean square error (RMSE) satisfies $\sqrt{\langle \delta {\phi }^2\rangle }>>1$ rad. This situation corresponds to a strong scattering regime, in which the amplitude scintillation index attains high values, reaching the saturation level of S₄ = 1. The power spectral density for amplitude scintillation shows a roll-off at higher spatial frequencies than the expected Fresnel spatial frequency $2\pi /\sqrt{\lambda L}$. At the same time, the wide range of spatial scales contributes equally significantly to the amplitude fluctuations. In particular, this range now includes scales associated with purely refractive scattering and scintillation.

The conditions described above may be well met in the low-latitude ionosphere after local sunset at ionospheric heights and under high solar activity conditions (Balan et al., 2018; Li et al., 2020). In particular, radio signals propagating through the edges of equatorial plasma bubbles (EPBs) exhibit intense amplitude scintillation (Sousasantos et al., 2022; Seechai et al., 2023; Spogli et al., 2023). We apply the phase gradient technique to simulate of such intense amplitude scintillation at the edges of EPBs.

We use the night disk scan (NI1) observations obtained by the Global-scale Observations of the Limb and Disk (GOLD) mission (Eastes et al., 2017). The data of 135.6 radiance values in Rayleighs of post-sunset atmospheric emission (atomic oxygen) have been combined into a single image from two channels of the GOLD imaging unit. Over the low-latitude region, the resulting picture can provide a good image of the EPB contours. With the approach adopted here, we do not reconstruct the entire shape of each plasma bubble, but rather its boundaries within the southern and northern equatorial ionization crests. This simplification is partly justified because the difference in electron density values at the edges of the EPBs is greatest in these regions, as are the corresponding values for the electron density gradients. To determine the boundary between the plasma bubbles and the surrounding ionosphere within the crest regions, we use the GOLD radiance images to generate the corresponding binary masks. The idea of using such masks is similar to that used in Section 4.1, cf. with the binary mask in Figure 2b. We also calculate the gradient fields for the obtained masks in a similar way as was done in Section 4.1. We set the radiance value of 120 Rayleighs as the threshold value for the mask construction.

Another gradient field needed for phase gradient screen simulations is again provided by the Swarm’s NeGIX index. The index is well suited for detecting plasma bubble boundaries, especially when one of the tandem satellites enters the region of depleted electron density. Figure 7a shows the situation when the Swarm satellites A and C cross different regions of the same plasma bubble or different plasma bubbles. It can also be seen that there are situations where one satellite detects an EPB while the other does not, or vice versa. As the satellites enter such regions, the direction of the NeGIX vector changes to the opposite one, cf. Figure 7b. The increase in magnitude of this vector is typical in the presence of EPBs. We also note that the FORMOSAT-7/COSMIC plasma density measurements are also an attractive source of information on EPB gradients (Zakharenkova et al., 2023) that can be used for scintillation modeling. In addition, these measurements provide information on plasma bubble structuring in a direction complementary to that along the Swarm orbit.

Figure 7 Profiles of the electron density along the orbits of Swarm A and C satellites (a) and the corresponding variation of the NeGIX (b), cf. also with Figure 6e. The difference in the records and the corresponding variability of the gradient index show that the satellites crossed different regions of the EPBs. The maximum magnitude of the NeGIX vector shown in (b) is 11571 e cm−3 km−1.

We simulate the scintillation as follows. We use the Global Ionospheric Scintillation Model (GISM) (Béniguel & Hamel, 2011) to simulate the background scintillation levels in the equatorial region. This model is based on the multiple-phase screen approach and is able to account for the strong scattering regime typical for the region under consideration. However, due to its climatological nature, the GISM is not able to reproduce patchy scintillation patterns caused by EPBs. To model such an effect, we use an additional phase gradient screen and simulate scintillation due to steep gradients. The additional screen is used only in the regions of the plasma bubble boundaries that are determined from the GOLD irradiance images as described above. The phase gradient screen consists of the conventional phase screen generated by using the parameters given in Table 2 and the gradient part consisting of the product of NeGIX with the gradient field of the binary mask extracted from the GOLD data. The numerical coefficient in front of the product term is hardly related to any physical parameter, and it is chosen by fitting the simulation results to empirical scintillation values, which yielded the value of 2.8 × 10⁻⁵ in arbitrary units.

The results of the simulations with the GISM and with the phase gradient screen are then superimposed and shown in Figure 6 for the post-sunset ionosphere over Africa, the Atlantic Ocean, and South America during the period November 3–6, 2022. The empirical amplitude scintillation indices were obtained from GNSS stations in Fortaleza, Brazil (−3.74° N, −38.58° E), São Paulo, Brazil (−23.54° N, −46.63° E), and on Tenerife, Spain (28.48° N, −16.32° E). Figures 6a–6c show several passes of Swarm A and C satellites over this region during the hours after sunset on November 3. These figures illustrate the behavior of the NeGIX gradient index as EPBs begin to develop. In Figure 6a⁴ a small increase of NeGIX can be seen in the crests region, which coincides with the emission region from the GOLD images. A similar situation can be seen in Figure 6b, where the behavior of the gradients is more regular. Figure 6c shows that the regular NeGIX behavior changes to a more chaotic one when approaching the northern ionization crest. In this case, the Swarm tracks cross several plasma bubbles. The largest gradients are seen in the region of the northern crest, where the difference between depleted and ambient electron densities is largest. As expected, the simulation results show a rapid increase in the scintillation index in this region.

Figures 6d–6f show examples of moderate (0.4 < S₄ < 0.6) and strong (S₄ ≥ 0.6) amplitude scintillation caused by plasma bubbles. Again, high values of the empirical scintillation indices correlate well with the regions characterized by strong gradients, especially at locations where the gradient direction reverses abruptly, i.e., near the edges of EPBs. Scintillation modeling with phase gradient screens is able to reproduce this behavior. The proposed simulation approach can be further refined if the shape of the EPBs is known. Figure 8 shows the excerpt of Figure 6e with simulated scintillation indices along the Swarm track using three methods: the GISM model (a), the GISM and phase gradient screen method with the binary mask (b), and with a more complex mask that aims to account for the plasma bubble structures between two crest regions (c). The latter mask, shown as a gray background in Figure 8c, uses the additional threshold value for atmospheric irradiance (10 Rayleighs) in addition to the value of 120 Rayleighs as used in Figure 8b. The coarse spatial resolution of the GOLD L1C images does not allow reliable determination of the boundaries between plasma bubbles in the region between the northern and southern crests. However, even the implementation of such a coarse mask as used in Figure 8c can improve the ability of the phase gradient screen approach to simulate EPB-related scintillation.

Figure 8 Comparison of methods for simulating equatorial scintillation in the presence of EPBs: (a) GISM simulation yielding the background scintillation values, (b) and (c) GISM simulations with the additional phase gradient screens. The difference between (b) and (c) is related to the masks used (shown as gray areas on the background). While (b) uses a simple binary mask that determines the irradiance regions with a single threshold level of 120 Rayleigh, the mask in (c) uses two thresholds at 120 and 10 Rayleigh in order to incorporate some EPB morphology between two crests of ionization.

5 Conclusions

This study focuses on the modeling of refractive scattering to scintillation phenomena in an irregular ionosphere. This contribution is not negligible when the ionospheric inhomogeneities and large-scale structures are characterized by very steep spatial gradients of electron density or the TEC. We show that the next-order correction to the geometric-optical approximation should be included in the theoretical consideration of trans-ionospheric radio wave propagation. This contribution is of the second order of smallness with respect to the fluctuating part of the ionospheric refractive index δn and is approximately equal to the scalar product of the spatial gradient of the electron density at the topside of the scintillation-producing ionospheric layer and the spatial gradient of the TEC. This dependence can explain the reported good correlation between scintillation indices and the gradient-based proxy indices such as ROTI, RODI, AATR, GIX, etc. In fact, some of these indices are expressed in terms of the TEC gradients (either spatial or temporal), e.g., GIX, or in terms of the electron density gradients, e.g., RODI, and the dependence on both types of indices is present in the correction term of the geometric optics approximation. In the present article we have used the electron gradient index NeGIX derived from the Swarm in-situ measurements of the electron density to assist simulations of scintillation events in low- and high-latitude regions.

Using this correction, we extend the conventional phase screen algorithm for simulating radio wave propagation in a randomly inhomogeneous ionosphere according to the so-called phase gradient screen approach. While the standard phase screen, placed along the wave propagation path at a certain distance from the signal receiver, simulates random phase advance or retardation due to scattering on ionospheric irregularities, the next-order correction incorporated in the phase gradient screen also takes into account the random refractive bending of the propagating signal rays. The ray bending results in a small tilt of the local wavefront, and the superposition of several such tilts results in an additional corrugation of the overall signal wavefront. This additional effect mainly contributes to the increase of the phase scintillation index, but in the case of strong scattering, it can also lead to an increase of the amplitude scintillation.

The primary application area of the phase gradient screen approach is in the high latitude ionosphere, where the refractive scintillation is greater than the diffractive scintillation due to scattering on steep spatial gradients. The resulting difference between the phase and amplitude fluctuation levels is usually referred to as “phase without amplitude” scintillation. For this situation, the weak scattering approximation is applicable, so that the RMSE of the phase fluctuations of the random phase screen is small ($\sqrt{\langle \delta {\phi }^2\rangle }<<1$ rad) and a single-phase screen approach can be used to simulate the wave propagation. On one hand, if steep gradient structures have scales larger than the Fresnel scale $\sqrt{\lambda L}$, they will not affect amplitude fluctuations. On the other hand, the range of scales contributing to phase scintillation is broader and includes irregularities with sizes $>\sqrt{\lambda L}$ corresponding to steep gradient structures. As a result, it is sufficient to use the conventional phase screen approach for modeling of amplitude scintillation, while the additional phase gradient screen correction may be required for simulating phase scintillation. For illustrative purposes, we considered scintillation events over the European sector before and after the onset of the geomagnetic storm on April 24, 2023. Simulations and empirical data suggest that the enhancement of phase scintillation occurs at the locations associated with the polar auroral emission and auroral particle precipitation regions.

For the wave propagation in the low-latitude ionosphere, the weak scattering approximation is not always valid. In this case, one could use either a single phase screen with $\sqrt{\langle \delta {\phi }^2\rangle }>>1$ rad or multiple phase screens with small phase variances, as implemented, e.g., in the GISM model. For this strong or multiple scattering regime it turns out that scales larger than the Fresnel scale are found to contribute to amplitude scintillation. In the presence of steep gradient ionospheric structures, both amplitude and phase scintillation may require phase gradient screen correction. We have illustrated this by considering amplitude scintillation enhanced by the presence of equatorial plasma bubbles in the post-sunset ionosphere over the low-latitude geographic sector. The GISM model, enhanced by a properly designed phase gradient screen, is able to reproduce the background scintillation levels modulated by the ionization crest regions as well as the patchy structures of enhanced scintillation caused by scattering on highly irregular plasma on the edges of plasma bubbles. To simulate plasma-bubble-associated scintillation we have used in-situ Swarm data, which have very sparse spatial coverage of the region of interest. However, our results show that the incorporation of the semi-empirical plasma bubble model into the scintillation modeling algorithm with the phase gradient screen component is quite promising and would pave the way towards more accurate and reliable simulation and prediction capabilities of intense equatorial scintillation.

Acknowledgments

This work is partially funded by ESA via the Swarm DISC Subcontract SW‐CO‐DTU‐GS‐133. Dmytro Vasylyev acknowledges the programmatic funding of the German Aerospace Center within the HiGAIN project. The editor thanks Weijia Zhan and an anonymous reviewer for their assistance in evaluating this paper.

Data availability statement

Electron density measurements acquired by Swarm A and C have been taken from the publicly available Swarm Data Access at https://swarm-diss.eo.esa.int. The data correspond to Langmuir Probe (LP) measurements of the Electric Field Instrument (EFI) with 2 Hz resolution and are available at the following link: https://earth.esa.int/eogateway/missions/swarm/data. Scintillation data has been accessed via the IMPC scintillation indices archive data at https://impc.dlr.de/products/ionospheric-perturbations/local-scintillation-measurements. Scintillation data from Saõ Paulo are taken from the electronic Space Weather upper atmosphere (eSWua, http://eswua.ingv.it) portal managed by the INGV Upper atmosphere physics and radiopropagation Working Group. TEC gradient data have been accessed via the ESA Space Weather Portal (https://swe.ssa.esa.int/techtide-federated). GOLD UV data are available at the GOLD Science Data Center (https://gold.cs.ucf.edu/data).

Appendix A: Solution for the eikonal function using the perturbation technique

This appendix summarizes, for ease of reference, the perturbation technique for determining the phase of an electromagnetic wave propagating in a random medium. Since we are interested in the contributions to phase fluctuations due to random refraction on medium inhomogeneities, the approximation of geometric optics could be used. Here we closely follow the considerations elaborated in Tatarskii (1971).

We start with the homogeneous wave equation in a medium with a given refractive index $\tilde{n}$

$$\Delta U\left(\mathit{r}\right)+k_0^2 {\tilde{n}}^2\left(\mathit{r}\right)U\left(\mathit{r}\right)=0,$$(A1)

where Δ is the Laplace operator and k₀ is the wavenumber of the signal wave in vacuum. The refractive index $\tilde{n}\left(\mathit{r}\right)$ is considered a random function of the spatial coordinates r with mean value $\langle \tilde{n}\rangle$. For further convenience, we rewrite (A1) in a slightly different form

$$\Delta U\left(\mathit{r}\right)+k^2 n^2\left(\mathit{r}\right)U\left(\mathit{r}\right)=0,$$(A2)

where $k=k_0\langle \tilde{n}\rangle$ is the wave vector in the medium and $n=\tilde{n}/\langle \tilde{n}\rangle$ is the scaled refractive index. The fluctuating part of the refractive index is denoted by $\delta \tilde{n}=\tilde{n}-\langle \tilde{n}\rangle$ and the corresponding scaled variable is δn = n − 1. The wave function can then be written as

$$U\left(\mathit{r}\right)=A\left(\mathit{r}\right) e^{\mathrm{ik\theta }\left(\mathit{r}\right)}.$$(A3)

Here A and θ are slowly varying functions⁵. In the framework of geometric optics the amplitude A is expanded into a series of the small parameter 1/k, i.e., A = A₀ + k⁻¹A₁ + k⁻²A₁ +…. Substituting this expansion together with equations (A3) into (A1) and collecting equal powers of k, one obtains from the expression for the coefficient of k² the following eikonal equation:

$$(\nabla \theta {)}^2-n^2(\mathit{r})=0.$$(A5)

This first-order nonlinear partial differential equation can be transformed into a set of the following linear ordinary differential equations (Tatarskii, 1971; Kravtsov & Orlov, 1990)

$$\frac{\mathrm{d}\mathit{r}}{\mathrm{d}s}=\mathcal{l},$$(A6)

$$\frac{\mathrm{d}\mathcal{l}}{\mathrm{d}s}=\mathcal{l}\times \left[\left(\nabla \ln n\right)\times \mathcal{l}\right],$$(A7)

$$\frac{\mathrm{d}\theta }{\mathrm{d}s}=n.$$(A8)

Here s is the parametric variable, so that $\mathrm{d}s=\sqrt{{\sum }_{i=1}^3\mathrm{}(\mathrm{d}{r}_i{)}^2}$ is interpreted as the line element along the ray curve r = r(s) and ℓ = n⁻¹∇θ is the unit vector tangent to this curve. The eikonal function is obtained by solving equation (A8) and is given as a line integral

$$\theta =\int \mathrm{}n\left(\mathit{r}\right(s\left)\right)\mathrm{d}s$$(A9)

along the curve r(s). The ray curve, in turn, is determined by solving the equations (A6) and (A7). The following discussion is devoted to the perturbation approach to solving these equations for the case of wave propagation in a randomly inhomogeneous medium.

For a randomly inhomogeneous medium we have n(r) = 1 + δn(r) and equation (A7) can then be rewritten as

$$\frac{\mathrm{d}\mathcal{l}}{\mathrm{d}s}=\mathcal{l}\times \left[\left(\nabla \delta n\right)\times \mathcal{l}\right].$$(A10)

Defining the small parameter ν = 〈δn²〉 ≪ 1, we expand the variables r and ℓ in a series in powers of ν as

$$\mathit{r}={\mathit{r}}_0+\nu {\mathit{r}}_1+{\nu }^2{\mathit{r}}_2+\dots ,$$(A11)

$$\mathcal{l}={\mathcal{l}}_0+\nu {\mathcal{l}}_1+{\nu }^2{\mathcal{l}}_2+\dots$$(A12)

Substituting these expansions in the equations (A6) and (A10) we get

$$\begin{array}{cc}\frac{\mathrm{d}}{\mathrm{d}s}{\mathcal{l}}_0=0, & \frac{\mathrm{d}}{\mathrm{d}s}{\mathit{r}}_0={\mathcal{l}}_0,\\ \frac{\mathrm{d}}{\mathrm{d}s}{\mathcal{l}}_1={\mathcal{l}}_0\times \left[\nabla \frac{\delta n}{\nu }\times {\mathcal{l}}_0\right], & \frac{\mathrm{d}}{\mathrm{d}s}{\mathit{r}}_1={\mathcal{l}}_1,\\ \dots & \end{array}$$(A13)

It follows from the first set of equations that the vector ℓ₀ does not depend on s and that r₀ = sℓ₀. We use this observation by solving the second set of equations in equation (A13) together with the first-order approximation for the argument of the function δn, i.e., δn(r) ≈ δn(sℓ₀), to obtain

$${\mathcal{l}}_1={\mathcal{l}}_0\times \left[\left(\frac{1}{\nu }\underset{0}{\overset{s}{\int }}\mathrm{}\nabla \delta n\left(s^{\mathrm{\prime }}{\mathcal{l}}_0\right)\mathrm{d}{s}^{\mathrm{\prime }}\right)\times {\mathcal{l}}_0\right],$$(A14)

$${\mathit{r}}_1={\mathcal{l}}_0\times \left[\left(\frac{1}{\nu }\underset{0}{\overset{s}{\int }}\mathrm{}(s-s^{\text{'}})\nabla \delta n\left(s^{\text{'}}{\mathcal{l}}_0\right)\mathrm{d}{s}^{\text{'}}\right)\times {\mathcal{l}}_0\right].$$(A15)

In these equations and in the following we use the following shorthand notation for refractive index gradients:

$$\nabla \delta n\left(s{\mathcal{l}}_0\right)\equiv {\nabla }_{\mathit{r}} \delta n\left(\mathit{r}\right){|}_{\mathit{r}\equiv s{\mathcal{l}}_0},$$(A16)

i.e., the gradient field is calculated first and then the gradient magnitude and direction are selected at the point r = sℓ₀.

The solutions for ℓ and r can be obtained from the equations (A14) and (A15) as

$$\mathcal{l}\left(s\right)={\mathcal{l}}_0+{\mathcal{l}}_0\times \left[\left(\underset{0}{\overset{s}{\int }}\mathrm{}\nabla \delta n\left(s^{\text{'}}{\mathcal{l}}_0\right)\mathrm{d}{s}^{\text{'}}\right)\times {\mathcal{l}}_0\right]+\dots ={\mathcal{l}}_0+\underset{0}{\overset{s}{\int }}\mathrm{}{\nabla }_{\perp }\delta n\left(s^{\text{'}}{\mathcal{l}}_0\right)\mathrm{d}{s}^{\text{'}}+\dots ={\mathcal{l}}_0+\delta {\mathcal{l}}_1\left(s\right)+\dots ,$$(A17)

$$\mathit{r}\left(s\right)=s{\mathcal{l}}_0+{\mathcal{l}}_0\times \left[\left(\underset{0}{\overset{s}{\int }}\mathrm{}(s-s^{\text{'}})\nabla \delta n\left(s^{\text{'}}{\mathcal{l}}_0\right)\mathrm{d}{s}^{\text{'}}\right)\times {\mathcal{l}}_0\right]+\dots =s{\mathcal{l}}_0+\underset{0}{\overset{s}{\int }}\mathrm{}(s-s^{\text{'}}){\nabla }_{\perp }\delta n\left(s^{\text{'}}{\mathcal{l}}_0\right)\mathrm{d}{s}^{\text{'}}+\dots =s{\mathcal{l}}_0+\delta {\mathit{r}}_1\left(s\right)+\dots ,$$(A18)

where ∇_⊥ = ∇ − ℓ₀(ℓ₀ · ∇) is the transverse gradient to the initial direction ℓ₀. For the simplification of equations (A17) and (A18) we used the vector algebra relation a × (b × c) = (a · c)b − (a · b)c together with the property⁶ |ℓ₀| = 1. Equation (A18) shows that the ray leaving the point r = 0 can be represented in the first approximation as a vector pointing in the direction of the unit vector ℓ₀. The second approximation gives a certain shift δr₁(s) transverse to the direction of ℓ₀.

We truncate the expansion (A18) at the second term and substitute it in equation (A9) to obtain

$$\theta (s{\mathcal{l}}_0+\delta {\mathit{r}}_1(s\left)\right)=s+\underset{0}{\overset{s}{\int }}\mathrm{}\delta n(s^{\mathrm{\prime }}{\mathcal{l}}_0+\delta {\mathit{r}}_1(s^{\mathrm{\prime }}\left)\right) \mathrm{d}{s}^{\mathrm{\prime }}\approx s+\underset{0}{\overset{s}{\int }}\mathrm{}\delta n\left(s^{\mathrm{\prime }}{\mathcal{l}}_0\right) \mathrm{d}{s}^{\mathrm{\prime }}+\underset{0}{\overset{s}{\int }}\mathrm{}\nabla \delta n\left(s^{\mathrm{\prime }}{\mathcal{l}}_0\right)·\delta {\mathit{r}}_1\left(s^{\mathrm{\prime }}\right) \mathrm{d}{s}^{\mathrm{\prime }}.$$(A19)

This expression includes the deviation of the ray from the straight line, since the spatial argument has a small component δr₁ transverse to the initial direction ℓ. If ones is interested in the eikonal determined along the straight line sℓ₀, i.e., in the initial direction of propagation, one can write

$$\theta \left(s{\mathcal{l}}_0\right)=\theta \left(\right[s{\mathcal{l}}_0+\delta {\mathit{r}}_1\left(s\right)]-\delta {\mathit{r}}_1(s\left)\right)\approx \theta (s{\mathcal{l}}_0+\delta {\mathit{r}}_1(s\left)\right)-\delta {\mathit{r}}_1\left(s\right)·\nabla \theta (s{\mathcal{l}}_0+\delta {\mathit{r}}_1(s\left)\right)=s+\underset{0}{\overset{s}{\int }}\mathrm{}\delta n\left(s^{\mathrm{\prime }}{\mathcal{l}}_0\right) \mathrm{d}{s}^{\mathrm{\prime }}+\underset{0}{\overset{s}{\int }}\mathrm{}\nabla \delta n\left(s^{\mathrm{\prime }}{\mathcal{l}}_0\right)·\delta {\mathit{r}}_1\left(s^{\mathrm{\prime }}\right) \mathrm{d}{s}^{\mathrm{\prime }}-\delta {\mathcal{l}}_1\left(s\right)·\delta {\mathit{r}}_1\left(s\right).$$(A20)

Consider the last two summands on the right hand side of (A20) and denote them as δθ₂. Then, using equations (A17), (A18) and the relevant vector algebra relations, we obtain

$$\delta {\theta }_2\left(s{\mathcal{l}}_0\right)=\left\{\underset{0}{\overset{s}{\int }}\mathrm{}\mathrm{d}{s}_1\underset{0}{\overset{s_1}{\int }}\mathrm{}\mathrm{d}{s}_2(s_1-s_2)-\underset{0}{\overset{s}{\int }}\mathrm{}\mathrm{d}{s}_1\underset{0}{\overset{s}{\int }}\mathrm{}\mathrm{d}{s}_2(s-s_2)\right\}\{\nabla \delta n\left(s_1{\mathcal{l}}_0\right)·\nabla \delta n\left(s_2{\mathcal{l}}_0\right)-\left({\mathcal{l}}_0·\nabla \delta n\left(s_1{\mathcal{l}}_0\right)\right)\left({\mathcal{l}}_0·\nabla \delta n\left(s_2{\mathcal{l}}_0\right)\right)\}=\left\{\underset{0}{\overset{s}{\int }}\mathrm{}\mathrm{d}{s}_1\underset{0}{\overset{s_1}{\int }}\mathrm{}\mathrm{d}{s}_2(s_1-s_2)-\underset{0}{\overset{s}{\int }}\mathrm{}\mathrm{d}{s}_1\underset{0}{\overset{s}{\int }}\mathrm{}\mathrm{d}{s}_2(s-s_2)\right\}{\nabla }_{\perp }\delta n\left(s_1{\mathcal{l}}_0\right)·{\nabla }_{\perp }\delta n\left(s_2{\mathcal{l}}_0\right).$$(A21)

An alternatively useful formula for δθ₂ as the function of the transverse displacement δr₁ can be written as

$$\delta {\theta }_2=\underset{0}{\overset{s}{\int }}\mathrm{}{\nabla }_{\perp }\delta n\left(s^{\text{'}}{\mathcal{l}}_0\right)·\delta {\mathit{r}}_1\left(s^{\text{'}}\right) \mathrm{d}{s}^{\text{'}}-\underset{0}{\overset{s}{\int }}\mathrm{}{\nabla }_{\perp }\delta n\left(s^{\text{'}}{\mathcal{l}}_0\right) \mathrm{d}{s}^{\text{'}}·\delta {\mathit{r}}_1\left(s\right),$$(A22)

where we have used the fact that δr is perpendicular to the initial direction of propagation ℓ, and hence the scalar product of it with the gradient ∇δn reduces to ∇_⊥δn · δr. Relating the eikonal (A20) to the phase gained by propagation through the irregular medium of thickness s, we obtain

$$\phi =\mathrm{k\theta }\left(s{\mathcal{l}}_0\right)=\mathrm{ks}+k\underset{0}{\overset{s}{\int }}\mathrm{}\delta n\left(s^{\text{'}}{\mathcal{l}}_0\right) \mathrm{d}{s}^{\text{'}}+k\underset{0}{\overset{s}{\int }}\mathrm{}{\nabla }_{\perp }\delta n\left(s^{\text{'}}{\mathcal{l}}_0\right)·\delta {\mathit{r}}_1\left(s^{\text{'}}\right) \mathrm{d}{s}^{\text{'}}-k\underset{0}{\overset{s}{\int }}\mathrm{}{\nabla }_{\perp }\delta n\left(s^{\text{'}}{\mathcal{l}}_0\right) \mathrm{d}{s}^{\text{'}}·\delta {\mathit{r}}_1\left(s\right).$$(A23)

Since the first term in this expression is constant, it does not contribute to the calculation of either phase or amplitude scintillation indices.

Appendix B: Averaged Eikonal contribution 〈δθ₂〉

Let us calculate the mean value of δθ₂. We first consider the correlation tensor of the refractive index gradients, R_∇δn, whose components are related to the correlation function of the refractive index fluctuations, R_δn, as follows

$${\left[{\mathit{R}}_{\nabla \delta n}({\mathit{r}}^{\mathrm{\prime }},{\mathit{r}}^{\mathrm{″}})\right]}_{i,j}=\langle {\nabla }_i\delta n\left({\mathit{r}}^{\mathrm{\prime }}\right){\nabla }_j\delta n\left({\mathit{r}}^{\mathrm{″}}\right)\rangle =\frac{{\partial }^2}{\partial r_i^{\mathrm{\prime }}\partial r_j^{\mathrm{\text{'}\text{'}}}}{R}_{\delta n}({\mathit{r}}^{\mathrm{\prime }},{\mathit{r}}^{\mathrm{″}}).$$(B1)

We also use the notation ρ = r′ − r″ = (ξ₁, ξ₂, ξ₃)^T and assume that the medium is locally homogeneous, so that the correlation function R_δn depends only on ρ = |ρ|. To work with new variables we use

$$\begin{array}{cc}\frac{\partial \rho }{\partial r_i^{\mathrm{\prime }}}=\frac{{\xi }_i}{\rho }={\mathcal{l}}_i^0& \frac{\partial \rho }{\partial r_i^{\mathrm{\text{'}\text{'}}}}=-\frac{{\xi }_i}{\rho }=-{\mathcal{l}}_i^0,\end{array}$$(B2)

where ${\mathcal{l}}_i^0$ is the ith component of the unit vector ℓ₀ = ρ/ρ. Using the chain rule and noting that ρdρ = ∑_iξ_idξ_i we get

$${\left[{\mathit{R}}_{\nabla \delta n}\left(\rho \right)\right]}_{i,j}=-\left[{\delta }_{\mathrm{ij}}\frac{1}{\rho }\frac{\mathrm{d}}{\mathrm{d}\rho }+{\mathcal{l}}_i^0{\mathcal{l}}_j^0\left(\frac{{\mathrm{d}}^2}{\mathrm{d}{\rho }^2}-\frac{1}{\rho }\frac{\mathrm{d}}{\mathrm{d}\rho }\right)\right]R_{\delta n}\left(\rho \right).$$(B3)

For the averaged quantity of interest, we find

$$\langle {\nabla }_{\perp }\delta n\left(s_1{\mathcal{l}}_0\right)·{\nabla }_{\perp }\delta n\left(s_2{\mathcal{l}}_0\right)\rangle =\sum _{i,j}\mathrm{}({\delta }_{\mathrm{ij}}-{\mathcal{l}}_i^0{\mathcal{l}}_j^0){\left[{\mathit{R}}_{\nabla \delta n}\left(\rho \right)\right]}_{i,j},$$(B4)

where ρ = ℓ₀(s₁ − s₂) such that ρ = s₁ − s₂. Substituting equation (B3) into (B4) we get (Chernov, 1960; Tatarskii, 1971)

$$\langle {\nabla }_{\perp }\delta n_1{\nabla }_{\perp }\delta n_2\rangle =-\sum _{i,j}\mathrm{}({\delta }_{i,j}-{\mathcal{l}}_i^0{\mathcal{l}}_j^0{)}^2\frac{1}{\rho }\frac{\mathrm{d}}{\mathrm{d}\rho }{R}_{\delta n}(\rho )=-2\frac{1}{\rho }\frac{\mathrm{d}}{\mathrm{d}\rho }{R}_{\delta n}(\rho ).$$(B5)

Finally, to emphasize that the derived relation is valid not only for spatially homogeneous random media, but also for the more general case of locally homogeneous random media, we rewrite equation (B5) in terms of the refractive index structure function

$$D_{\delta n}\left({\mathit{r}}_1,{\mathit{r}}_2\right)=\langle {\left[\delta n\left({\mathit{r}}_1\right)-\delta n\left({\mathit{r}}_2\right)\right]}^2\rangle =2\left[R_{\delta n}\left({\mathit{r}}_1,{\mathit{r}}_1\right)-R_{\delta n}\left({\mathit{r}}_1,{\mathit{r}}_2\right)\right]$$(B6)

as

$$\langle {\nabla }_{\perp }\delta n_1{\nabla }_{\perp }\delta n_2\rangle =\frac{D_{\delta n}^{\mathrm{\prime }}\left(\rho \right)}{\rho },$$(B7)

where prime denotes the derivative.

The mean value of the eikonal contribution can then be calculated from equation (A21) as

$$\langle \delta {\theta }_2\left(s{\mathcal{l}}_0\right)\rangle =\underset{0}{\overset{s}{\int }}\mathrm{}\mathrm{d}{s}_1\underset{0}{\overset{s_1}{\int }}\mathrm{}\mathrm{d}{s}_2{D}_{\delta n}^{\mathrm{\prime }}(s_1-s_2)-\underset{0}{\overset{s}{\int }}\mathrm{}\mathrm{d}{s}_1\underset{0}{\overset{s}{\int }}\mathrm{}\mathrm{d}{s}_2\frac{(s-s_2)}{(s_1-s_2)}{D}_{\delta n}^{\mathrm{\prime }}(s_1-s_2)=\underset{0}{\overset{s}{\int }}\mathrm{}\mathrm{d}{s}_1\underset{0}{\overset{s_1}{\int }}\mathrm{}\mathrm{d}{s}_2\left[1-\frac{(s-s_2)}{(s_1-s_2)}\right]D_{\delta n}^{\mathrm{\prime }}(s_1-s_2)-\underset{0}{\overset{s}{\int }}\mathrm{}\mathrm{d}{s}_1\underset{s_1}{\overset{s}{\int }}\mathrm{}\mathrm{d}{s}_2\frac{(s-s_2)}{(s_1-s_2)}{D}_{\delta n}^{\mathrm{\prime }}(s_1-s_2).$$(B8)

Changing the order of integration and integrating one obtains (Keller, 1960; Tatarskii, 1971; Mahalov & McDaniel, 2019)

$$\langle \delta {\theta }_2\rangle =\underset{0}{\overset{s}{\int }}\mathrm{}\left[\left(s-s_1\right)-\frac{1}{2}\frac{\left(s^2-s_1^2\right)}{s_1}\right]D_{\delta n}^{\mathrm{\prime }}\left(s_1\right)\mathrm{d}{s}_1-\frac{1}{2}\underset{0}{\overset{s}{\int }}\mathrm{}(s-s_1{)}^2\frac{D_{\delta n}^{\mathrm{\prime }}\left(s_1\right)}{s_1}\mathrm{d}{s}_1=-\underset{0}{\overset{s}{\int }}\mathrm{}(s-s_1{)}^2\frac{D_{\delta n}^{\mathrm{\prime }}\left(s_1\right)}{s_1}\mathrm{d}{s}_1.$$(B9)

The major contribution to the last integral comes from the values s₁ → 0. Thus, one can approximate this integral as

$$\langle \delta {\theta }_2\rangle \approx -s^2\underset{0}{\overset{s}{\int }}\mathrm{}\frac{D_{\delta n}^{\mathrm{\prime }}\left(s_1\right)}{s_1}\mathrm{d}{s}_1.$$(B10)

For r > r₀, where r₀ is the characteristic correlation radius for the refractive index fluctuations, the function $D_{\delta n}^{\text{'}}\left(r\right)$ rapidly approaches zero. Since in the thin random medium model the medium thickness s is of the order of r₀, we can extend the upper limit of integration in equation (B10) to infinity without loss of generality, cf. (Keller, 1960; Tatarskii, 1971). In this case one can relate the mean eikonal increment 〈δθ₂〉 to the variance of the angle-of-arrival fluctuations ${\sigma }_{\alpha }^2$ as (Keller, 1960; Tatarskii, 1971)

$$\langle \delta {\theta }_2\rangle =-s{\sigma }_{\alpha }^2,$$(B11)

where α is the angle of rotation between the wavefronts of the signal arriving at two closely spaced receiving detectors.

For the Shkarofsky spectral model (9), the angle-of-arrival fluctuations can be derived as

$$\begin{array}{cc}{\sigma }_{\alpha }^2=\langle \delta n^2\rangle \sqrt{\pi }\frac{\mathrm{\Gamma }\left(1-\frac{p}{2}\right)}{2^p\mathrm{\Gamma }\left(\frac{p-1}{2}\right)}{\left(\frac{{\kappa }_0}{{\kappa }_m}\right)}^{p-1}{\kappa }_{m}s,& p<2.\end{array}$$(B12)

The quantity ${\sigma }_{\alpha }^2$ and hence the mean phase component 〈δθ₂〉 are vanishingly small in the ionosphere under quite conditions when $\sqrt{\langle \delta n^2\rangle }\sim 10^{-6}-10^{-7}$ since κ₀ < κ_m and s ∼ 1/κ_m. Moreover, for the Kolmogorov type of the spectrum, i.e., for p = 5/3, it has been shown by Tatarskii (1971) that 〈δθ₂〉 is vanishingly small when the amplitude fluctuations at the receiver site are small, i.e. S₄ ≪ 1. We also note that in the view of equation (A22) the approximation 〈δθ₂〉 = 0, resulting from ${\sigma }_{\alpha }^2=0$, is formally equivalent to 〈∇_⊥δn · δr₁〉 = 〈∇_⊥δn〉 · 〈δr₁〉 = 0, i.e., to the assumption that ∇_⊥δn and δr₁ are statistically independent random variables with zero means.

Appendix C: Contributions of irregularities at different heights to the diffractive scintillation and to the refractive tilt errors

This appendix discusses the difference between the refractive phase scintillation associated with the tilt error and the diffractive scintillation. In the case of diffractive scintillation the major contribution to the scintillation levels comes from the scattering on irregularities of the Fresnel scale ℓ_F. This scale is determined by the condition that the diffraction angle λ/ℓ_F coincides with the value ℓ_F/L, i.e., the angular diameter of the irregularity for the receiver position, see the Figure 9a. Here λ is the wavelength and L is the height of the scattering irregularity. The Fresnel scale is then ${\mathcal{l}}_F=\sqrt{\lambda L}$. The irregularity height L can be determined from the amplitude scintillation spectra, namely by identifying the Fresnel spatial frequency κ_F = 2π/ℓ_F as the roll-off frequency in the spectrum. Such a method gives, for example, an estimate of 350–400 km for L for L-band signals, which roughly corresponds to the ionospheric F2 peak heights (Kintner et al., 2007; Muella et al., 2013).

Figure C1 Geometric constructions for scattering on the Fresnel-scale irregularities (a) and on the irregularities of larger scales contributing to the wavefront tilt (b), (c). Se the text for more explanations.

Now let us consider the refractive phase fluctuations associated with the tilt of the wavefront caused by scattering on an irregularity of size ℓ > ℓ_F. If we think of the random medium as a collection of blobs or irregularities of refractive index variations, then we can consider such structures as random lenses. This is a very simplified view, but it can give some insight into the refractive phase scintillation associated with the tilt error. Consider how a lens of diameter ℓ and the refractive index δn affects the propagation of an incoming bunch of signal wave rays. Here, the quantity δn is the stochastic component of the medium refractive index n. The focal length of this lens is given by F = ℓ/δn and depending on the sign of δn the lens tends to focus or defocus the incoming rays. For definiteness, let us choose δn > 0 so that the lens is a focusing one. For the ionosphere, assuming ℓ is in the range 500 m–1 km, the focusing distance is F ∼ 10⁷–10⁸ m, so the lens is a very weak one with a very long focal length.

We are interested in what happens when the lens is tilted at the randomly varying angle α ≪ 1 with respect to the case where the incoming signal ray direction coincides with the lens axis. The lens tilt is equivalent to the situation where the incoming wavefront is tilted with respect to the initial wavefront. Now consider two observation points A and B at a small distance d from each other, see Figures 9b and 9c. The phase difference between these two points δS, i.e., the random phase increment resulting from the tilt error, can be related to the tilt angle α as α = λδS/2πd. The mean square fluctuation of the angle α for the Shkarofsky spectral model for irregularities has already been given by equation (B12).

Suppose two receivers R1 and R2 are located at the distances L₁ and L₂ > L₁ from the lens. Experiencing the tilt, the rays are deflected from the initial direction to the distances d₁ and d₂ for R1 and R2, respectively. In terms of our previous notations the deflection distance d_i, i = 1, 2 is the abbreviated notation for δr₁ (L_i), cf. equation (A18). It is obvious that d_i = L_i d/F such that the ratio of the variances for deflection distances is

$$\frac{\langle ({\mathrm{d}}_2{)}^2\rangle }{\langle ({\mathrm{d}}_1{)}^2\rangle }={\left(\frac{L_2}{L_1}\right)}^2\frac{{\sigma }_{\alpha ,2}^2}{{\sigma }_{\alpha ,1}^2},$$(C1)

where ${\mathit{\sigma }}_{\mathit{\alpha },\mathit{i}}^{\mathbf{2}}$ is the angle-of-arrival variance with statistical parameters of the medium at hights L_i, i = 1, 2. In this way, we account for the more general situation where both receivers R1 and R2 have the same location, but sense the “lens” tilt effect from two different layers in the ionosphere. Assuming that the statistical properties of the medium do not change significantly with the altitude one can see from the relation (C1) that the upper layers of the ionosphere, say at the altitude L₂, have more influence on the phase fluctuations caused by the random wavefront tilt than the lower layers, say at the altitude L₁. This situation is in contrast to the diffractive type of scintillation where the scintillation producing layer is roughly coincident with the ionospheric F2 peak height and some additional remapping is required to use the Swarm data for scintillation modeling (Imam et al., 2023).

Appendix D: Gradient of electron density product NeGIX

In this appendix we give a brief summary of the new Swarm product called the electron density (Ne) Gradient indeX (NeGIX). It is based on the in-situ electron density measurements from the Swarm satellites Alpha and Charlie (Swarm A and C).

The estimation of the gradient components is performed by using the combination of the electron densities obtained between two satellites or using the values obtained by the same satellite. The electron density values in such combinations are distinguished by the notation Ne_j, where the index j refers to the spatial position P_j where the measurement was made at the time t_j. For the calculation of the NeGIX, the combination of points is used, for which the measurement times are taken within the fixed time window of 8 seconds. For the set of all measured values obtained within this time window one defines the value

$$\left|\nabla N{e}_{\mathrm{ij}}\right|=\frac{|N{e}_i-N{e}_j|}{\Delta s_{\mathrm{ij}}},$$(D1)

where Δs_ij is the spatial distance between the points P_i and P_j. The direction of the axis between the points P_i and P_j is called “dipole” by Jakowski & Hoque (2019). For N available measurements within the mentioned time window, the maximum number of possible dipoles is N_D = N(N − 1)/2.

The value |∇Ne_ij| is the absolute value of the electron density gradient for ionospheric structures of scale Δs_ij. Using the information about the relative position of points P_i and P_j we can reconstruct the gradient vector as

$$\nabla N{e}_{\mathrm{ij}}=\left|\nabla N{e}_{\mathrm{ij}}\right|(\sin {\delta }_{\mathrm{ij}} \cos {\delta }_{\mathrm{ij}}{)}^{\mathrm{T}}=(\nabla \mathrm{Ne}{x}_{\mathrm{ij}} \nabla \mathrm{Ne}{y}_{\mathrm{ij}}{)}^{\mathrm{T}},$$(D2)

where ∇Nex and ∇Ney are the meridional and zonal components of the gradient and δ_ij is the azimuth angle for the corresponding dipole.

For a set of dipoles for the specified time window we can define an average gradient value within that window as

$$\langle \nabla \mathrm{Ne}\rangle =\frac{1}{N_D}\sum _{\{i,j\}}\mathrm{}\nabla N{e}_{\mathrm{ij}}=(\langle \nabla \mathrm{Nex}\rangle \langle \nabla \mathrm{Ney}\rangle {)}^{\mathrm{T}},$$(D3)

where the summation is performed over all pairs of P_i, P_j that form N_D dipoles. The equation (D3) defines the NeGIX index for the Swarm. The conventional units of the NeGIX index are [cm⁻³ km⁻¹]. Figure D1 shows an example of the range of values of the NeGIX and its typical behavior on the day and night sides of the Earth.

Figure D1 The values of the total NeGIX index obtained during and after the geomagnetic storm of April 23–24, 2023 for the longitudinal sector correspond to those in Figures 4 and 5. Both the ascending (Earth’s day side) and descending (night side) parts of the Swarm A and C orbits were considered and each figure contains several passages of the satellites. The ascending passages on April 23 and the descending passages on April 24 took place under the most disturbed conditions (Kp index = 8). The NeGIX during these passages shows in general higher values.

References