Modeling of ionospheric scintillation

– A signal, such as from a GNSS satellite or microwave sounding system, propagating in the randomly inhomogeneous ionosphere, experiences chaotic modulations of its amplitude and phase. This effect is known as scintillation. This article reviews basic theoretical concepts and simulation strategies for modeling the scintillation phenomenon. We focused our attention primarily on the methods connected with the random phase screen model. For a weak scattering regime on random ionospheric irregularities, a single-phase screen model enables us to obtain the analytic expression for phase and intensity scintillation indices, as well as the statistical quantities characterizing the strength of scintillation-related fades and distortions. In the case of multiple scattering, the simulation with multiple phase screens becomes a handy tool for obtaining these indices. For both scattering regimes, the statistical properties of the ionospheric random medium play an important role in scintillation modeling and are discussed with an emphasis on related geometric aspects. As an illustration, the phase screen simulation approaches used in the global climatological scintillation model GISM are brie ﬂ y discussed.


Introduction
The scintillation phenomenon is familiar to everyone who has looked at the sky during cloudless nights.By looking for some period of time, one can observe the remarkable twinkling of stars.These random changes in the apparent star brightness or its position are due to the random distortion of light coming from a star caused by refractive index irregularities in a turbulent troposphere.The fluctuation in brightness or, in other words, in intensity is one example of the scintillation phenomenon.The pulsation of star images (star dancing) is another manifestation of this phenomenon known as angular scintillation.
If we were able to sense the electromagnetic radiation at radio wavelength, we would similarly perceive scintillation of distant radio sources of celestial origin and from artificial satellites.In this case, scintillation is caused by random scattering on irregularities of electron density of the electrons and protons in the Earth ionosphere, solar wind, or in the interstellar medium.As the electron density varies as a random function of location, the associated refractive index changes similarly.Hence, as a plane wavefront of the incoming signal passes through the region where the refractive index is low, the velocity of radio waves will be high, and the wavefront will advance farther.With a higher refractive index, the wavefront is delayed.
In combination, these effects will result in a corrugated wavefront and phase fluctuations.Some portions of the wavefront will be convex in the direction of propagation, and some of them will be concave.At some distance from the scattering region, the partial waves with such concave or convex wavefronts will interfere constructively or destructively, producing a complex intensity pattern.This article will focus on ionospheric scintillation phenomena, i.e., amplitude or intensity and phase fluctuations.Amplitude scintillation is more intense in equatorial regions, while phase scintillation is more frequent and stronger at high latitudes (Aarons, 1982;Liu et al., 2012;Tsai et al., 2017).Amplitude scintillation normally occurs when the Fresnel dimension of the propagating radio signal is of the order of irregularity scales in the ionosphere.This condition is met for the post-sunset equatorial irregularities formed in equatorial F-spread and the equatorial plasma bubbles (Moorthy et al., 1979;Patel et al., 2009;Yokoyama & Stolle, 2017).The phase scintillation is primarily caused by steep ionospheric density gradients and irregularities formed in the cusp and auroral precipitation regions and polar cap patches (Spogli et al., 2009;Jiao & Morton, 2015).Fluctuations in phase in these regions are additionally enhanced due to geometric effects, i.e., the relative position of field-aligned irregularities and the signal path (Forte & Radicella, 2004).The occurrence of ionospheric scintillation is also strongly time-dependent.It varies daily, seasonal, and yearly and is more frequent during post-sunset hours at low Topical Issue -Ionospheric plasma irregularities and their impact on radio systems latitudes (Basu et al., 1996;Huang et al., 2014a;Jin et al., 2018).During vernal and autumn equinoxes, there is a higher probability of scintillation (Muella et al., 2013;Jiao & Morton, 2015).Finally, the strength and occurrence of phase and intensity scintillation correlate with solar and geomagnetic activity at high latitudes (Abdu et al., 1998;Souto Fortes et al., 2015;Guo et al., 2017;Béniguel, 2019).
Intense phase and intensity fluctuations may have a strong negative impact on modern technical infrastructure, the robustness of some life-critical services, and the quality of scientific data acquired in remote-sensing measurement campaigns.Intensity scintillation causes the signals from a satellite to fade below the average level.If the level of fading exceeds some prescribed margin of a receiver, the signal becomes masked by the noise and cannot be reliably retrieved.Phase scintillation induces a frequency shift that is responsible for cycle slips in GNSS signals (Knight & Finn, 1998;Horvath & Crozier, 2007).When this shift exceeds the bandwidth of the phase lock loop, the satellite tracking might fail, and in some cases, the loss of lock with the satellite might follow (Conker et al., 2003;Aquino et al., 2005;Sreeja et al., 2011).For positioning, navigation, and timing (PNT) applications, such as precise point positioning (PPP), this effect can impair the functionality of this service, putting vital systems at risk.
In astronomical and remote sensing applications, scintillation may degrade the contrast and resolution of images as well as incorporate artifacts and increase the noise level in retrieved data.In radio astronomy, ionospheric and interplanetary scintillations lead to image wander, degradation of image contrast, pulse broadening of pulsar signals, angular broadening of point sources, and thus degradation of the resolving power of radio telescopes (Strom et al., 2001).In applications such as GNSS reflectometry, intensity scintillation has the greatest impact on the degradation of retrieval performance in altimetric and scatterometric applications (Camps et al., 2016).As the algorithms for the sounding of Earth and planetary atmospheres are usually based on ill-posed inversion problems (Leitinger, 2001), the introduction of additional sources of noise due to scintillation may considerably influence the performance of retrieval methods.In occultation measurements, strong refractive scintillation impacts the measured values of bending angle and contributes to the residual ionospheric error (Hinson, 1986;Carrano et al., 2011;Ludwig-Barbosa et al., 2020).The remarkable effect of intensity and phase scintillation on synthetic aperture radar (SAR) missions is the "azimuth streaking" notable as a spatial modulation in the position of optimal azimuth registration (Gray et al., 2000;Carrano et al., 2012;Sato et al., 2021).Other effects of scintillation that influence the performance of SAR are defocusing (Belcher, 2008;van de Kamp et al., 2009), signal decorrelation (Ishimaru et al., 1999), degradation of coherence properties Béniguel & Hamel (2011), and random phase perturbations (Rogers et al., 2014).
However, it should be noted that scintillation has not only negative effects.The knowledge of the scintillation level makes it possible to localize irregularities of the electron density in the ionosphere (Sokolovskiy, 2000), deduce the climatological trends for their occurrence (Siefring et al., 2011), retrieve the properties of ionospheric, interplanetary, and interstellar medium (Rickett, 1977;Bhattacharyya et al., 1992;Shishov et al., 2010), or even infer the structure of the local magnetic fields in planetary ionospheres (Hinson & Tyler, 1982;Hinson, 1984).In combination with other instrumentation, scintillation monitoring becomes a powerful tool for investigating the formation and evolution of ionospheric structures such as spread-F and equatorial plasma bubbles (Kelley et al., 2002).
Due to the aforementioned impact of scintillation on various important services and on remote sensing data products, the theoretical modeling and computer simulation of these phenomena were attracting the attention of researchers for more than seventy years (for review, see e.g.Ratcliffe, 1956;Gurvich & Tatarskii, 1975;Fante, 1975a;Crane, 1977;Yeh & Liu, 1982;Kravtsov, 1992;Priyadarshi, 2015).The first theoretical ideas on how to model scintillation were closely related to the emergent statistical theory of turbulence (Kolmogorov, 1941) and the development of the theory of random functions and corresponding methods of spectral analysis (see e.g.Rice, 1944Rice, , 1945)).The former theory gave an understanding of how one can describe complex random systems such as turbulence or randomly inhomogeneous media, while the theory of random functions provided the rigorous mathematical apparatus for dealing with signals modulated randomly in phase and amplitude.The major theoretical milestones in such developments were connected with the need for a description of optical light propagation in tropospheric turbulence or sound waves in a turbulent ocean (Tatarski, 2016).To such landmarks belong the introduction of the concept of locally homogeneous random media (Obukhov, 1949), the development of Rytov's theory of smooth perturbations (Obukhov, 1953), the derivation of field moment equations (Shishov, 1968), and asymptotic methods for solving such equations (Gochelashvily & Shishov, 1975;Yakushkin, 1975;Fante, 1975b).These concepts allow one to establish a rigorous mathematical apparatus for dealing with wave propagation in a random medium and obtain formulas of scintillation indices for multiple cases important in practice.
In connection with the ionospheric propagation, the first theoretical attempts to explain scintillation were associated with the diffraction theory (Ratcliffe, 1956).The propagation of electromagnetic waves in a random ionosphere has been proposed to be modeled via a thin sheet of material that modulates the wave amplitude and phase.The transmitted wave is then represented as a set of uniform infinite plane waves traveling in different directions.The values of corresponding directional cosines are considered to be governed by the specific angular spectrum (Booker & Clemmow, 1950).At some distance from the modulating sheet of material, the interference of plane waves forms a specific diffraction pattern in intensity.If the modulation strength of the thin sheet is a random function of position, the phase shifts of plane waves emanating from it, as well as the formed diffraction pattern are random functions of the position (Booker et al., 1950).In these initial approaches, the randomness of scattering directions, i.e., of the directional cosines values, is the driver of random modulation for both phase and amplitude directly behind the screen.However, it was soon realized that for most cases of radio wave propagation in the ionosphere, the amplitude modulation factor could be neglected, and the spectrum of random phase modulation is sufficient to model scintillation (Hewish, 1951).This concept of the phase-modulated sheet, called simply the phase screen, has become an extremely handy tool for studying and modeling ionospheric phase and amplitude scintillation.
Among further advances and milestones in phase screen modeling of scintillation, one can distinguish the incorporation of spatial anisotropy of field-aligned ionospheric irregularities in phase screen generation (Briggs & Parkin, 1963;Singleton, 1970a), the introduction of multiple-phase screens for a description of extended random media (Fejer, 1953), and the replacement of Gaussian spectra for phase screen generation with more realistic ones with power-law dependence (Rufenach, 1972).Further, Uscinski (1985) has shown that there exists a close resemblance between the methods of moment equations for electromagnetic wave amplitudes and the theory of multiple-phase screens.The popularity of phase screen simulations of scintillation was also growing due to advances in numerical mathematics.After developing the fast Fourier transform (FFT) algorithm (Cooley & Tukey, 1965) and the split-step algorithm for solving wave equations with a variable refractive index (Hardin et al., 1973), the FFT-based phase screen generation became spread out as the numerical method for both intensity and phase scintillation simulations.
The present article reviews the basic methods and theories used in the modeling of ionospheric scintillation.Instead of reviewing the ideas and results for particular scintillation models, we have the intention to provide the basic theoretical framework for describing and modeling electromagnetic field propagation in a random ionosphere that is applicable to the majority of the models listed in Table 1.The principal focus of this article is on the scintillation modeling with phase screens, namely the single-and multiple-phase screen approaches.Section 2 presents the theoretical methods for describing the propagation of electromagnetic waves in a disordered medium.
Here we derive the wave equation for an electromagnetic field propagating in the random ionospheric media.The solution to this equation can be obtained by effectively superposing the wave propagation in a vacuum together with the random modulation of its phase.Particularly, the statistical properties of this phase modulation are deduced from the power spectral densities of electron density fluctuations.Several popular models of such spectra are summarized in Section 3. Section 4 outlines the influence of geometric parameters of the communication link and scintillation-causing irregularities on the strength of the scintillation.We discuss in some detail the assumptions of flatand spherical geometries in the problem of scattering on random ionospheric inhomogeneities.Geometric parameters such as zenith and azimuth angles of the link, shape, and direction of anisotropic ionospheric irregularities and propagation factors are discussed.The analytic results for scintillation indices that include these effects are given in Section 5.In particular, the indices obtained in the flat-geometry approximation serve as a basis of the WBMOD scintillation model.The provided results in the spherical-geometry case can be useful for the consistent extension of the WBMOD applicability to the arbitrary slant communication links.The extension of these results to the case of multiple scattering is given in Section 6, where we discuss the FFT-based multiple phase screen (MPS) simulation algorithm.The Global Ionospheric Scintillation Model (GISM) is used to exemplify the scintillation simulation via the MPS.Finally, some conclusions have been drawn in Section 7.

Theoretical models of scintillation
When propagating in a random medium, an electromagnetic wave undergoes chaotic distortions of its phase and amplitude (or, alternatively, its power or intensity).In order to analyze, predict and effectively mitigate the degradation of the transmitted signal encoded in such a wave, the probability distribution functions (PDFs) for random phase and amplitude are needed.In general, the construction of such PDFs requires knowledge of all their moments or cumulants.In practice, it appears that moments up to the second-order are of importance and could alone serve as proxy parameters characterizing affects of random media on the propagating signal.This restriction of the interest to low-order moments is driven by two factors.Firstly, for a strongly irregular medium, one would expect that the PDF of interest approaches the Gaussian distribution as a consequence of the limiting theorems of probability theory (Papoulis & Pillai, 2002).In this case, knowledge of the first and second-order moments is sufficient to determine the total PDF.Secondly, the most important information is generally contained in low-order moments.For example, the mean-field provides information on the way that the incident wave is attenuated on passing through the medium.The second field moments are related to the visibility of fringes of two interfering signals, yielding information on the angular spread and power of the propagating signal.The second intensity moment (fourth field moment) is useful for the characterization of intensity fluctuations of the transmitted wave.
For this reason, the following statistical moments and their combinations are of primary interest in studying the effects of a randomly inhomogeneous ionosphere on the propagation of radio and microwave signals: where E is the field amplitude, I is the field intensity, S is the phase departure of the wave detected by the observer.The statistical averaging (ensemble-or time-averaging) here and in the following is denoted as h...i.Conventionally, the oneminute time averages are used for the scintillation metrics defined by equations ( 1) and (2).The proxy quantities S 4 and r S are called the amplitude (or intensity) and phase scintillation indices, correspondingly.As can be seen from the definition, both indices provide information on how strong fluctuating intensity and phase differ from their regular or mean values.
Restricting our attention to the analysis of the second-order statistical moments (1), ( 2) is one of the simplifications made while treating scintillation phenomena from the theoretical perspective.In the following sections, we will make a number of further simplifications.We will neglect the bending of signal rays in the ionosphere and consider the depolarization effects caused by the propagation in the ionosphere to be negligibly small.For simplicity, the spherical shape of the Earth is assumed, but the flat Earth approximation is discussed in the appropriate context.The ionospheric plasma would be considered to be collisionless, its random variations would be treated by statistical means that again rely on second-order statistical quantities such as autocorrelation functions.The simplified model for autocorrelation properties of electron density fluctuations would be utilized for incorporating the anisotropy effects in scintillation phenomena due to scattering on field-aligned irregularities.The bulk of hte random ionosphere would be treated effectively as a set of thin random phase screens, whose temporal variation would be incorporated into the spatial variation with the help of the so-called Taylor hypothesis.In view of all these assumptions and simplifications, the resulting theoretical methods are quite capable of gaining insights into the nature of scintillation phenomena and to be used in computer simulations of scintillation fading levels.
In this section, we outline the theoretical foundations of signal wave propagation in a random ionospheric medium and discuss some approximations and assumptions used.In particular, we derive the wave equation for a signal wave propagating in a randomly inhomogeneous medium, which is the central equation in the forthcoming analysis.The solution of this equation is essential for the calculation of the scintillation indices, as will be shown in the following sections.

Wave propagation in random medium
Let us consider a plane wave propagating in z direction.This wave impinges a random medium that extends over the range 0 z Dz, cf. Figure 1.In this section, we restrict our consideration to the normal incidence case that will be generalized to the oblique incidence case in Section 4. The electromagnetic field of this wave we write as Eðr; tÞ ¼ E 0 e iðkzÀxtÞ ; z < 0; where r = (x y z) T is the three-dimensional spatial coordinate vector, x is the angular frequency, k = x/c is the free-space wave number, c is the speed of light.The quantities E 0 and f(r) are the time-independent field amplitudes for waves propagating in a vacuum and in a random medium, respectively.
D. Vasylyev et al.: J. Space Weather Space Clim. 2022, 12, 22 If e is the dielectric permittivity of the random medium, its general expression can be written as where hei is the mean dielectric permittivity and de is the dimensionless random modulation factor possessing the property hdei = 0.With the use of Taylor's hypothesis (Taylor, 1938) the time dependence of the dielectric permittivity can be incorporated into the spatial coordinate.Taylor argued that if the velocity fluctuations of a random medium are small compared to the mean flow velocity, the temporal variation at any fixed point in space can be considered as a passage of a "frozen" spatial structure past the observer at the mean flow velocity v drift .In this case, the time dependence can be incorporated in the spatial coordinate, i.e., e(r, t) = e(r À v drift t,0).For convenience, in the following, we will simply write the spatial dependence as e(r).In many practical situations, the Taylor hypothesis can be applied with relatively mild restrictions (Tennekes, 1975) and is suitable for analyzing ionospheric plasma turbulence (Mikkelsen & Pécseli, 1980;Dyrud et al., 2008;Rino, 2011).
For electromagnetic waves with frequencies much larger than the plasma frequency (i.e.larger than 1-1.5 MHz) propagating in the collisionless plasma, the background dielectric permittivity can be written as (Rawer, 1993) heðrÞi ¼ e 0 1 À e 2 hN e ðrÞi e 0 m x 2 where e 0 is the vacuum permitivity, x p is the plasma frequency, e is the elementary charge, N e is the electron density and m is the electron mass.From this equation follows that the presence of free electrons modifies the free-space dielectric permittivity and makes it dependent on the frequency of the propagating electromagnetic wave.In the following we use the refractive index nðrÞ ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffi eðrÞ=e 0 p instead of the dielectric permittivity.
For a weakly inhomogeneous random ionosphere (Rino, 2011) the depolarization effects can be neglected such that each polarization component f p (r), p = 1, 2 of the field f(r) can be treated separately.For simplicity in the following, we omit the index p associated with certain polarization components.If we assume that the signal wave has a beam-like form with a definite propagation direction, which in our case is along the z axis, we can separate the fast and slowly varying field components as f(r) = u(r)e ikz .The fast varying component e ikz will not be of much interest in the following.Indeed, this term is eliminated in the definition of the intensity scintillation index, and in phase measurements it is usually treated as the phase reference; thus, it is also not relevant for a calculation of the phase scintillation index.The wave equation for the slowly varying amplitude u(r) can then be written in the so-called paraxial approximation as 2ik ouðr ?; zÞ oz þ Á ?uðr ?; zÞ þ 2k 2 dnðr ?; zÞ uðr ?; zÞ ¼ 0; ð6Þ where D \ = o 2 /ox 2 + o 2 /oy 2 is the component of the Laplace operator transversal to the propagation direction and r \ = (x y) T is the corresponding transversal coordinate vector.In the paraxial approximation, we have neglected the term with second-order derivatives with respect to z since for the directed beam, it is much smaller than the term kou/oz.The fluctuating part of the refractive index in equation ( 6) is where r e % 2.818 Â 10 À15 m is the classical electron radius and k = 2p/k is the wavelength of the electromagnetic wave.
The wave equation ( 6) is the starting point for a derivation of statistical quantities defined in equations (1) and (2).In this connection, various approximative schemes have been developed for the solving of (6), such as the geometric optics approximation (Tatarskii, 1971;Wheelon, 2004), the method of smooth perturbations, or Rytov approximation (Clifford, 1978;Wheelon, 2003), the Huygen-Kirchhoff method (Ratcliffe, 1956;Banakh & Mironov, 1979), the method of field moments (Tatarskii, 1971;Liu et al., 1974;Prokhorov et al., 1975), and Figure 1.Geometry of propagation through a thin (a) and extended (b) random medium.These regimes correspond to a situation when the medium thickness is comparable to or larger than the typical scale of medium inhomogeneities.D. Vasylyev et al.: J. Space Weather Space Clim. 2022, 12, 22 the method of functional integrals (Tatarskii et al., 1992;Zavorotny et al., 1992), to name just a few.

Split-step approach and the concept of a random phase screen
The popular solution method of equation ( 6) is the so-called split-step algorithm (Hardin et al., 1973;Tappert & Hardin, 1975).It is based on the observation that the wave equation is a combination of the terms responsible for wave propagation in a vacuum [the first and second terms in the L.H.S of Eq. ( 6)] and in a random medium with the fluctuating refractive index dn [the first and the third terms in the L.H.S of Eq. ( 6)].Now, if the spatial scale along the z axis, over which the refractive index fluctuation changes considerably, is denoted by ', two situations might occur.If the thickness of the random medium Áz % ', cf. Figure 1a, the medium slab can be effectively replaced with a thin randomly phase changing screen.In the case of an extended random medium with Áz > ', cf. Figure 1b can be effectively represented as a set of thin phase screens placed in a vacuum.We refer to the first case as the weak scattering regime and the second case as the multiple scattering regime.The combination of the split-step method with the model of single or multiple phase screens is a popular tool to solve equation ( 7), which we discuss in some detail in Section 6.1.
In order to see why the concept of the random phase screen appears while solving equation ( 6), let us consider the propagation along the infinitesimally small distance dz ( ' in the medium.The distance is chosen in such a way that dn(r \ , z + dz) dn(r \ , z), i.e., the fluctuating part of the refractive index, does not change sufficiently.While propagating this distance, the field gains a phase increment k dn dz and can be written as uðr ?; z þ dzÞ ¼ uðr ?; zÞe ik dn dz % uðr ?; zÞ 1 þ ik dn dz À 1 2 ðk dn dzÞ 2 þ ::: : By forming the combination we see that it is exactly the third term in the L.H.S. of equation ( 6) up to the minus sign.In this way, the medium-related part of the derivative ou/oz cancels the third term in (6), and only the vacuum propagation wave equation is left over.We can thus interpret the paraxial equation as the propagation in a vacuum superimposed with the random phase modulating effect due to the spatially variable refractive index fluctuation dn.
The cumulative phase gain while propagating in a random medium of thickness Dz is The second-order statistical properties of random phase fluctuations are inferred from the autocorrelation function, which can be written for a spatially homogeneous medium as q dn ðr 1;? À r 2;? ; zÞdz; ð11Þ where for any (real) quantity X the autocorrelation is usually defined as The second line in equation ( 11) is obtained by changing the integration of coordinates according to z = z 0 À z 00 , 2Z = z 0 + z 00 and the subsequent integration over Z.The limits of the integration were extended to infinite values since the correlation function q dn decays quickly, as the argument exceeds the correlation length ', which in turn does not exceed the medium thickness Dz.
The relation of autocorrelation functions for phase and refractive index fluctuations (or, equivalently, of the electron density fluctuations) provided by equation ( 11) suggests that the knowledge of statistical properties of the random medium is useful for generating physically sound phase screens.It appears that for phase screen generation, it is more convenient to operate with the Fourier images of the autocorrelation functions.The latter quantities, known as the power spectral densities (PSDs), and underlying analytical models are discussed in the next section.

Statistical description of ionospheric random medium
The autocorrelation function for an arbitrary real quantity X, as defined in equation ( 12), has the Fourier image which is the three-dimensional PSD of process X.The dimensionality of the PSD refers to the three-dimensional vector of spatial frequencies j = (j x j y j z ) T .Examples of quantity X that are of interest in the context of wave propagation in the ionosphere are the electron density fluctuations, dN e , or, equivalently, the reflective index fluctuations, dn, the random phase increment due to the transmission through a medium slab, du, the phase departure of the wave, S, observed by a ground observer, and the fluctuations of the amplitude dE/hEi or of the log-amplitude v = log(1 + dE/hEi).In the following, we will substitute the corresponding quantity instead of X in the specified contexts.
In practical situations, not all of the information on all components of the vector j is usually available.For example, if two or more ground stations record the diffraction pattern of a signal from a satellite, the correlation analysis can be performed to infere the properties of ionospheric irregularities (Briggs, 1984;Costa et al., 1988;Kriegel et al., 2017).The statistical D. Vasylyev et al.: J. Space Weather Space Clim. 2022, 12, 22 quantity inferred is then defined in a plane rather than in the three-dimensional space.Similarly, the satellite or rocket measurements usually yield one-dimensional information along their flight track (Inhester et al., 1990).In all these situations, the full three-dimensional spectrum (13) is not accessible, but rather its two-or one-dimensional analogies can be inferred.One therefore defines the two-dimensional PSD as (Tatarskii, 1971) and the one-dimensional one as (Kovasznay et al., 1949) The last relation in equation ( 15) is valid if the quantity X is spatially isotropic, i.e. its variation does not depend on some preferable direction in space.
The isotropy condition tends to be very restrictive and cannot be fully met in the ionosphere, where the electron density irregularities tend to align along the geomagnetic lines of force.The spatial homogeneity condition used in equation ( 11) is less restrictive and means that the correlation functions depend on the relative distance between two spatial points, where the quantities of interest are recorded, rather than on the exact point coordinates.In general cases, this holds if we are interested in spatial separations smaller than the characteristic sizes of ionospheric irregularities.Since the medium thickness, Dz, is larger than such sizes, it is useful to introduce the concept of local homogeneity.This concept follows from the observation that although the quantity X(r) might not be spatially homogeneous, the difference X(r) À X(r 0 ) as the detrended process might possess this property.So instead of the autocorrelation functions, the locally spatially homogeneous process X is characterized by the structure function (Wheelon, 2004;Rino, 2011;Tatarski, 2016) In particular, by using equations ( 11) and ( 14) the following relationship can be established between the structure functions of phase and refractive index fluctuations: See equation (17) bottom of the page Thus, the local homogeneity of refractive index fluctuations implies the local homogeneity of random phase increments for the propagating wave.
In the following, we will primarily use the structure functions instead of the autocorrelation functions.In this way, we emphasize that the random medium under consideration is not restricted to being spatially homogeneous, and the obtained formulas are valid in less restricted cases of local homogeneity.Additionally, as will be discussed in Section 6.1, the phase structure function ( 17) is a convenient metric to check if the relevant parameters in a computer simulation of scintillation with phase screens have been properly chosen.
We now list the most widely used spectral models for modeling ionospheric irregularities.For now, we will concentrate merely on the functional dependence of these models and discuss the inclusion of anisotropy effects in the next section.

Gaussian model
The simplest idea is to assume that the refractive index (electron density) fluctuations are correlated according to Gaussian law with some characteristic correlation radius r 0 .The corresponding spectra are The spectra are properly normalized such that holds true.The Gaussian correlation function and corresponding spectral models were popular in the early stages of theoretical descriptions of scintillation in the ionosphere (Bramley, 1954;Bowhill, 1961;Briggs & Parkin, 1963;Liu et al., 1974) and had their validity for characterizing the random ionosphere under highly disturbed conditions (Fejer, 1953;Bramley, 1954;Mercier, 1962;Uscinski, 1978;Ishimaru, 1978;Priyadarshi, 2015).

Power-law models
After the demonstration by Rufenach (1972) that the ionospheric F region irregularities exhibit power-law spectra, multiple measurements have confirmed that this dependence is quite universal for various ionospheric conditions and geographic latitudes (Singleton, 1974;Crane, 1976;Umeki et al., 1977;Franke & Liu, 1983;Bhattacharyya & Rastogi, 1985;Basu et al., 1988a).The three-dimensional power-law PSD in the general form can be written as where Q[j] is a quadratic form of spatial frequency, p (3) is the three-dimensional spectral index, c is an optional scalar parameter, and f(j) is some function to assure the convergence of the integrated spectrum (Voitsekhovich, 1995).The quadratic form incorporates the spatial anisotropy of ionospheric irregularities placed in the geomagnetic field, the subject to be addressed in Section 4. The power law of type (20) signifies that there exists some nonlinear process that allows the inflow of the energy, e.g.due to solar heating or wind shears, to be redistributed between the inhomogeneities.This process is of a cascade type where the energy is redistributed among the different scales while larger inhomogeneities are split into smaller ones and are called the Richardson energy cascade. 5An illustration of such a process is the evolution of turbulent flows.In the developed turbulence, the energy-containing inhomogeneities of the largest size L 0 , referred to as the outer scale, consist of stretching and bending under the action of inertial forces that results in their splitting into less energetic and smaller inhomogeneities.This process continues until the inner scale ' 0 , i.e., the smallest inhomogeneity scale, is reached, the viscous forces start to dominate the inertial ones, and the energy is drained by dissipation.
Since the spectra ( 14) and ( 15) are of reduced dimensionality compared to (20), the corresponding two-and one-dimensional spectral indices would differ from the tree-dimensional one according to It follows that depending on the type of experiment used for obtaining the spectral characteristics of a random medium, the specific spectral index can be deduced, and the other two indices can be calculated by virtue of equation ( 21).In the following, we adopt the convention of using the one-dimensional spectral index and denote it as p = p (1) .This choice is dictated by the convenience in comparison with the vast majority of experimental data and results of a numerical simulation.However, one should keep in mind that only the three-dimensional spectra have a direct physical interpretation as they exhibit a true  Kintner & Seyler (1985) scaling behavior of the underlying nonlinear cascade process.Table 2 summarizes the typical values of spectral indices for particular ionospheric instability phenomena.

von Karman model
The popular spectral models for phase screen simulation are the von Karman PSDs where j 0 = 2p/L 0 is the spatial frequency associated with the outer scale L 0 , C(x) is the gamma function, and K n ðxÞ is the modified Bessel function of the second kind (Macdonald function).The normalization constants in ( 22) are obtained from equation ( 19).The corresponding structure functions for the refractive index and phase fluctuations are Expression ( 24) is obtained by substituting (22) in equation ( 17) and integrating.

Shkarofsky model
The von Karman spectra do not depend on the spatial frequency j m / 1=' 0 associated with the inner scale and, thus, do not account for the change of spectral properties in the highfrequency dissipation region, j > j m .The spectra proposed by Shkarofsky (1968) are useful if this behavior is of importance, see e.g.(Liu & Yeh, 1978;Knepp, 1983;Rino, 2011).The normalized Shkarofsky spectra are The corresponding structure functions of refractive index and phase fluctuations read as In deriving equation ( 28) the relation ( 17) has been used.

Two-component spectral models
The last spectral model we would like to mention is the twocomponent spectrum typical for complex systems with two cascade processes involved in forming the refractive index structures.Examples of such behavior are the direct and inverse cascades in the two-dimensional turbulence (Davidson, 2004).The two-component behavior of the spectrum might also be related to the formation of specific geometric forms of irregularities, such as the shocklike structures in equatorial bottom side spread-F (Hysell et al., 1994b).In any case, there exists an additional break at frequency j b such that j 0 < j b < j m .The two-component spectra have been studied by many authors, see e.g.Refs.Shishov (1974Shishov ( , 1976)); Carrano & Rino (2015), and here we derive the normalized spectra following the ideas of Shkarofsky (1968) as Here 2 F 1 (a, b; c; x) is the hypergeometric function and p 1 and p 2 are two spectral indices that correspond to spatial frequency regions j 0 < j j b and j b < j < j m , respectively.Table 2 provides some typical spectral indices for ionospheric processes characterized by two-component spectra.
It is worth noting that, in contrast to spectra in equation ( 25), the spectra of different dimensionality (29), (30), and (31) do not relate to each other by means of relations ( 14) and (15).Each spectrum is derived by assuming a certain dependence on the spatial frequency with exponents following from the dimensionality reasoning.The one-dimensional spectrum (31) was proposed by Shkarofsky (1968), and the two-dimensional spectrum (30) was used Franke et al. (1984), Franke & Liu (1985) for modeling scintillation in equatorial regions.The expressions for the corresponding structure functions are barely useful for practical calculations, and their expressions will be omitted.
To summarize, each spectral model considered above has its specific advantage from a practical perspective.The Gaussian spectrum is attractive since it allows one to reduce the computation complexity when dealing with the computation of scintillation indices, cf.Ref. Uscinski (1982).Moreover, as has been noted above, this spectrum is applicable to a highly disturbed random ionosphere.The Shkarofsky spectrum has universal applicability in the whole range of spatial frequencies of interest.The two-component spectrum can be considered as the generalization of the Shkarofsky model to the case when the electron density exhibits dependence on two spectral indices.In situations when the two-component spectra relevant for scintillation modeling might occur in structures associated with plasma bubbles (Basu et al., 1983;Muralikrishna et al., 2007;Yokoyama, 2017b), they are also connected with the nonlinear convection (Hysell et al., 1994a), reversed flow events Jin et al. (2019), and polar cap patches Spicher et al. (2015), to name just a few.Finally, the von Karman model is a handy alternative to the Shkarofsky spectrum in the phase screen computer simulations.The corresponding algorithms are based on the discretization of both the spatial domain, where the phase screen is defined, and the spatial frequency domain of the corresponding PSD.The dependence of the spectrum on the inner scale ' 0 , i.e., on the wave frequency j m , is not important if the grid spacing of the corresponding phase screen is chosen to be greater than ' 0 .The latter condition is usually met in the practice where the speed of computations is essential.
The discussed Gaussian and power-law spectra and the related structure functions are given in Figure 2. From Figure 2a, it is evident that the von Karman (22), Shkarofsky (25), and two-component spectra ( 29)-(31) behave very similarly at small spatial frequencies and start to deviate at the spectral break, j b , and at the inner scale, j m , frequencies.For a two-component spectrum, the change of the spectral slope p 1 in the range j 0 < j < j b to the slope p 2 in the range j b < j < j m is evident.Figure 2b shows the structure functions for phase and refractive index fluctuations for several spectral models.The general behavior of the structure functions is to saturate at spatial separations exceeding the largest irregularity scale given by the outer scale L 0 = 2p/j 0 .The interesting observation is that although the structure functions D dn behave differently for the models considered, the phase structure functions D du for the Gaussian, von Karman, and Shkarofsky spectra are almost identical.This result is quite surprising considering that the spectral characteristics of the Gaussian spectrum are quite distinct from those of the power-law spectra, cf. Figure 2a.Thus, although the specific spectral components of the Gaussian PSD contribute to fluctuations dn in a distinct way in comparison to the corresponding spectral components of the power-law PSDs, the statistical properties of the resulting phase fluctuations are almost identical.As a consequence, the simulation of the scintillation indices via the random phase screen technique, discussed in Section 6, yields similar results for power-law and Gaussian phase screens, provided that the correlation radius of the Gaussian PSD is set to be equal to the outer scale, i.e. when r 0 = L 0 = 2p/j 0 .
Figure 2. Power spectral densities (a) and the structure functions (b) for several models discussed in the text.The spectral indices are chosen to be: p = p 1 = 1.3 and p 2 ¼ 3:8.The one-, two-, and three-dimensional spectra are normalized such that hdn 2 i ¼ 1.The different power-law models show similar behavior at the outer scale frequency, j 0 , but deviate at the spectral break, j b , and the inner scale, j m , spatial frequencies.The correlation radius in the Gaussian spectrum is r 0 ¼ 2p=j 0 .D. Vasylyev et al.: J. Space Weather Space Clim. 2022, 12, 22 4 Geometric parameters in the modeling of scintillation The scintillation strength depends on the relative positions of the observer and sender as well as on the geometric form of scintillation-producing irregularities crossing their field of view (Parkin, 1967).Expressed in terms of local coordinates, the first aspect corresponds to the dependence of scintillation indices on the zenith and azimuth angles.Additionally, after scattering within the ionospheric layer, the signal propagates some distance in free space until being detected by the ground observer.This free space propagation should also be accounted for when calculating the scintillation.The second aspect is considered in the modeling of scintillation by incorporating the anisotropy of ionospheric irregularities into the power spectral density of the refractive index fluctuations.

Dependence on zenith and azimuth angles
To characterize the position of ionospheric irregularities in space, we use the zenith and the azimuth angles in the following considerations.We should differentiate these angles for some spatial point along the communication path where the incoming wave is scattered by irregularity and the corresponding variables of the ground observer.We denote the zenith and azimuth angles for the former case as h and /, correspondingly, while for the latter case, they are denoted as h o and / o .The line element ds along the slant path at the point with height z above the ground is related to the corresponding vertical length element dz as where h is the zenith angle with respect to the local vertical direction of the given point.If one neglects the finite curvature of the Earth or the finite curvature of the ionospheric shell, the zenith angle h would coincide with the observer's zenith angle h o for all points along the straight propagation path.
If the finiteness of the curvature is taken into account, the following relation can be established under the assumption of spherical geometry: Here R È is the Earth radius, z is the altitude of the given point.
The distinction between the zenith angles h and h o is important as some authors prefer to use the coordinates associated with the point along the propagation path, referring to equation (32) (Butler, 1951;Ellison & Seddon, 1952;Booker, 1958;Rino, 1979a), while some other authors prefer to perform recalculations according to equation (33) in order to get the scintillation index in terms of the observer's coordinates (Briggs & Parkin, 1963;Wand & Evans, 1975).The flat geometry approximation h = h o can be used under the condition of small values of zenith angles.
In the approximation of flat geometry, the zenith angle does not depend on the position along the propagation path.Thus, the integration of equation ( 32) for the slant path yields Z z 0 ds ¼ z sec h ¼ z sec h o .For the spherical geometry, the corresponding path length is obtained from equation ( 33) where the slant path correction factors expressed in terms of the dimensionless variable Depending on which coordinates are adopted, either the factor S(h o ) or the factor S 0 (h) should be used in the calculation of the slant path length.
In similar manner we should distinguish the azimuth angle of the observer, / o , and the corresponding angle at some point along the slant path.The relationship between these two quantities involves the geographic coordinates, namely the latitude U and the longitude K, of the corresponding points and reads as This equation suggests that in general the azimuth angles /, / o differ from each other and only in the case of U % U o and K % K o one can set / % / o .The usage of the coordinate sets {h, /} or {h o , / o } is dictated merely by practical convenience.For example, in the case of the GNSS applications when the signal is recorded by the ground receivers, the appropriate coordinates are h o and / o .For remote sensing applications with SAR satellites, the choice of h, / is more appropriate for studying the signal scattering in the ionospheric layer.In this case it is also convenient to express these variables in terms of the local coordinates of the transmitting satellite using relations similar to equations ( 33) and (36).

Dependence on geometric parameters of irregularities
The dependence of scintillation indices on zenith and azimuth angles is closely related to the geometrical parameters of the anisotropic scintillation-produced irregularities (Singleton, 1970b(Singleton, , 1973)).Moreover, as the irregularities are usually fieldaligned, the additional rotation on the geomagnetic dip and declination angles should be performed for proper accounting of the relative position of a scintillation-causing irregularity and propagation path.
Some information on the geometrical structure of ionospheric irregularities can give, for example, the classical correlation analysis of diffraction patterns observed by spatially separated ground observers (Briggs, 1968).This analysis reveals that the diffraction patterns could be analyzed if one assumes that the surfaces of constant correlation of refractive index D. Vasylyev et al.: J. Space Weather Space Clim. 2022, 12, 22 fluctuations are the ellipsoids with the major semi-axes a, b, c, cf. Figure 3a.Following the Singleton convention (Singleton, 1970a), we use the coordinates r, s, t such that the minor ellipsoid semi-axis c is along the t-axis and the semi-axes a and b are along the axes s and r respectively, cf. Figure 3b.Finally, it is convenient to associate the minor semi-axis with the reference correlation radius r 0 .Then other semi-axes are expressed through r 0 via the dimensionless scaling factors a, b ! 1.Specifically, the equation (Briggs & Parkin, 1963;Singleton, 1970a;Moorcroft & Arima, 1972) defines the surface of constant correlation, i.e. the surface on which q dn and D dn attain the constant values.Alternatively, the surfaces (37) with a ?1/a, b ?1/b, r 0 ?1/r 0 represent the surfaces in the reciprocal domain where U dn attain constant values.
In the spatial domain, it is convenient to relate the Singleton coordinates r, s, t to those in the conventional coordinate system x, y, z.If we choose x, y, z to be a NED (North, East, Down) coordinate system, the Singleton and NED coordinates can be related by a series of proper rotations, as shown in Figure 3.The rotation angles are: the magnetic declination d and inclination (dip) w angles as well as the tilt angle c that characterizes the tilt of a bunch of magnetic field lines relative to the observation plane.These geometric considerations result in the following definition of the quadratic form of spatial frequencies that enters the general definition of three-dimensional PSD in equation ( 20) (Rino & Fremouw, 1977;Rino, 1979aRino, , 1982aRino, , 2011)): Here b C is a symmetrical matrix with the following elements: In the limiting case d = 0 these formulas coincide with the results of Rino & Fremouw (1977).
If the statistical properties of phase fluctuations are of interest, for example, in connection with the phase screen modeling, then, in the view of definition (17), the quadratic form to be inserted in the spectrum U dn is This result is valid if the propagation direction coincides with the z axis of the NED coordinate sysem, i.e., for zero zenith angle.For small zenith angles, the flat geometry approximation can be used, i.e., one can assume that the ionospheric layer is plan-parallel and the Earth's surface is approximately flat.In this case, the corrections to equation ( 45) due to the finiteness of the zenith angle can be incorporated by substituting j z ¼ À tan ha k;? Á j ?as has been shown in Rino & Fremouw (1977).The quadratic form can be written as where are the coefficients of the quadratic form expressed in terms of the observer azimuth / and zenith h as long as the approximation (32) is applicable.The vector a k;? ¼ ðcos / sin /Þ T is the unit vector along the transverse direction to the wave vector k such that a k;? Á j ?¼ j x cos / þ j y sin /.The relevant geometrical aspects of the model ( 46) are shown in Figure 4a.We also note that the assumption of the flat geometry is equivalent to the assumption that h = h o as discussed in Section 4.1.Some improvement of this approximation could serve the hybrid model, in which one uses the quadratic form (46) but defines the dependence between the zenith angles of the observer and the point along the slant path via the relation (36).Such modification allows one to avoid the divergencies associated with sec h at large zenith angles.
Finally, the following quadratic form can be derived for an arbitrary value of the zenith angle: where Here the local coordinate system x; y; z (North, East, Down) originates at height z above the observation plane XY.The coordinate system x 0 ; y 0 ; z 0 is chosen in such a way that the z 0 axis is pointed along the wave vector k while the x 0 and z 0 axes lie in the plane formed by k and the down axis z.In turn, the direction of the wave vector is defined in terms of the zenith angle h and the azimuth angle /.The models (a) and (b) also assume the flat and spherical geometries, respectively.This results in different definitions of the slant ranges and transformation rules for h, / into the observer's zenith, h o , and azimuth / o , angles (see the text for more details).
D. Vasylyev et al.: J. Space Weather Space Clim. 2022, 12, 22 Here the spatial frequencies j 0 x , j 0 y are reciprocal coordinates of the coordinates x 0 , y 0 defined in the plane transversal to the propagation direction, i.e., to the wave vector k.The model ( 49) assumes the spherical symmetry of the Earth and ionospheric shell.This assumption, to which we refer as the spherical geometry approximation, is not biased by the projection transformation of the kind used in ( 46), and, thus, is consistent for arbitrary values of zenith angles.

Propagation factors
So far, we have considered the propagation through a random ionospheric medium of thickness Dz, cf.equation ( 3).In practical situations, there is a finite distance from the layer of such medium and the observation plane placed at z > Dz in the coordinate system chosen, cf.Figures 1a and 1b.The propagation over such a distance can be considered to happen in a vacuum.The propagation factors take this effect into account (Clifford, 1978;Tatarski, 2016).
Here we consider two models corresponding to equations ( 46) and ( 48) with the note that the obtained results for the zero zenith angle are the same as obtained within the model ( 45).We also consider a weak scattering case that corresponds to placing a single phase screen at z = 0.For the spectrum based on the model ( 46) the correlation functions for the phase departure, S, and for its log-amplitude, v = log(1 + dE/hEi), can be obtained as

See equation (50) bottom of the page
where R \ = (X Y) T is the vector defined in the observer plane such that R ?¼ r ?À z tan ha k;? , cf. Figure 4 and In a similar way for the spectral model based on equation ( 48) we obtain

See equation (52) bottom of the page
where S 0 ðhÞ is given by equation ( 35) and the variable j 0 ? is defined in the plane transversal to the propagation direction while R ? is transversal to the observer's vertical direction, see Figure 4.Here b T is the transformation matrix from R ? to the coordinates that are transversal to the propagation direction and have the origin at the ground.The explicit form of b T is not relevant for further discussions.
The propagation factors in equations ( 50) and ( 52) for the case of plane waves propagating through a thin phase changing screen read as (Yeh & Liu, 1982) In the so-called approximation of smooth perturbations or Rytov approximation the propagation factors are defined as These results have been obtained for plane wave propagation and the extension to the spherical wave can be found in Tatarskii (1971).
Inspecting the functional dependencies ( 53) and ( 54) as met in equations ( 50) and ( 52), one can see that the log amplitude propagation factor serves as the low-pass filter that inhibits the spatial frequencies higher than the Fresnel frequencies

1=
ffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi kz sec h p or 1= ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi kz S 0 ðhÞ q depending on the model.This is the reason why P v is sometimes referred to as the Fresnel filter.The action of the function has no profound effect while multiplying with the power-law spectra and in approximation of weak scattering one can set P S % 1.The phase scintillation index is readily obtained from (50) or from ( 52) via the relation That is the variance of the phase fluctuations observed on the ground, which is approximately equal to the variance of phase variations in the medium as if its thickness is replaced with the propagation distance z from the scattering layer to the ground.Another useful approximation in the limit of weak scattering is the relation of the scintillation index S 4 to the log amplitude variance.One obtains where the factor of four appears since the intensity scintillation index is expressed through the fourth power of the amplitude.This result is known as a thin phase screen approximation (Rino, 1979a) or the Rytov approximation (Obukhov, 1953;Yeh & Liu, 1982;Wheelon, 2003;Tatarski, 2016) depending on the form of the used propagation factor.Specifically, the propagation factor (53) corresponds to the thin-screen model, whereas the factor (54) corresponds to Vasylyev et al.: J. Space Weather Space Clim. 2022, 12, 22 the Rytov approximation.For r v < 0:4 this formula can be extended to the form S 2 4 % expð4r 2 v Þ À 1 that follows from the observation that the intensity statistics in this regime are governed by the log-normal distribution (Clifford & Yura, 1974).

Scintillation indices: weak scattering
For the weak scattering regime a single phase screen model can be applied for the calculation of scintillation indices according to equations ( 55) and ( 56).The relative simplicity of this case allows one to obtain these indices in analytical form.In this section we summarize these results for the models corresponding to equations ( 46) and ( 48).
If we use the von Karman spectrum ( 22) the first model will yield where A, B, C are defined in ( 47).Here the factor ab appears due to the normalization requirement of the PSD function and the parameter is called the strength of the scattering structure or simply the structure constant of the irregular medium.

Flat geometry approximation
We first give the results for the flat geometry approximation elaborated by Rino (1979a), see also Figure 4a.Substituting equation ( 57) in ( 50) and performing the integration we obtain the following phase scintillation index in the approximation of small zenith angles: where is called the geometric term (Rino, 1979a;Fremouw, 1980a;Rino, 1982a).The geometric factor describes the enhancement of scintillation due to the alignment of the propagation path along the major axis of field-aligned irregularities.This enhancement effect is notable in auroral regions and primarily affects the signals from polar orbiting satellites with a lesser effect on the signals from GNSS satellites (Forte & Radicella, 2004).We also note that the product of the structure function with the typical altitude of the scattering layer, i.e., the expression C s z, is used as a convenient characteristic of how highly the ionospheric layer is irregular.Thus, this factor, called the height integrated strength, modulates the strength of the scintillation (Fremouw & Larsinger, 1981;Strangeways et al., 2013;Mannix et al., 2017).
The result (59) can be extended to the case when the twocomponent spectrum ( 29) with included anisotropy for irregularities is used instead of (57).In this case, the phase scintillation index is found to be r where 2 F 1 (a, b; c; x) is the hypergeometric function of Gauss.By deriving this expression, it has been assumed that the anisotropy coefficients for spatial frequencies j > j b are proportional to the corresponding coefficients for the interval j 0 j j b .This is equivalent to the assumption that the axial ratios a and b are equal in both intervals, but the correlation lengths r 0;1 r 0 and r 0;2 might be different.Under such an assumption, the scintillation index (61) turns out to be independent of r 0;2 .In the limit of one component spectrum, p 1 ¼ p 2 , j b ¼ j 0 , equation ( 61) reduces to the result (59).We now insert equations ( 57) and ( 53) in equation ( 50) and use the definition (55) to obtain the intensity scintillation index in the limit of weak scattering.In order to simplify the integration, we use the observation for the ionospheric scintillation j 2 0 Z ( 1, where Z is defined in equation ( 51).As a consequence, the outer scale wavelength j 0 can be set to zero (Rino, 1979a).With this simplification, the intensity scintillation index is obtained as Here the combined geometry and propagation factor is where P k ðxÞ is the Legendre function of order k and a 1 ¼ sin h cos /, a 2 ¼ sin h sin / are the components of the unit vector a k;? .By derivation of equation ( 62) it has been assumed that 0 < p < 4. We also note that equation ( 62) is fully equivalent to the formulas originally derived in Rino (1979a) but is given in a more simplified and symmetrical form.The analytic results (62) and ( 59) have been used as a basis to create the climatological scintillation model WBMOD (Secan et al., 1995(Secan et al., , 1997)).
The physical meaning of equation ( 62) has been given by Rino (1979a) as follows.The strength of the amplitude scintillation is controlled by the product C s Á ð1= ffiffiffi Z p Þ Àp , i.e. it is proportional to the variance of the refractive index fluctuations according to equation ( 58) and to the factor related to the effective Fresnel radius ffiffiffi Z p .The latter contribution can be interpreted as a one-dimensional power spectral density of ionospheric irregularity evaluated at the spatial wavenumber corresponding to the Fresnel radius.D. Vasylyev et al.: J. Space Weather Space Clim. 2022, 12, 22

Extension to a spherical geometry case
For large zenith angles, the effects connected with the finite curvature of the ionospheric layer become important.Figure 4b shows the proper scattering geometry for this case.As we have mentioned already, the model ( 48) equipped with the appropriate definition of the local coordinate system along the communication path, cf.equations ( 33) and (36), is more suitable for a consistent calculation of scintillation indices.We first observe that the scintillation indices are obtained from the corresponding autocorrelation functions with the spatial argument set to zero.Comparing the definitions (50) and ( 52) for R ?¼ 0 reveals that the expression ( 59) is valid for arbitrary zenith and azimuth angles if one replaces the parameters A, B, C with the corresponding primed parameters defined in equations ( 49).The same is applicable for the calculation of S 4 with the exception that in the primed coordinates the function h defined in ( 51) is simplified to hðj 0 ?Þ ¼ j 2 ?and one should use the variable Z 0 instead of Z.This modification results in the following expression for the intensity scintillation index: The difference in expressions for the propagation factors J and J 0 relies on the absence of the specific projection transformation connected with the flat geometry assumption inherent to the model (62).Figures 5a and 5b show the geometric terms for the flatand spherical-geometry models, i.e. the factors G and G 0 defined in equations ( 60) and (64), respectively.Both geometric factors account for geometric enhancement connected with the purely geometric effects of the relative position of the propagation path and the ionospheric irregularity.If the propagation path is aligned along the major semi-axis of the surface of constant autocorrelation, the geometric factors exhibit maximal values.One can observe the shift of the enhancement region for G 0 at large values of the observer zenith angle h o .This is because G 0 in contrast to G accounts for the discrepancy between h o and the zenith angle h at the scattering point.The difference in the levels of enhancement factors in Figures 5a and 5b is connected to the difference in the definition of slant ranges for flatand spherical-geometry models.When substituting the geometric factors in the corresponding expressions for the phase scintillation indices, this difference is washed out as seen on Figures 5c and 5d .This is because the factor G, while being substituted in equation ( 59), is multiplied with the factor sec h o , which grows fast for large zenith angles.The net effect of this multiplication results in comparable phase scintillation values for the flatand spherical-geometry models.
Figures 6a and 6b show the scintillation indices corresponding to flatand spherical-geometry models for the zero azimuth angle at the scattering point, i.e., for / = 0°.Figures 6a and 6b show the corresponding models for the case / = 90°.At large zenith angles, one observes the shift of the scintillation maximum due to the discussed geometric enhancement effect.Additionally, one can observe that the scintillation indices calculated within the flat geometry approximation tend to diverge at h ?90°.This happens due to the divergent character of the secant function at large values of the argument.As expected, both models yield similar values of scintillation indices for small zenith angles.We also note the interpretation of specific values of the dip angle w and their relation to the geographic position of the observer in Figures 5 and 6 should be performed with the proper model of the geomagnetic field.
To our best knowledge, the results obtained for the spherical geometry case are new.This extension has a number of important applications.For example, in some reflectometric applications, the signal is sent under small elevation angles relative to the reflecting point on the Earth surface.In this situation, the analysis of scintillation-associated fades of the sounding signal requires taking into account the effects due to the finite curvature of the Earth.Similarly, the spherical geometry approximation is suitable for limb sounding and occultation measurements for the analysis of the scintillation-related noise contributions.It should be noted, however, that the obtained formulas ignore the refractive ray bending of the propagating signal.This effect can still be incorporated in the above formalism by replacing the observer zenith angle h o with the apparent zenith angle (Weisbrod & Anderson, 1959;Mahan, 1962), which value is determined by the refractive properties of the ionosphere.We omit any details on these calculations here and note that the ray tracing algorithms incorporated in scintillation models such as the GISM is responsible for an accurate accounting of the refractive ray bending effects.

Scintillation indices: multiple scattering
In the previous section, we considered the model of a single thin phase screen.This model makes it possible to describe the weak scattering regime and corresponds to the case of the random medium whose size is comparable to the typical correlation length of irregularities, cf. Figure 1a.The idea of a phase screen simulation of scintillation phenomena can be adopted as well to the extended random medium, cf. Figure 1b, where the medium thickness exceeds the correlation length of ionospheric inhomogeneities.The propagation in an extended medium is accompanied by multiple scattering and is effectively modeled by propagation through a set of thin phase screens.The current section provides an overview of the fast Fourier transform (FFT) assisted multiple screen simulation method.

Multiple phase screens approach
In this section, we consider the method of multiple phase screens for simulating electromagnetic wave propagation through a random medium.As we have already mentioned, the crucial part of this method is to generate random phase screens such that their statistical properties would resemble those of the random refracting medium.We therefore start with the power spectral density for refractive index fluctuations U dn ðj x ; j y ; j z Þ.This spectrum depends on the individual onedimensional wave numbers j x ¼ 2p=x, j y ¼ 2p=y, j z ¼ 2p=z with x, y, z being associated with the typical scale sizes of the medium in the x, y, and z direction, correspondingly.We choose the z axis to lie along the propagation path direction and denote the transversal wave number as usually6 with j ?¼ ðj x j y Þ T and the corresponding spatial vector as r ? .As we have already discussed in Section 4, in such a coordinate system it is sufficient to know the function U dn ðj x ; j y ; j z ¼ 0Þ with an appropriate model of the irregularity anisotropy to generate the phase screen at the point z ¼ 0. By means of the definition (17), we can use the two-dimensional spectral density for this task by virtue of the relation F dn ðj ?; zÞ ¼ 2pdðzÞU dn ðj ?; j z ¼ 0Þ: As the next step we define the two-dimensional power spectral density for phase fluctuations as where dz is some distance along the propagation path that is smaller or equal to the thickness Áz of the random medium slab.The particular models of the spectrum F dn are given by equations ( 22), ( 25) and can be easily extended to the case of anisotropic irregularities by methods presented in Section 4. Equation ( 66) can now be used for generating the single random screen where the location is set to the position z s ¼ z þ dz=2, i.e. in the middle of the path between z and z þ dz.The screen plane is perpendicular to the propagation direction z and the transversal vectors j ?and r ?lie within this plane.It is a well-known property of random functions that the process governed by the specific power spectral density can be generated through the linear filtering of Gaussian white noise with the square root of the spectrum, followed by an inverse Fourier transform (Rice, 1945;Jeruchim et al., 2000).The reason why the square root appears is that the power spectral density is connected with the signal power, i.e. the second power of the signal.Thus, we can generate a continuous random phase screen realization using the relationship where h G ðj ?Þ is a hermitian complex Gaussian process with zero mean and unit variance.
In numerical simulations the integral in ( 67) is approximated by the summation duðx; y; z s Þ ¼ X jx X jy ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi F du ðj x ; j y ; z s Þ q hðj x ; j y Þe iðjxxþjy yÞ Áj x Áj y : The variables x and y are now considered to be at points determined by the M Â N grid with spacings Áx and Áy, respectively.It is clearly that x ¼ mÁx, y ¼ nÁy with m, n being the integer indices.It is a usual convention to take the center of coordinates in the center of the phase screen grid, such that ÀM=2 m M=2 and ÀN =2 n N =2.Similarly, we choose the variables j x and j y to be considered as coordinates of points on a complementary grid, such that j x ¼ m 0 Áj x , j y ¼ n 0 Áj y and m 0 , n 0 are integer indices.
The discrete noise process hðj x ; j y Þ in equation ( 68) can be expressed in terms of the discrete version of the process h G of equation (67) as where is the discrete complex Gaussian noise with zero mean and unit variance and g r , g i its real and imaginary components.Due to the hermiteness of h, we have the following for the complex conjugate value: h Ã G ðm 0 ; n 0 Þ ¼ h G ðÀm 0 ; Àn 0 Þ.The scaling factor 1= ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi Áj x Áj y p is introduced so that the correlation function of the process h corresponds to a discrete two-dimensional delta function.In view of all the limitations connected with the generation of true random numbers in a computer simulation, we are still going to refer to the generated phase screen (68) as a random one.We now rewrite equation ( 68) such that the summation runs over the indices m 0 and n 0 .If the sizes of the screen along the x and y axes are L x and L y , respectively, we can write for the grid spacings: Áx ¼ L x =M and Áy ¼ L y =N .It is convenient in the following to rescale the spatial frequencies as f ¼ j=2p and discretize the frequencies as f x ¼ m 0 Áf x , f y ¼ n 0 Áf y .As the definition (68) involves the discrete Fourier transform, it is optimal to use the fast Fourier algorithm and choose M and N both to be powers of two.Moreover, using the Nyquist sampling theorem we chose the spacings in the reciprocal domain as Áf x ¼ 1=L x and Áf y ¼ 1=L y .Performing the appropriate change of coordinates in equation ( 68) we obtain Here instead of using the coordinate z s , describing the position of the screen, we incorporated index s in the definition of the phase.This index now refers to the specific phase screen, since as the next step we will consider that we have a set of l phase screens located at positions z s , s ¼ 1; 2; :::; l.An example of a generated random phase realization together with the corresponding PSD is given in Figure 7.
The anisotropy of the ionospheric irregularities is exhibited in the elongation of the power spectral density along the j x coordinate, cf. Figure 7a.The corresponding random phase screen in Figure 7b exhibits patterns elongated in the complimentary direction.Figure 8a compares the phase structure function obtained from the (isotropic) random phase screen simulation and the theoretical one given by equation ( 24).The criterion of such good agreement of the simulation with the theory serves as a verification that the phase screen parameters are chosen appropriately, there are no aliasing artifacts, etc.
We now turn to the multiple phase screen (MPS) simulation algorithm (Knepp, 1983;Grimault, 1998;Knepp & Nickisch, 2009).The idea of MPS simulations is based on the split-step method proposed in Refs.(Hardin et al., 1973;Tappert & Hardin, 1975) and mentioned already in Section 2. The basic idea is that the solution of the paraxial wave equation ( 6) can be modeled as a result of propagation over a set of successive regions filled with the vacuum.These free-space regions are considered to be separated by random phase changing screens in such a way that each phase screen redefines the initial conditions for the vacuum propagation in the region between this screen and the next sequential screen.The solution of the free-space wave equation is obtained by solving the corresponding equation in the Fourier domain and the subsequent application of the inverse Fourier transformation into the spatial domain.The FFT algorithm is used to perform these Fourier transformations.
Let us look at this procedure in some detail.We chose the location of the source at ðx 0 ; y 0 ; z 0 Þ and again align the z axis along the propagation direction.The position of the jth phase screen we specify to be at z j þ dz j =2 with dz j ¼ z jþ1 À z j and j ¼ 1; 2; :::; l, with l being the number of phase screens.The electromagnetic wave propagation is given by the paraxial equation ( 6) which we split in a vacuum propagation part and  Vasylyev et al.: J. Space Weather Space Clim. 2022, 12, 22 the part corresponding to propagation in the infinitely thin random phase screen: See equations (72), ( 73) bottom of the page where z s À z 0 is the distance from the receiver to the transmitter and indices "vac" and "med" refer to the vacuum and infinitely thin medium propagation parts, respectively.The geometry of the MPS simulation is summarized in Figure 8b.We now solve the wave equation ( 72) for propagation in a vacuum with z 2 ðz j þ dz j =2; z jþ1 þ dz jþ1 =2Þ in order to find the field at the end of this interval where ũmed is the Fourier image of the field amplitude formed at the ith phase screen.Here we have also used the identity . Similar expressions to (74) are obtained for the end intervals ½z 0 ; z 1 þ dz 1 =2Þ and ðz l þ dz l =2; z r with appropriate expressions for ũmed .The boundary condition for the ith screen, ũmed , in equation ( 74) is found in turn by taking the Fourier transform of the solution of equation ( 73) at z ¼ z j þ dz j =2, that is See equation (75) bottom of the page where in the last line we have used the assumption that the screen is infinitely thin.The solutions ( 74) and ( 75) show that the field at z jþ1 þ dz jþ1 =2 is determined by the field formed by the phase screen at z j þ dz j =2 propagating the distance ðdz j þ dz jþ1 Þ=2 in the vacuum behind this screen.Similarly, the field in z j þ dz j =2 is determined by the field generated by the ðj À 1Þth phase screen propagating the distance ðdz jÀ1 þ dz j Þ=2 in the vacuum.Thus the field at the given location can be traced backward by such a split-step procedure up to the source of the electromagnetic wave.These steps in solving the paraxial wave equation can be easily rewritten in the discrete form when we define the variables r ?and j ? on the grids defined above.Equation (74) then reads as where ũmed ðm 0 ; n 0 ; jÞ ¼ FFTðu vac ðm; n; jÞ exp½iduðm; n; jÞÞ ð77Þ follows from equation ( 75).Here FFT and FFT À1 are operations of direct and inverse fast Fourier transform, respectively.The phase fluctuation defined on the grid for i-th phase screen, duðm; n; jÞ, is given by equation (71).By iterating the procedure of the field calculation given in equations ( 76), (77) over 2ik all phase screens, the transmitted field at the receiver site z r can be finally obtained for the specified initial field at z 0 .
The moments of the electromagnetic field required for intensity scintillation index calculations, and the variance of phase fluctuations at the receiver, i.e., phase scintillation index, are obtained by the following procedure.Firstly, the propagation through the set of phase screens is simulated multiple times.For each simulation run the phase screens are generated from scratch.Secondly, the values of transmitted amplitude and phase are collected and used to construct such quantities as the intensity, squared intensity, and squared phase.Finally, the obtained quantities are then averaged over an ensemble of simulation runs and the corresponding scintillation indices are calculated.
We conclude this section by noting that the outlined method of phase screen generation based on the FFT algorithm is not the only possible approach.This method is most commonly used and allows one to generate phase screens with given statistical properties.One should point out that the accuracy and consistency of the FFT method have a close relationship with the number of points in the reciprocal spatial domain at which the PSD is sampled.The insufficient number of sample points may result in the introduction of unphysical artifacts and aliasing effects.Some strategies and algorithms for choosing the optimal sample points, grid spacings, or forms of computational meshes are discussed in detail in Knepp (1983), Lane et al. (1992), Schmidt (2010), Burckel & Gray (2013), Jia et al. (2015).They are based primarily on the use of the abovementioned Nyquist criterion, on the proper estimate of the receiver region of interest, or on the addition of low-frequency information using subharmonics.The correct reconstruction of the phase structure function could serve as a marker that the sampling grid is chosen properly, cf. Figure 8a.In this respect, the alternative phase screen generated approaches that do not require additional checks might be of interest for the simulation of ionospheric scintillation.We mention here the Zernike polynomial-based phase screen generation method (Roddier, 1990;Carbillet & Riccardi, 2010), the random midpoint displacement-based fractal method (Zhai et al., 2015), the multiscale method (Beghi et al., 2011) or the method based on local principal component analysis (Beghi et al., 2013).

Example of multiple phase screen simulations: the GISM model
As an illustration of the MPS algorithm, we will discuss the Global Ionospheric Scintillation Model (GISM) (Béniguel, 2002;Béniguel & Hamel, 2011;Béniguel, 2019).The background ionospheric model for the GISM is the NeQuick 2 model (Hochegger et al., 2000;Nava et al., 2008), which provides the value of electron density for a specified spatial point in the ionosphere.The refractive effects due to propagation in the ionosphere are accounted for by the ray-tracing algorithm incorporated in the GISM.Thus, in general cases, the propagation path is not a straight line but rather a complex ray curve consisting of multiple segments.The number of segments and the length of each segment are determined adaptively by the ray-tracing algorithm on the basis of the relative position of the transmitter and the receiver.The phase screens are positioned along such a signal ray so that the screens are concentrated at the maximum of the F layer slant profile determined along the ray.The number, l, and spacings dz j , between phase screens are also determined adaptively.
The generation of phase screens in the GISM is based on the von Karman spectrum, F dn for refractive index fluctuations, cf.equation ( 22).The strength of fluctuations in the spectrum for jth screen, i.e. the quantity hdn 2 i j ¼ r 2 e k 2 k À2 hdN 2 e i j , is determined by assuming that the electron density has a constant level of fluctuations ffiffiffiffiffiffiffiffiffiffiffiffiffi ffi hdN 2 e i j q =hN e i j ¼ 0:05.This fluctuation level is assumed to be the same for all phase screens, while the mean electron density, hN e i, varies in accordance with the NeQuick model.
The multiple phase screen simulation within the GISM is organized according to the procedure discussed in detail in Section 6.1.In situations when the anisotropy of irregularities does not sufficiently influence scintillation strength, the multiple phase screens are generated along only one spatial coordinate.Such a simplification considerably saves computational resources and results in MPS with one-dimensional phase screens.In this case, the random phase is generated on a onedimensional grid, such that it can be written as a one-dimensional discrete Fourier transform where L is the screen size and the one-dimensional spectral function W du is related to the corresponding function for refractive index fluctuations (22) as W du ¼ k 2 dz 2 j W dn .Analogous equations are obtained for the field amplitudes u vac ðm; jÞ and ũmed ðm 0 ; jÞ.
For a specified communication scenario, e.g. for orbiting satellite-ground or satellite-satellite links, the GISM outputs a table with information on the coordinates of the transmitter, the RMS value of intensity, and phase scintillation indices, and their margins.The latter gives the information on the probability for an occurrence of calculated scintillation indices.We also mention that the current version of the GISM is optimized for a calculation of the scintillation indices in low latitude regions.Some examples of modeled scintillation indices obtained with the GISM are shown in Figure 9.In particular, Figure 9  The examples of scintillation sky-maps for the observer in Douala/Cameroon (04 03 0 N, 009 41 0 E) are shown in (c) and (d) for S 4 and r S , respectively.One can clearly see that the scintillation level increases with the increase of the zenith angle if the line of sight is toward the direction where the total electron content is high.

Conclusions
The theoretical methods for modeling ionospheric scintillation phenomena have reached the state of maturity during the last several decades.This progress is connected with the development of sophisticated numerical techniques based on the formalism of phase screens.A set of phase screens effectively describes the action of random ionospheric disturbances on propagating electromagnetic waves.Combining free space propagation together with the transmission through phase screens, one can solve the wave equation for wave amplitude as well as infer the corresponding random phase increments.Finally, the scintillation indices of phase and amplitude are obtained by accumulating statistics of simulated intensity and phase fluctuations of the transmitted wave.
The present article reviews these developments in some detail.The emphasis has been placed on the basic concepts, methods, and approximations of ionospheric scintillation modeling based on phase screen formalism.The most important concept of such modeling is that we can simulate the random ionospheric medium only in terms of statistical quantities, e.g. its structure functions and the power spectral densities for electron density (refractive index) fluctuations.The complexity of such a random medium is evident by inspecting the empirical spectra of electron density fluctuations.These spectra show the power-law dependence on spatial frequencies, i.e. quantities that are proportional to reciprocal scales of irregularities causing scintillation.The power-law signifies that the continuum of various scales is involved in the random modulation of amplitude and phase of the transmitting signal.Accounting for the scattering of signal waves at each such scale would be too complex, and this fact makes the statistical approach the most useful and central in scintillation modeling.The article reviews several popular models of power spectral densities such as Gaussian, von Karman, Shkarofsky, and two-component spectra.
As the ionospheric irregularities tend to align along the geomagnetic field lines, these spectra in general, exhibit anisotropy.D. Vasylyev et al.: J. Space Weather Space Clim. 2022, 12, 22 A standard approximation is that such field-aligned irregularities are characterized by an ellipsoidal form for a surface of constant autocorrelation of electron density fluctuations.For realistic scintillation modeling, it is important to know the parameter set of the specified spectral model, i.e., the axial ratios of the correlation ellipsoid, geometric parameters of its position relative to the communication link, the value of the spectral index, the strength of scattering structures, the outer scale of the medium, etc.We review the dependence of scintillation indices on these parameters.In particular, in the case of weak scattering, it is possible to establish such dependencies in analytic form.We have skipped, however, a detailed discussion of the morphology of ionospheric irregularity parameters.
The objective of discussing some key ideas in ionospheric scintillation modeling has led to the omission of several aspects.Scintillation phenomena should be treated in a specified operational context.This means that depending on the type of the technical system, measuring or detecting apparatus, infrastructure object of interest, etc., only specific aspects of scintillationrelated disturbances gain their importance.Depending on the type of particular infrastructure element, the measurement device, or signal postprocessing algorithms, the naturally occurring random fluctuations of phase and amplitude may be partially mitigated or even enhanced.For example, GNSS choke ring antennas are designed to optimally reduce multipath effects while tracking low-elevation satellite signals, while the postprocessing of signals from satellites performs the detrending of geometric factors due to orbital motion.In reflectometric applications at large zenith angles, the receiver detects the multipath contributions that might enhance the scintillation levels.In SAR imaging applications, the signal crosses the random ionospheric medium two times: while propagating in the forward direction as well as after the reflecting from the ground.Moreover, part of the signal can be backscattered and contribute to the ambient noise detected by an antenna.Thus, one expects the enhancement of scintillation due to the specifics of SAR observations.As we can see, all these applications have their own specifics and scintillation indices should be calculated consistently and in agreement with each use case or user requirement under consideration.An example of such case-specific scintillation products could serve the generation of synthetic scintillation data sets Humphreys et al. (2009); Ghafoori & Skone (2014); Rino et al. (2018).Still, these specific products are based on the theoretical foundations of scintillation modeling presented in this article.
Due to practical importance the scintillation forecast remains an important driver for modeling.The climatological models such as WBMOD or GISM are capable of roughly predicting scintillation levels.For example, short-term prediction can be achieved by observing the time series for the model input such as the F 10:7 solar flux, and predicting their future values, e.g., using the exponential smoothing methods Snyder et al. (2002), the autoregressive and moving average models (Box & Jenkins, 1976), the BATS and TBATS models (De Livera et al., 2012), to name a few.The methods of machine learning and artificial neuronal networks are useful for these demands as well (Qi & Zhang, 2008;Makridakis et al., 2018).
The climatological models tend to provide the averaged indices which might deviate from processes that are highly variable in time and space, which are responsible for scintillation occurrences such as those due to extreme space weather.
Therefore, the alternative forecasting approaches have been developed, e.g. based on the direct analysis and assimilation in the climatological models of recorded time series of scintillation indices and other geophysical parameters (Groves et al., 1997;Basu et al., 2002;Rezende et al., 2010;Redmon et al., 2010;McGranaghan et al., 2018) or on the usage of the background physics-based models for electron density (Nugent et al., 2021).In this connection, the space-based high-resolution sensing appears to be critical for accurate characterization and forecasting of scintillation conditions (Kelly et al., 2014).
Several steps might help us to establish reliable and accurate scintillation modeling, monitoring, and prediction services.With every year, more and more data is accumulated, allowing us to gain better insights into ionospheric conditions over solar cycles.The number of satellites and ground stations scanning the Earth ionosphere is also steadily increasing.The combination of data from different remote sensing sources, such as GNSS, SAR, and geostationary satellites, contribute to our better understanding of the volumetric spatial variability of the ionosphere.The development and deployment of the highresolution electron density monitoring instruments such as EISCAT3D, global ultraviolet imagers (GUVIs), and special sensor ultraviolet spectrographic imagers (SSUSIs), is attractive for an analysis of ionospheric irregularities with an improved spatial and temporal resolution.This, in turn, would facilitate our better understanding of the influence of ionospheric irregularities on radio-wave propagation in the disturbed ionosphere.We expect that the mentioned advances would facilitate the development of scintillation models comparable in their capabilities with the current meteorological weather prediction models.

Figure 3 .
Figure 3.The ellipsoidal surface with semi-axes a, b, c, is the surface of the constant autocorrelation function of refractive index fluctuations.For expressing the equation of the surface in local geographic coordinates, the transformation from North-East-Down coordinate system x; y; z to the Singleton coordinates r; s; t is used.The transformation consists of: (a) the rotation on the magnetic declination angle / followed by the rotation on the magnetic inclination or dip angle w; (b) the rotation on angle c that characterizes the tilt of the ionospheric layer relative to the observer plane.In order to visualize the orientation of the plane within which the geomagnetic field has constant magnitude the multiple field lines are shown in the sketch (b).

Figure 4 .
Figure 4. Cross-sections (hatched areas) of ellipsoids of constant autocorrelation of electron density fluctuations with the plane associated either with the slab of irregularities (a) or with the wave propagation direction (b).Here the local coordinate system x; y; z (North, East, Down) originates at height z above the observation plane XY.The coordinate system x 0 ; y 0 ; z 0 is chosen in such a way that the z 0 axis is pointed along the wave vector k while the x 0 and z 0 axes lie in the plane formed by k and the down axis z.In turn, the direction of the wave vector is defined in terms of the zenith angle h and the azimuth angle /.The models (a) and (b) also assume the flat and spherical geometries, respectively.This results in different definitions of the slant ranges and transformation rules for h, / into the observer's zenith, h o , and azimuth / o , angles (see the text for more details).

Figure 5 .
Figure 5. Geometric terms as the functions of zenith and magnetic inclination (dip) angles for: (a) flat geometry model, (b) spherical geometry model.The corresponding phase scintillation indices are shown on plots (c) and (d).The magnetic declination angle d, the azimuth angle /, and the tilt angle c are set to zero.The irregularity axial parameters are a ¼ 2, b ¼ 1:5.Other parameters are: k ¼ 30 cm, p ¼ 1:8, L 0 ¼ 5 km, z ¼ 350 km, hdn 2 i ¼ 1:3 Â 10 À13 .

Figure 6 .
Figure 6.Intensity scintillation indices for: (a) and (b) for the flatand spherical-geometry models and zero azimuth angle at scattering point, / ¼ 0 .Figures (c), (d) show the corresponding scintillation indices for / ¼ 90 .Two different values of the angle / are chosen in order to illustrate the variability of the scintillation index with azimuth.Other calculation parameters are the same as for Figure 5.

7.
Example of a generated two-dimensional random phase screen (b) using the anisotropic power spectral density F du (a).Pictures show only the central portion of the MPS grids in spatial and reciprocal spatial domains.The grid spacings of plots (a) and (b) are Áj x ¼ Áj y ¼ 2:5 Â 10 À3 m À1 and Áx ¼ Áy ¼ 9:7 m, respectively.The number of the shown grid vertices in plot (b) is reduced for a better visibility.The anisotropy coefficients are: A

Figure 8 .
Figure 8.(a) Simulated phase structure functions are compared with the theoretical von Karman models for isotropic and anisotropic random medium cases.The cross-section along one spatial dimension of this function is shown for convenience.(b) Geometry of multiple phase screens simulation method.

Figure 9 .
Figure 9. Examples of scintillation plots generated with the GISM.The total electron content values in TEC units calculated by the NeQuick2 background ionospheric model for 15 October 2015 at 20:00 UT (a).The corresponding post sunset distribution of vertical S 4 index shown on a globe for signal at the L band (b).The sky-maps as functions of zenith and azimuth angles for: (c) intensity scintillation index S 4 , (d) phase scintillation index r S .The location of the ground observer is indicated by a star on (b).
(a)  shows the global map of the total electron content at 20:00 UT on 15th October 2015.The global distribution of the post-sunset intensity scintillation index is shown in (b).

Table 2 .
Diversity of one-dimensional spectral indices for typical nonlinear processes and instabilities in the ionosphere.The abbreviations HL and E refer to high-latitude and equatorial regions, respectively.The ranges of indices, p 1 and p 2 , of two component spectra are given as (p 1 ) (p 2 ).