Internal Gravity Wave Turbulence in the Earth’s Ionospheric F-Layer

1. Introduction

Acoustic-gravity waves (AGWs) play crucial roles in characterizing ionospheric plasma dynamics. They are essential for accurately predicting atmospheric behaviors under varying meteorological conditions in the Earth’s atmosphere, the stratified layers of the solar atmosphere, and the magnetospheric environments of other planets [1,2,3,4,5,6]. A low-frequency branch of AGWs is typically known as internal gravity waves (IGWs), which have frequencies in the range ${10}^{-4}{\mathrm{s}}^{-1}<$$\omega$$<1.7\times {10}^{-2}{\mathrm{s}}^{-1}$ and wavelengths of approximately 10 km. IGWs propagate only in stable stratified fluids, when the typical Brunt–Väisälä frequency ${\omega }_g$ is positive. They play a vital role in transporting wave energy to charged particles in the Earth’s atmosphere and ionosphere. The nonlinear evolution of these waves and their interactions with fluid density perturbations are fundamental to both geophysical and applied research in atmospheric convection and the generation of atmospheric turbulence [7].

Nonlinear interactions between finite-amplitude AGWs in the Earth’s atmosphere lead to both resonant and nonresonant energy transfer between high- and low-frequency components of the wave spectrum, with resonant processes playing a significant role in rotating atmospheric conditions [8]. Analytical studies [9] exploring nonlinear couplings among acoustic-gravity modes with different wavelengths suggest the development of vortical structures. Gardner [10] described the vertical wavenumber spectrum as a scale-independent diffusivity that accounts for damping due to molecular viscosity, turbulence, and off-resonant wave-wave interactions. Similarly, a physically based gravity wave source spectrum parameterization for convection was implemented in a Global Circulation Model by Beres et al. [11]. Shaikh et al. [7] investigated the turbulent behavior of two-dimensional nonlinearly interacting IGWs and showed that the governing nonlinear equations for vorticity and density fluctuations in a nonuniform fluid lead to a dual cascade, resulting in an energy spectrum indicative of Kolmogorov-type turbulence.

The role of charged particles in the ionosphere and their effect on nonlinear solitary vortices associated with planetary Rossby waves was first systematically explored by Kaladze and Tsamalashvili [12], and later expanded upon in the works of Kaladze [13,14] and Kaladze et al. [15]. A detailed analysis of how these charged particles affect the linear and nonlinear propagation characteristics of acoustic gravity waves in the weakly ionized and electrically conductive D, E, and F layers of the ionosphere while accounting for the Coriolis force was provided in [16]. The latter proposed an extended two-dimensional model for the evolution of IGW-driven solitary vortices, including the effect of Pedersen conductivity, and concluded that Pedersen conductivity dissipates the wave energy due to the Ohm losses. The work also discussed how the Hall conductivity and the Coriolis force due to rotating fluids are irrelevant, but Pedersen conductivity can not be in the ionospheric F-layer. The theory was recently advanced in [17] employing analytical and numerical approaches to explore different forms of dipolar vortex structures and their time-dependent evolution, demonstrating how Pedersen conductivity and system nonlinearities influence their dynamics.

In this work, we aim to present a detailed numerical investigation of the turbulent properties of nonlinearly interacting IGWs in the ionospheric F-layer by the effects of Pedersen conductivity, thereby advancing further the nonlinear theory of IGW-driven solitary vortices as in Refs. [7,16,17]. Although Shaikh et al. [7] presented a similar study without the effects of Pedersen conductivity, the results may be relevant in the ionospheric D-layer. Here, we consider the Earth’s F-layer to analyze the turbulent spectrum resulting from dual cascading, along with the associated fluid diffusion influenced by large-scale structures, and compare the results without the effect of Pedersen conductivity, such as those in the D-layer or Ref. [7].

2. Basic Equations

The nonlinear dynamics of solitary vortices driven by IGWs or low-frequency AGWs in stable stratified ionospheric incompressible fluids can be described by the following coupled equations [16].

$${\displaystyle \frac{\partial }{\partial t}}\left({\nabla }^2\psi -{\displaystyle \frac{\psi }{4H^2}}\right)+J(\psi ,{\nabla }^2\psi )=-{\displaystyle \frac{\partial \chi }{\partial x}}-{\displaystyle \frac{{\sigma }_p{B}_0^2}{{\rho }_0}}\left({\displaystyle \frac{{\partial }^2\psi }{\partial z^2}}-{\displaystyle \frac{\psi }{4H^2}}\right),$$

(1)

$${\displaystyle \frac{\partial \chi }{\partial t}}+J(\psi ,\chi )={\omega }_g^2{\displaystyle \frac{\partial \psi }{\partial x}},$$

(2)

where ${\nabla }^2\equiv \left({\partial }^2/\partial x^2\right)+\left({\partial }^2/\partial z^2\right)$, $\psi (x,z,t)(\sim e^{-z/H}\psi )$ is the velocity potential (stream function), $\chi (x,z,t)$ is the normalized fluid density perturbation, H is the fluid density scale height, ${\omega }_g=\sqrt{g/H}$ is the Brunt–Väisälä frequency with g denoting the constant gravitational acceleration, $J(a,b)=(\partial a/\partial x)(\partial b/\partial z)-(\partial a/\partial z)(\partial b/\partial x)$ denotes the Jacobian of two functions a and b, ${\rho }_0$ is the constant background fluid density, $B_0$ is the constant geomagnetic field, and ${\sigma }_p$ is the Pedersen conductivity. In the absence of the latter, Eqs. (1) and (2) reduce to those in Ref. [18].

Before, proceeding to study Eqs. (1) and (2) numerically for the turbulent wave dynamics of IGW-driven solitary vortices, we review some results as in Ref. [16], namely, the energy conservation, an estimation for the Pedersen conductivity, and the dispersion relation and the damping or growth rate of instability. From Eqs. (1) and (2), it is straightforward to obtain the following energy conservation equation [16].

$${\displaystyle \frac{\partial E}{\partial t}}=-{\displaystyle \frac{{\sigma }_p{B}_0^2}{{\rho }_0}}\int \int \left[{\left({\displaystyle \frac{\partial \psi }{\partial z}}\right)}^2+{\displaystyle \frac{{\psi }^2}{4H^2}}\right]dxdz,$$

(3)

where E is the total energy, given by,

$$E=\int \int \left[{\displaystyle \frac{1}{2}}{(\nabla \psi )}^2+{\displaystyle \frac{\psi }{8H^2}}+{\displaystyle \frac{{\chi }^2}{2{\omega }_g^2}}\right]dxdz.$$

(4)

Since the integral on the right-hand side of Eq. (3) is positive definite, the negative sign before it indicates that the wave energy dissipates with time by the influence of the Pedersen conductivity.

Typical forces influencing IGW dynamics include the pressure gradient force, gravity, the Coriolis force, and the $\mathbf{J}\times \mathbf{B}$ force. The current $\mathbf{J}$ can be associated with parallel, perpendicular (Pedersen), or Hall conductivity. Pressure gradient and gravity are prerequisites for AGWs. The Coriolis and $\mathbf{J}\times \mathbf{B}$ forces are not always relevant in all three ionospheric layers (D, E, F). For details, see Ref. [16]. In the D-layer, the $\mathbf{J}\times \mathbf{B}$ force (and, thus, Pedersen conductivity) is negligibly small compared to the Coriolis and other forces. In the E layer, this force matters mainly due to the Hall current. In contrast, in the F layer, the $\mathbf{J}\times \mathbf{B}$ force, driven by Pedersen conductivity, dominates, even over the Coriolis force. Therefore, in the model equations (1) and (2), the relevant forces are the pressure gradient, $\mathbf{J}\times \mathbf{B}$ (from the Pedersen current), and gravity.

In what follows, we estimate the ratio ${\sigma }_p{B}_0^2/{\rho }_0$, associated with the Pedersen conductivity ${\sigma }_p$. As in Refs. [16,19], the conditions ${\omega }_{\mathrm{ce}}{\omega }_{\mathrm{ci}}\gg {\nu }_e{\nu }_{\mathrm{in}}$ and ${\omega }_{\mathrm{ci}}\gg {\nu }_{\mathrm{in}}$ hold good in the F-layer, where ${\omega }_{\mathrm{ce}}=e{B}_0/m_e$ and ${\omega }_{\mathrm{ci}}=Ze{B}_0/m_i$ are, respectively, the electron and ion-cyclotron frequencies and ${\nu }_e={\nu }_{\mathrm{ei}}+{\nu }_{\mathrm{en}}$ with ${\nu }_{\mathrm{ei}}$, ${\nu }_{\mathrm{en}}$, and ${\nu }_{\mathrm{in}}$ denoting, respectively, the electron-ion, electron-neutral, and ion-neutral collision frequencies. Under these conditions, the Pedersen conductivity can be approximated as: ${\sigma }_p\approx ({\nu }_e{\nu }_{\mathrm{in}}/{\omega }_{\mathrm{ce}}{\omega }_{\mathrm{ci}}){\sigma }_{\Vert }$, where ${\sigma }_{\Vert }$ is the parallel conductivity. Using the typical F-layer values [16], $n/N\sim {10}^{-5}-{10}^{-3}$, ${\nu }_e\sim {\nu }_{\mathrm{ei}}\sim {10}^3{\mathrm{s}}^{-1}$, ${\nu }_{\mathrm{en}}\sim {10}^2{\mathrm{s}}^{-1}$, ${\nu }_{\mathrm{in}}\lesssim 10{\mathrm{s}}^{-1}$, ${\omega }_{\mathrm{ce}}\sim {10}^7{\mathrm{s}}^{-1}$, and ${\omega }_{\mathrm{ci}}\sim 3\times {10}^2{\mathrm{s}}^{-1}$, we estimate, ${\sigma }_p{B}_0^2/{\rho }_0\approx {\nu }_{\mathrm{in}}(n/N)\sim {10}^{-6}-{10}^{-4}{\mathrm{s}}^{-1}$. Next, we consider the linear analysis to obtain the dispersion relation and the damping rate associated with the Pedersen conductivity of Eqs. (8) and (9). Assuming the perturbations to vary as plane waves of the form, $exp(i{k}_{x}x+i{k}_{z}z-i\omega t)$, where $\omega$ is the wave frequency and $k_x$ and $k_z$ are, respectively, the wave vector components along the x- and z-directions, and letting $\omega ={\omega }_0+i{\Gamma }_0$, where $|{\Gamma }_0|\ll {\omega }_0$, we obtain the following dispersion relation [16]

$${\omega }_0^2={\displaystyle \frac{k_x^2{\omega }_g^2}{k^2+1/4H^2}},$$

(5)

and the damping rate, given by,

$${\Gamma }_0=-{\displaystyle \frac{1}{2}}{\displaystyle \frac{{\sigma }_p{B}_0^2}{{\rho }_0}}{\displaystyle \frac{k_z^2+1/4H^2}{k^2+1/4H^2}},$$

(6)

where $k^2=k_x^2+k_z^2$. From Eqs. (5) and (6), it is evident that the IGW mode has the characteristics similar to low-frequency ion-acoustic waves in plasmas, and the damping indeed occurs due to Pedersen conductivity and the damping rate is of the order of the ratio ${\sigma }_p{B}_0^2/{\rho }_0$, the key parameter in the model equations.

To begin the simulation study in Sec. Section 3, it is imperative to make Eqs. (1) and (2) dimensionless by redefining the variables as follows:

$$\begin{array}{cc}\hfill & (x,z)\to (x,z)/l_0,t\to t/t_0,\psi \to \psi /{\psi }_0,\hfill \\ \hfill & \chi \to \chi /{\chi }_0,{\chi }_0={\psi }_0/t_0{l}_0,H^{\prime }=H/l_0,\hfill \\ \hfill & \alpha ={\psi }_0{t}_0/l_0^2,{\omega }_g^{\prime }={\omega }_g{t}_0,\eta ={\sigma }_p{B}_0^2{t}_0/{\rho }_0,\hfill \end{array}$$

(7)

where $l_0$ and $t_0$ denote, respectively, the characteristic length and time scales. The dimensionless parameters $\alpha$, $\eta$, and ${\omega }_g^{\prime }$, respectively, correspond to the Jacobian nonlinearities, the normalized Pedersen conductivity, and the normalized frequency. Thus, Eqs. (1) and (2) reduce to

$${\displaystyle \frac{\partial }{\partial t}}\left({\nabla }^2\psi -{\displaystyle \frac{\psi }{4H^{\prime 2}}}\right)+\alpha J(\psi ,{\nabla }^2\psi )=-{\displaystyle \frac{\partial \chi }{\partial x}}-\eta \left({\displaystyle \frac{{\partial }^2\psi }{\partial z^2}}-{\displaystyle \frac{\psi }{4H^{\prime 2}}}\right),$$

(8)

$${\displaystyle \frac{\partial \chi }{\partial t}}+\alpha J(\psi ,\chi )={\omega }_g^{\prime 2}{\displaystyle \frac{\partial \psi }{\partial x}}.$$

(9)

3. Simulation Results and Discussion

We numerically study the nonlinear interactions of the large-scale velocity field potential (stream function) and small-scale density perturbations, as described by Eqs. (8) and (9). To this end, we develop a two-dimensional spectral code for the numerical integration of these equations [7,20]. The initial conditions consist of isotropic fluctuations with no imposed mean fields constructed using random amplitudes and phases in Fourier space. Spatial discretization is achieved in our simulations through a discrete Fourier representation of the turbulent fields. Numerical stability is maintained by applying proper dealiasing techniques to eliminate spurious Fourier modes. We use the fourth-order Runge-Kutta (RK4) scheme to solve the equations over time, with a time step, $dt=0.001$. The computational domain contains $512\times 512$ modes in two dimensions. Furthermore, ${\psi }_0\sim l_0^2/t_0$, ${\chi }_0\sim l_0/t_0^2$ so that $\alpha \sim 1$, and $\eta \sim ({10}^{-2}-{10}^{-1})$ depending on the magnetic field strength in the F-layer. Also, we take $H^{\prime }\sim 0.1$ and ${\omega }^{\prime }\sim 0.8$ throughout the simulation.

We note that the turbulent dynamics of IGWs, described by Eqs. (8) and (9) with $\eta =0$, are known in the literature [7] and may apply to the ionospheric D-layer. We reproduce the results for $\eta =0$ to compare with those in the F-layer. As observed earlier in Ref. [7] for $\eta =0$, the perturbed velocity potential forms large-scale structures. In contrast, density fluctuations cascade to smaller spatial scales. The coexistence of both large- and small-scale structures is the key feature of two-dimensional turbulent systems. Specifically, two-dimensional hydrodynamic turbulence in an incompressible fluid conserves both energy and mean-squared vorticity. Under external forcing, these conserved quantities exhibit a dual cascade. The forward cascade generates small-scale structures, while the inverse cascade leads to large-scale structures. Such dynamics typically appear in driven turbulence simulations, as in the present work. In these, specific modes are randomly excited in spectral space, redistributing energy and preserving physical invariants in k-space [7]. Similar features appear for $\eta =0$ (see left panels of Figure 1 and Figure 2). In the F-layer ($\eta =0.14$), the initially randomly propagating IGWs interact and generate large-scale velocity potential structures in the early phase (see right panels of Figure 1 at different times). As the system evolves, these structures break into intermediate-scale features due to the Pedersen conductivity. Meanwhile, density perturbations develop into smaller-scale structures (see right panels of Figure 2 at different times). Increasing values of the parameter $\eta$ also cause a faster decay of wave amplitudes, emphasizing the damping effect on the vortex dynamics.

During the decay of wave turbulence, energy in large-scale eddies cascades down to smaller scales, establishing a statistically steady inertial regime characterized by the forward transfer of one of the system’s invariants. Decaying turbulence frequently results in the emergence of coherent structure as it relaxes, which causes nonlinear interactions to become inefficient once they saturate. IGWs appear to exhibit turbulence similar to two-dimensional hydrodynamic turbulence, where density and velocity potential fluctuations correspond to fluid energy and vorticity, respectively [7]. Subplot (a) of Figure 3 illustrates the evolution of energy associated with the density and potential field modes (The case of $\eta =0$). Once the nonlinear interactions saturate, the turbulent energy stops growing and remains almost constant throughout the simulations. However, after nonlinear interactions saturate, the transfer of energy among turbulent modes in spectral space becomes inefficient, and the energy transfer rate becomes insignificant. Higher k-values dominate the energy, so dissipation by Pedersen conductivity is expected to dissipate the turbulent energy effectively. In contrast, subplot (b) of Figure 3 illustrates a significant decay of the turbulent energy. This behavior is primarily driven by dissipation associated with Pedersen conductivity ($\eta \ne 0$). As a result, higher Pedersen conductivity leads to a more rapid decay of turbulent energy.

The inertial range spectrum should contain characteristic features of strongly anisotropic turbulent flows. As such, the energy spectrum resulting from IGW turbulence under the effect of Pedersen conductivity may differ from the classical Kolmogorov form. To examine the spectral characteristics of IGW turbulence, the inertial range energy spectrum is presented in Figure 4. The inertial range energy spectrum exhibits a Kolmogorov-like $k^{-5/3}$ scaling when Pedersen conductivity is absent [$\eta =0$, Subplot (a)]. In the F-layer, when $\eta$ is increased to $0.14$, the spectrum becomes significantly steeper, with a spectral index close to $k^{-2.40}$ [Subplot (b)]. This index deviates from the usual Kolmogorov-like ($k^{-5/3}$) scaling, indicating that, in the F-layer, the effect of Pedersen conductivity leads to weaker turbulence compared to the D-layer or Ref. [7].

Finally, we estimate the turbulent transport coefficient arising from the self-consistent dynamics of small- and large-scale turbulent fluctuations. An effective electron diffusion coefficient due to momentum transfer can be expressed as

$$D_{\mathrm{eff}}={\int }_0^{\infty }\langle \mathbf{v}(\mathbf{r},t)·\mathbf{v}(\mathbf{r},t+t^{\prime })\rangle d{t}^{\prime },$$

(10)

where $\mathbf{v}$ is the fluid velocity, the angular bracket denotes spatial averaging, and ensemble averages are normalized to unit mass [7]. In our simulations, we compute the effective diffusion coefficient, $D_{\mathrm{eff}}$, from both the velocity field and its fluctuations to quantify turbulent transport associated with large-scale flows in the F-layer, for different values of Pedersen conductivity. As previously noted in Ref. [7], and as also demonstrated in our simulations [See Figure 5, subplot (a)], the effective cross-field transport remains low during the initial phase due to the Gaussian nature of field perturbations. As turbulence increases and large-scale structures begin to form, the diffusion coefficient grows rapidly, indicating stronger transport. Moreover, when the system reaches a steady state, the nonlinearly coupled modes form stationary structures, and $D_{\mathrm{eff}}$ eventually becomes saturated. The primary cause of the enhanced cross-field transport is the emergence of large-scale structures within turbulence dominated by IGWs. We also examine the fluid transport by considering the influence of the Pedersen conductivity. The effective diffusion coefficient, $D_{\mathrm{eff}}$, is derived from the cross-correlation of the particle drift velocity, as described by equation (10). As shown in Figure 5 (b), $D_{\mathrm{eff}}$ gradually levels off as the system transits into a steady-state regime with reduced energy transfer by the influence of $\eta$.

4. Conclusion

We have investigated the turbulent dynamics of IGW-driven vortex turbulence in the conductive ionosphere, focusing on the F-layer, where the Pedersen current plays a dominant role over the parallel and Hall currents, and compared the results with those without this effect, e.g., relevant in the D-layer. Our analysis shows that, in contrast to the D-layer or where the effect of the Pedersen current is insignificant, the large-scale perturbed velocity potential and small-scale density fluctuations in the F-layer exhibit a non-Kolmogorov-like energy spectrum. We also observe that the formation of intermediate-scale structures in the F-layer can lead to significantly lower cross-field diffusion than that in the D-layer. The system’s energy loss occurs due to magnetic viscosity, influenced by Pedersen conductivity in the F-region. This energy dissipation can lead to a reduction in the vortex structure [17].

To conclude, the turbulent dynamics of IGWs could help infer atmospheric and space weather conditions, understand the coupling of magnetospheric and ionospheric perturbations, and predict disruptions to radio signals.