Development and technical paper 10 Feb 2021
Development and technical paper  10 Feb 2021
Azimuthal averaging–reconstruction filtering techniques for finitedifference general circulation models in spherical geometry
 ^{1}CAS Key Laboratory of Geospace Environment, School of Earth and Space Sciences, University of Science and Technology of China, Hefei, China
 ^{2}Mengcheng National Geophysical Observatory, University of Science and Technology of China, Hefei, China
 ^{3}CAS Center for Excellence in Comparative Planetology, Hefei, China
 ^{4}Department of Earth Sciences, the University of Hong Kong, Pokfulam, Hong Kong SAR, China
 ^{5}High Altitude Observatory, National Center for Atmospheric Research, Boulder, CO, USA
 ^{6}Applied Physics Laboratory, Johns Hopkins University, Laurel, MD, USA
 ^{1}CAS Key Laboratory of Geospace Environment, School of Earth and Space Sciences, University of Science and Technology of China, Hefei, China
 ^{2}Mengcheng National Geophysical Observatory, University of Science and Technology of China, Hefei, China
 ^{3}CAS Center for Excellence in Comparative Planetology, Hefei, China
 ^{4}Department of Earth Sciences, the University of Hong Kong, Pokfulam, Hong Kong SAR, China
 ^{5}High Altitude Observatory, National Center for Atmospheric Research, Boulder, CO, USA
 ^{6}Applied Physics Laboratory, Johns Hopkins University, Laurel, MD, USA
Correspondence: Jiuhou Lei (leijh@ustc.edu.cn)
Hide author detailsCorrespondence: Jiuhou Lei (leijh@ustc.edu.cn)
When solving hydrodynamic equations in spherical or cylindrical geometry using explicit finitedifference schemes, a major difficulty is that the time step is greatly restricted by the clustering of azimuthal cells near the pole due to the Courant–Friedrichs–Lewy condition. This paper adapts the azimuthal averaging–reconstruction (ring average) technique to finitedifference schemes in order to mitigate the time step constraint in spherical and cylindrical coordinates. The finitedifference ring average technique averages physical quantities based on an effective grid and then reconstructs the solution back to the original grid in a piecewise, monotonic way. The algorithm is implemented in a community upperatmospheric model, the Thermosphere–Ionosphere Electrodynamics General Circulation Model (TIEGCM), with a horizontal resolution up to $\mathrm{0.625}{}^{\circ}\times \mathrm{0.625}{}^{\circ}$ in geographic longitude–latitude coordinates, which enables the capability of resolving critical mesoscale structures within the TIEGCM. Numerical experiments have shown that the ring average technique introduces minimal artifacts in the polar region of general circulation model (GCM) solutions, which is a significant improvement compared to commonly used lowpass filtering techniques such as the fast Fourier transform method. Since the finitedifference adaption of the ring average technique is a postsolver type of algorithm, which requires no changes to the original computational grid and numerical algorithms, it has also been implemented in much more complicated models with extended physical–chemical modules such as the Coupled Magnetosphere–Ionosphere–Thermosphere (CMIT) model and the Whole Atmosphere Community Climate Model with thermosphere and ionosphere eXtension (WACCMX). The implementation of ring average techniques in both models enables CMIT and WACCMX to perform global simulations with a much higher resolution than that used in the community versions. The new technique is not only a significant improvement in space weather modeling capability, but it can also be adapted to more general finitedifference solvers for hyperbolic equations in spherical and polar geometries.
Highlights.

The ring average technique is adapted to solve the issue of clustered grid cells in polar and spherical coordinates with a finitedifference method.

The ring average technique is applied to develop a $\mathrm{0.625}{}^{\circ}\times \mathrm{0.625}{}^{\circ}$ highresolution TIEGCM and more complicated geoscientific models with polar and spherical coordinates as well as finitedifference numerical schemes.

The highresolution TIEGCM shows good capability in resolving mesoscale structures in the ionosphere–thermosphere (I–T) system.
Mesoscale structures with a typical horizontal size of 100–500 km have gained more and more attention in research on the dynamics of the upperatmospheric system. A number of studies have been carried out to investigate these structures, including the formation and evolution of polar cap patches and tongues of ionization (Basu et al., 1995; Foster et al., 2005; Zhang et al., 2013), dynamics of ionospheric irregularities (Makela and Otsuka, 2012; Sun et al., 2015), variations of the polar thermospheric density anomaly (Crowley et al., 2010; Lühr et al., 2004), and the space weather effects of mesoscale electric field variability (Codrescu et al., 1995; Matsuo and Richmond, 2008; Zhu et al., 2018; Lotko and Zhang, 2018). These dynamic mesoscale structures have shown critical importance in understanding the physics of the solar–terrestrial system and in space weather predictions, which challenges the resolution and accuracy of numerical models of the upperatmospheric system in resolving these important mesoscale signatures.
Spherical or cylindrical coordinates are commonly used in solving geophysical problems, including the modeling of the upperatmospheric systems. As a workhorse for space weather research, a number of general circulation models (GCMs) for the coupled ionosphere–thermosphere (I–T) system have been developed based on spherical coordinates using finitedifference schemes (e.g., Richmond et al., 1992; FullerRowell et al., 1996; Ridley et al., 2006; Ren et al., 2009). However, restricted by the longitudinal grid resolution, current horizontal resolutions used in I–T GCMs are still insufficient for fully resolving mesoscale atmospheric structures, which are either marginal or subgrid. The latest released version of the community code, the Thermosphere–Ionosphere Electrodynamic General Circulation Model (TIEGCM), has a longitude–latitude resolution of $\mathrm{2.5}{}^{\circ}\times \mathrm{2.5}{}^{\circ}$. The Coupled Thermosphere–Ionosphere–Plasmasphere (CTIP) model has a latitude resolution of 2^{∘} and a longitude resolution of 18^{∘}, and the most recent version of the Global Ionosphere Thermosphere Model (GITM) has a flexible grid with a latitudinal resolution up to 0.3125^{∘}, but the typical longitudinal resolution remains 2.5^{∘} due to severe time step restrictions for globalscale calculations.
The major difficulty in increasing longitudinal resolution in sphericalgeometrybased GCMs is that the explicit time stepping is constrained by the clustering azimuthal cells near the pole due to the Courant–Friedrichs–Lewy (CFL) condition (Courant et al., 1928). A number of attempts have been proposed to address this coordinate singularity issue (e.g., Purser, 1988; Bouaoudia and Marcus, 1991; Williamson et al., 1992; Takacs et al., 1999; Fukagata and Kasagi, 2002; Prusa, 2018). To use a time step that is larger than the global minimum requirement from CFL conditions, one common method used in a spherical GCM is to employ a lowpass Fourier filter at polar latitudes, which removes nonphysical, highfrequency zonal waves generated due to numerical instability caused by the local violation of CFL conditions (e.g., Skamarock et al., 2008). Although the Fourier filter can maintain computational stability and permit a much larger temporal step, the applicability of the fast Fourier transform (FFT) filter method is problemdependent, which also creates barriers in moving models forward to finer spatial resolutions. Moreover, the linear filtering of zonal components generated through a nonphysical time step may decrease the accuracy of the model calculations near the polar region, which affects physical conservations of, e.g., mass, momentum, and energy that are essential for the longterm behavior of the GCM (Williamson and Browning, 1973).
Recently, Zhang et al. (2019) developed a new technique called the ring average method for hyperbolic equations to mitigate the CFL restrictions in spherical or polar geometry on the basis of the method originally proposed in the Lyon–Fedder–Mobarry (LFM) magnetohydrodynamics (MHD) simulations (Lyon et al., 2004). The method is a “postsolver” type of algorithm applied after solving all the physical quantities in the original spherical coordinates; thus, no modification to the numerical solver or the computational grid is required when applying the ring average. Test simulation results have shown the effectiveness of the ring average algorithm in increasing the time step by a factor of 100 while maintaining the fidelity of the solutions. The original ring average technique was developed for solving hyperbolic equations in spherical or polar geometry based on finitevolume schemes, which redistributes the solution azimuthally through a conservative averaging–reconstruction algorithm. The finitevolume version of the ring average technique not only releases the time step constraint in spherical geometry, but also keeps the conservative nature of finitevolume schemes to machine precision. In this paper, we adapt the ring average technique to finitedifference schemes for solving hyperbolic equations. Defined on an effective reduced polar grid, the finitedifference adaption of the ring average technique also conducts a postsolver step of averaging–reconstruction in each azimuthal ring to maintain the numerical stability and relax the severe computational time step constraint. To demonstrate the effectiveness of the finitedifference version of the ring average technique, we use solutions from both linear advection equations and the TIEGCM as test beds. The ring average algorithm enables the use of a highresolution TIEGCM, such as $\mathrm{0.625}{}^{\circ}\times \mathrm{0.625}{}^{\circ}$ in longitude and latitude, with reasonable time steps and minimal numerical artifacts. Further applications of the technique to the Coupled Magnetosphere–Ionosphere–Thermosphere (CMIT) model and the Whole Atmosphere Community Climate Model with thermosphere and ionosphere extension (WACCMX) are also addressed.
This paper is organized as follows: in Sect. 2, we describe the details of the model and the ring average technique to solve the problem of clustering of polar grid cells. A hydrodynamic convection experiment with the ring average technique has also been conducted to test the capability of the method. Section 3 shows the preliminary results of the highresolution TIEGCM with the ring average technique implemented and the further applications of the technique. The findings of this work are summarized in Sect. 4.
2.1 Ring average in the finitedifference form
An example of standard polar grids with a horizontal resolution of $\mathrm{2.5}{}^{\circ}\times \mathrm{2.5}{}^{\circ}$ (longitude × latitude) in the TIEGCM is shown in Fig. 1. It is evident that in Fig. 1a, the azimuthal (longitudinal) computational nodes in the standard polar grid are significantly clustered near the pole, even with 144 cells in the azimuthal direction, resulting in very “thin” cells with small azimuthal extensions, which restricts the explicit time step for the advection equations. This azimuthal clustering becomes even worse when grid resolution increases; the time step drops to $\mathrm{1}/\mathrm{4}$ while the grid resolution doubles, corresponding to an increase in computational resources of at least 32 times, which becomes expensive, especially for global simulations with high spatial resolutions, in order to resolve mesoscale structures.
The finitedifference adaption of the ring average algorithm is based on a similar averaging–reconstruction process over a predefined “effective” azimuthal grid as used in the finitevolume version of the algorithm. Figure 1b shows an example of an effective polar grid for applying the finitedifference ring average technique. In the polar grid shown in Fig. 1b, since the reconstructed solution is monotonic within each effective computational cell, a much larger time step is allowed compared to the original grid shown in Fig. 1a. As shown in Fig. 1b, the effective longitudinal grid resolutions have been reduced and are less clustered towards the pole. For the most inside (highest latitude) grids, the 144 azimuthal cells (Fig. 1a) have been grouped into 9 effective cells (chunks), with 16 original cells in each chunk. Moving away from the pole, more chunks are employed. As an example, the numbers of chunks from inside to outside shown in Fig. 1b in the effective grid are 9, 9, 18, 18, 36, 36, 72, 72, 72, and 72, allowing a relatively smooth transition in the size of the cells going radially outward. Note that the choice of the number of chunks in each ring is nonunique. Numerical tests with finitevolume solvers have shown that the computational solution, under both smooth and discontinuous flow conditions, is insensitive to small changes in the chunk configuration (Zhang et al., 2019).
We use the following example of solving the linear advection equation to illustrate the averaging–reconstruction process within each chunk. Consider the following linear advection equation of an incompressible fluid in the azimuthal direction as an example:
where v is the advection velocity, ρ is the density profile, and x is the azimuthal dimension (x∈[0 2π]) along one ring. Assuming the x direction is uniformly discretized into N_{total} computational cells with $\mathrm{\Delta}x=\frac{\mathrm{2}\mathit{\pi}}{{N}_{\mathrm{total}}}$, a centraldifference forward Euler form of Eq. (1) for density ρ in cell k between time n and n+1 is written as
where k denotes the index of an individual thin cell in the original azimuthal grid. Δt is the time step regulated by the CFL condition. Without a ring average type of treatment, the time step Δt is restricted by the fact that thin azimuthal cells cluster near the pole. The ring average technique takes the average solution within a chunk m that contains 2s cells in the original grid as shown in Fig. 2. Summing over the finitedifference form of Eq. (1) within chunk m gives
Then, summing over the k indices within chunk m, Eq. (3) becomes
where ${\mathit{\rho}}_{is+\frac{\mathrm{1}}{\mathrm{2}}}^{n}$ and ${\mathit{\rho}}_{i+s+\frac{\mathrm{1}}{\mathrm{2}}}^{n}$ are the left and right values on the boundary of chunk m, as indicated by the red triangles in Fig. 2. The lefthand side of Eq. (7) is basically the time rate of the change in terms of the chunk density ϱ_{m}:
where ${\mathit{\varrho}}_{m}^{n+\mathrm{1}}=\sum _{k=is+\mathrm{1}}^{i+s}{\mathit{\rho}}_{k}^{n+\mathrm{1}}$ and ${\mathit{\varrho}}_{m}^{n}=\sum _{k=is+\mathrm{1}}^{i+s}{\mathit{\rho}}_{k}^{n}$. If assuming smoothness of the solution, which applies to typical upperatmospheric flow conditions, and using a piecewise linear reconstruction for the two interface values at time level n,
the righthand side of Eq. (7) is in the form of a centraldifference approximation of the spatial derivative $\frac{\partial \mathit{\varrho}}{\partial x}$ in chunk m:
Equating Eqs. (8) and (11) and considering the fact that the ΔX in computing the chunk derivative is actually 2sΔx, we obtain
where ΔT=2sΔt. Equation (12) is in the same numerical differential form of the advection equation in terms of the chunk density ϱ in the effective grid within the same order of finitedifference approximation:
Equation (12) also suggests that in principle the ring average method is capable of using a time step that is approximately 2s times larger than the original Δt restricted by the thin cells (assuming the CFL condition is dominated by the azimuthal direction in the innermost ring). Note that the above derivation of the finitedifference version of the ring average algorithm is independent of the numerical schemes solving the linear advection equation (Eq. 1). Thus, the ring average algorithm requires no modifications to the existing hydrodynamic equations solved by GCMs. Furthermore, since the ring average algorithm is applied after all the variables are solved on the original spherical grid, it requires no changes to the existing computational grid.
In the reconstruction step, the above algorithm uses the piecewise linear method (PLM) to reconstruct solutions within each chunk for the next time step of the GCM calculations, resulting in secondorder accuracy. To achieve higher accuracy in the reconstruction step, a piecewise parabolic reconstruction method (PPM) (Colella and Woodward, 1984) may be used in the algorithm, which provides fourthorder accuracy for the reconstruction step. In the following section when applying the ring average algorithm in a GCM, we use both PPM and PLM for different variables. The criteria for using PPM or PLM here depend on their spatial gradient from the fluid calculations. For variables that have a relatively greater spatial gradient, we use the PPM to reach high accuracy and maintain stability; otherwise, the PLM is used for the calculations.
The algorithm shown in Fig. 3 illustrates the steps of applying the ring average technique using either PLM or PPM. The steps consist of chunk division, chunk averaging, and reconstruction. The averaging–reconstruction process (ring average) in this study is similar to Zhang et al. (2019), with modifications to the reconstruction method (PPM or PLM) adapted to finitedifference schemes. The detailed procedures of the ring average technique in this study are described as follows.
2.1.1 For variables using the PLM reconstruction

Step 1. Divide the azimuthal grid cells into chunks and pull data into the chunks.

Step 2. Calculate the average value F_{A}, left interface value F_{L}, and right interface value F_{R} at chunk m (m is the index of the chunk number in an azimuthal ring). F_{L} and F_{R} are the interface values in each chunk determined by the following parabola functions:
$$\begin{array}{}\text{(14)}& {\displaystyle}& {\displaystyle}{F}_{\mathrm{L}}=\left({F}_{m\mathrm{2}}+\mathrm{7}{F}_{m\mathrm{1}}+\mathrm{7}{F}_{m}{F}_{m+\mathrm{1}}\right)/\mathrm{12},\text{(15)}& {\displaystyle}& {\displaystyle}{F}_{\mathrm{R}}=\left({F}_{m\mathrm{1}}+\mathrm{7}{F}_{m}+\mathrm{7}{F}_{m+\mathrm{1}}{F}_{m+\mathrm{2}}\right)/\mathrm{12},\end{array}$$where F_{m−2}, F_{m−1}, F_{m}, F_{m+1}, and F_{m+2} are the average values F_{A} at chunks with an index of m−2, m−1, m, m+1, and m+2, respectively.

Step 3. Reconstruct the variables by interpolating the average data linearly in each chunk:
$$\begin{array}{}\text{(16)}& {F}_{k}=\left(\mathrm{1}{\displaystyle \frac{k}{N}}\right){F}_{\mathrm{L}}+{\displaystyle \frac{k}{N}}{F}_{\mathrm{R}},\end{array}$$where N is the number of cells within each chunk and k is the local index ranging from 1 to N.

Step 4. Redo the above procedures for the next azimuthal ring until the ring average is not needed.
2.1.2 For variables using the PPM reconstruction
The procedures in PPM are the same with PLM except for Step 3.
Step 3. Reconstruct the variables parabolically in each chunk using the following function:
where A, B, and C are constants representing the parabolic function, which connects F_{L} and F_{R}:
2.1.3 For vector variables using the PLM reconstruction and Fourier reduction
Step 0 in Fig. 3 corresponds to a Fourier expansion (reduction) step that is required for vector GCM variables in spherical coordinates before applying the ring average process. The main purpose of the Fourier reduction step is to maintain the direction of vectors after ring averaging, especially for the neutral meridional and zonal wind across the pole. Thus, only the second and higher Fourier components of the data in the azimuthal cell are smoothed using the ring average filter, while the zero and first Fourier components are kept unchanged. Here are the details of the Fourier expansion process.
Step 0. Calculate the Fourier components of the azimuthal data:
where i is the thin cell index along the azimuthal direction ranging from 1 to N_{total}, N_{total} is the total thin cell number in the azimuthal direction, F_{i} represents the second and higher Fourier components that will be later reconstructed, and A_{0}, A_{1}, and B_{1} are the zero and first Fourier coefficients:
The higher Fourier components F_{i} are pulled into chunks for the ring average processes.
Step 1–4. The steps are the same as the above PLMs, except that the reconstructed data, ${F}_{i}^{\prime}$, are brought back together with the first two Fourier components after the reconstruction:
2.2 Ring average for the advection equation
In this section, in order to illustrate the implementation of the ring average algorithm in a finitedifference code, we solve the twodimensional (2D) linear advection equation in the polar geometry as an example. The code used in the 2D linear advection solver is a main subroutine used in the ring average module for GCMs. This twodimensional advection test in polar geometry is also useful to demonstrate the effectiveness of the finitedifference ring average technique in handling a strong, narrow shear flow near the pole. A fourthorder centralfinitedifference scheme is used to solve the following mass continuity equation under the incompressible assumption:
where ρ is the density, and u is the timestationary flow velocity defined in polar coordinates (r, θ). The polar geometry of this test is defined with a resolution of 0.625^{∘} in both longitude and latitude, with 144 cells in the r direction uniformly distributed within (0, 1) and 576 cells in the θ direction uniformly distributed within (0, 2π).
Figure 4a shows the initial state (t=0) of ρ, ranging linearly along the y direction from a magnitude of 2 at the topside to 0.01 at the bottom side:
The timeindependent flow velocity u is set as a Gaussiandistributed shear flow towards −y and centered at x=0.15 with a peak velocity of −1 and a halfwidth of 0.01:
As simulation time evolves, a large density gradient occurs near the pole driven by the timestationary shear flow, with its pattern following the analytical distribution of the flow velocity u. Figure 4a–c show three snapshots of the density in the linear advection experiment at t=0, t=0.75, and t=1.5 using the finitedifference version of the ring average technique with PPM reconstructions, as described in Sect. 2.1. For comparison, Fig. 4d–f show the corresponding snapshots derived from another simulation using an FFT filter, and Fig. 4g–i show the results at the same simulation time calculated from the fourthorder finitedifference scheme without applying any filtering technique. In the simulation using the ring average, the number of averaging chunks in each azimuthal ring near the pole is set to be [18 18 18 18 36 36 36 36 72 72 72 72 144 144 144 144] from the first ring to the 16th ring, as indicated by the white circles in the top panels around 80^{∘} N. For the FFT filter, a Fourier expansion is applied in the azimuthal direction at each time step to the fluid density. Waves with frequencies that are higher than the cutoff frequencies are eliminated from the Fourier spectra of prognostic variables. The values of prognostic variables are then reconstructed through an inverse Fourier transform using the modified Fourier spectra. Each latitude grid has its own cutoff frequency, and the wave number to be cut off in this experiment near the pole is set to be [1 1 2 2 2 2 4 4 4 4 8 8 8 8 10 10], which is similar to the TIEGCM FFT filter spectrum (Wang, 1998).
As shown in Fig. 4b and c, density structures with a large spatial gradient flow across the pole as time progresses. Compared with the nonfilter case in Fig. 4h and i, no evident numerical instability or artificial structure occurred when applying the ring average technique. In contrast, the density structure using an FFT filter in Fig. 4d–f exhibits numerical oscillations in the radial direction, together with an artificial depletion of density near the pole. This density depletion is due to the nonconservative nature of the FFT method by truncating highspatialfrequency wave modes in a linear way. Figure 4j–l show a onedimensional comparison of the density profiles along x=0.15, with the region of averaging chunks denoted by yellow. The comparisons suggest that the density flow is not noticeably affected by the implementation of the ring average technique in the finitedifference solver. Note that the time step used after applying the ring average technique is 0.0001 s, which is 25 times larger than that used in the simulation without the ring average (dt=0.000004 s). Although the FFT filter can result in nonoscillatory solutions in the finitedifference solver, as shown by the onedimensional profiles in Fig. 4l, evident density oscillations occur near the pole due to numerical instability caused by the FFT method. The cutoff frequency of the FFT filter is casedependent and has a problem of mass loss compared to the ring average method. The advection experiment illustrates that the ring average technique is capable of relaxing the severe time step constraint and resolving large density gradients when passing through the clustered grid cells near the pole.
2.3 Ring average for GCMs
We use the NCAR TIEGCM to demonstrate the effectiveness of the ring average technique in resolving mesoscale upperatmospheric structures. TIEGCM is a physicsbased 3D global model that solves the coupled equations of momentum, energy, and continuity for neutral and ion species in the upperatmospheric I–T system using a fourthorder and centered finitedifference scheme to evolve the advection terms on each pressure surface with a staggered vertical grid (Qian et al., 2014; Roble et al., 1988; Richmond et al., 1992). The TIEGCM utilizes a spherical coordinate system fixed with respect to the rotating Earth, with geographic latitude and longitude as the horizontal coordinates and the pressure surface as the vertical coordinate. The following is a brief introduction of the basic equations in the TIEGCM.
The thermospheric energy equation is
with temperature T_{n}, time t, the vertical coordinate $Z=\mathrm{ln}({P}_{\mathrm{0}}/P)$, the pressure P, and the reference pressure P_{0}; g is gravity, K_{T} is the molecular thermal conductivity, C_{P} is the specific heat per unit mass, H is the pressure scale height, K_{E} is the eddy diffusion coefficient, ρ is the atmospheric mass density, U is the horizontal neutral velocity with the zonal and meridional components u_{n} and v_{n}, w is the vertical velocity defined by $w=\mathrm{d}Z/\mathrm{d}t$, R^{*} is the universal gas constant, $\stackrel{\mathrm{\u203e}}{m}$ is the mean atmospheric mass, and Q and L are the heating and cooling rates. The mean molecular mass $\stackrel{\mathrm{\u203e}}{m}$ is determined by
where Ψ and m represent the mass mixing ratio and the molecular mass, respectively, for the three thermospheric major species O_{2}, O, and N_{2}.
The zonal momentum equation is expressed as
and the meridional momentum equation is
where λ and φ represent the geographic latitude and longitude, respectively. R_{E} is the radius of the Earth, μ is the viscosity coefficient, which is the sum of eddy and molecular viscosity coefficients, f is the Coriolis parameter, ϕ is the geopotential, H is the pressure scale height, v_{i} and u_{i} are the meridional and zonal E×B ion drift velocities, and λ_{xx}, λ_{xy}, λ_{yx}, and λ_{yy} are the iondrag tensor coefficients. The TIEGCM “vertical velocity” $w=\mathrm{d}Z/\mathrm{d}t$ is determined by solving the continuity equation:
The real vertical velocity is obtained by first integrating the continuity equation (Eq. 29) over Z to get w and then multiplying w by the neutral pressure scale height to get the right unit.
The thermospheric major species in the TIEGCM include O_{2}, O, and N_{2}. The continuity equation for the mass mixing ratio of O_{2} and O is given by
where $\stackrel{\mathrm{\u0303}}{\mathrm{\Psi}}=({\mathrm{\Psi}}_{{\mathrm{O}}_{\mathrm{2}}}$, Ψ_{O}), τ is the diffusion timescale equal to 1.86×10^{3} s, ${m}_{{\mathrm{N}}_{\mathrm{2}}}$ is the molecular mass for molecular nitrogen, T_{00}=273 K is the standard temperature, $\stackrel{\mathrm{\u0303}}{\mathit{\alpha}}$ is the matrix operator of the diffusion coefficients, K(Z) is the eddy diffusion coefficient, and $\stackrel{\mathrm{\u0303}}{S}$ and $\stackrel{\mathrm{\u0303}}{R}$ are the production and loss terms for these two species. The diagonal matrix operator L has elements of the form
where i=1 and i=2 denote O_{2} and O, respectively. The N_{2} mass mixing ratio is determined by
The minor species in the TIEGCM are N(^{4}S), N(^{2}D), and NO. The timescale of N(^{4}S) is relatively short and is thus considered to be in photochemical equilibrium. N(^{4}S) and NO have longer lifetimes, so the transport effects must be taken into account. The governing equation for these two species is
where
where $\stackrel{\mathrm{\u0303}}{\mathrm{\Psi}}=({\mathrm{\Psi}}_{\mathrm{NO}}$, ${\mathrm{\Psi}}_{\mathrm{N}{(}^{\mathrm{4}}\mathrm{S})})$, $\stackrel{\mathrm{\u0303}}{A}$ is the vertical molecular diffusion coefficient, and $\stackrel{\mathrm{\u0303}}{S}$ and $\stackrel{\mathrm{\u0303}}{R}$ are the production and loss terms for each species. Terms in $\stackrel{\mathrm{\u0303}}{E}$ represent the effects of gravity, thermal diffusion, and the frictional interaction with the major species on the vertical profiles of these two species. $\stackrel{\mathrm{\u0303}}{F}$ is a matrix operator for the frictional interactions, $\stackrel{\mathrm{\u0303}}{\mathit{\alpha}}$ is the thermal diffusion coefficient, and $\stackrel{\mathrm{\u0303}}{m}$ is the molecular mass for the two minor species.
The ions of the ionosphere in the TIEGCM include O^{+}, ${\mathrm{O}}_{\mathrm{2}}^{+}$, NO^{+}, N^{+}, and ${\mathrm{N}}_{\mathrm{2}}^{+}$, and the electron density is calculated by chemical equilibrium of these ions. All major ionospheric ions except O^{+} are assumed to be in photochemical equilibrium; thus, their densities can be calculated simply by balancing the loss and production rates. The O^{+} density is determined not only by O^{+} loss and production but also by transport due to E×B drifts, neutral winds, and fieldaligned ambipolar diffusion. The O^{+} continuity equation can be expressed as
where n is the O^{+} density, Q is the production rate, L is the loss rate, and ▿⋅(nV) is the transport term. The ion velocity V is given by
where V_{‖} (caused by ambipolar diffusion and neutral winds) and V_{⟂} (caused by E×B drifts) are the parallel and perpendicular velocities with respect to the magnetic field line, respectively, b is a unit vector along the magnetic field, ν_{in} is the ionneutral collision frequency, g is the acceleration due to gravity, ρ_{i} is the ion mass density, P_{i} and P_{e} are the ion and electron pressures, respectively, U_{n} is the neutral velocity, B is the magnetic field, and E is the electric field.
By assuming a thermal quasisteady state, the electron energy equation is
with I the geomagnetic dip angle, K^{e} the electron thermal conductivity parallel to the magnetic field, ∑Q_{e} the sum of all local electron heating rates, and ∑L_{e} the sum of all local cooling rates.
For the electrodynamics, i.e., the “neutral wind dynamo process”, TIEGCM assumes steadystate electrodynamics with a divergencefree current density J for longer timescales:
where σ_{P} and σ_{H} are the Pedersen and Hall conductivities, respectively. ${\mathbf{J}}_{\left\right}$ and J_{M} are the ohmic component of current density parallel to the magnetic field and the nonohmic magnetospheric component, respectively.
The ionospheric convection pattern for computing the plasma advection velocity V_{⟂} at high latitudes is specified by either the Heelis et al. (1982) or the Weimer (2005) empirical model, while at the bottom boundary the migrating tides are specified using the Global Scale Wave Model (Hagan and Forbes, 2002, 2003). The current standard version of TIEGCM (TIEGCM v2.0) provides two spatial resolution options: (1) $\mathrm{5}{}^{\circ}\times \mathrm{5}{}^{\circ}$ in horizontal geographic latitude–longitude grid and $\mathrm{1}/\mathrm{2}$ scale height in the vertical direction or (2) $\mathrm{2.5}{}^{\circ}\times \mathrm{2.5}{}^{\circ}$ in horizontal geographic latitude–longitude grid and $\mathrm{1}/\mathrm{4}$ scale height in the vertical direction.
In this study, the ring average technique is implemented in the TIEGCM v2.0 to solve the issue of clustering grid cells near the poles in the development of a highresolution version of the TIEGCM. This technique is applied as a postprocessing treatment of the fluid variables including oxygen ion density O^{+}, neutral temperature T_{n}, thermospheric compositions Ψ, and meridional, zonal, and vertical winds (U_{n}, V_{n}, w) at each time step, with different reconstruction methods (PPM or PLM) for different variables (Table 1). Due to the use of MPI parallelization in the TIEGCM in supercomputers, the ring average technique firstly collects the azimuthal data in the root thread, conducts the averaging–reconstruction process, and finally redistributes data into each MPI thread. Figure 5 illustrates ring average filters used in the main algorithms of the TIEGCM, including the thermosphere solvers in Equations (25)–(34), the ionosphere solver for O^{+} in Eqs. (35)–(39), and the dynamo solver for electrodynamic coupling in the Eq. (40). For neutral variables in the thermosphere solver, the ring average technique with the PLM reconstruction is utilized. Specifically, for the meridional and zonal neutral winds, the second and higher Fourier components are processed with the PLM ring average filter to maintain the direction of vectors across the pole, as displayed in Fig. 5. The oxygen ion (O^{+}) in the ionosphere usually has much sharper gradients than the neutral variables, e.g., tongue of ionization (TOI) structures; thus, the PPM is used in the reconstruction process to provide highorder accuracy and handle the larger local gradient. Meanwhile, to balance the numerical stability and computational speed, a subcycling technique, which has a smaller time step for O^{+} than neutral variables, has been applied in the O^{+} solver because the ions can move much faster than the neutrals with the E×B drifts, especially during major geomagnetic storms.
On the basis of the ring average technique, a new highresolution version of TIEGCM with a horizontal longitude–latitude resolution as high as $\mathrm{0.625}{}^{\circ}\times \mathrm{0.625}{}^{\circ}$ is developed. Table 2 lists the ring average setup used in different TIEGCM resolutions. The third column in Table 2 represents the number of averaging chunks in each azimuthal ring near the pole from the first innermost to the outermost rings. For example, in the first azimuthal ring near the pole of the $\mathrm{0.625}{}^{\circ}\times \mathrm{0.625}{}^{\circ}$ grid resolution, 64 longitude grids ($\mathrm{576}/\mathrm{9}=\mathrm{64}$) and 40 longitude degrees ($\mathrm{360}{}^{\circ}/\mathrm{9}=\mathrm{40}{}^{\circ}$) are grouped into a chunk. In the outermost filtered ring (around 71.25^{∘} latitude), one averaging chunk only contains two longitude grids. Table 3 summarizes the information on different spatial resolutions of the TIEGCM, including the current version of the $\mathrm{2.5}{}^{\circ}\times \mathrm{2.5}{}^{\circ}$ TIEGCM with the default FFT filter and the $\mathrm{2.5}{}^{\circ}\times \mathrm{2.5}{}^{\circ}$, $\mathrm{1.25}{}^{\circ}\times \mathrm{1.25}{}^{\circ}$, and $\mathrm{0.625}{}^{\circ}\times \mathrm{0.625}{}^{\circ}$ resolution TIEGCM with the ring average filter. As the resolution doubles, the time step decreases approximately linearly rather than quadratically. In practice, the $\mathrm{0.625}{}^{\circ}\times \mathrm{0.625}{}^{\circ}$ resolution of the code runs about 2 times faster than real time with 256 processors on the NCAR/CISL Cheyenne supercomputer system (12 h for a 1 d geomagnetic storm simulation), which is at fairly low computational cost for mesoscaleresolving global simulations. The preliminary results of the highresolution TIEGCM will be shown in the following section.
To show the capability of the new highresolution TIEGCM based on the ring average technique in resolving mesoscale I–T structures, we have simulated the ionospheric and thermospheric variations during the 17 March 2013 major geomagnetic storm as an example. Figure 6 displays the comparison of polar maps of electron densities between different filter techniques with the $\mathrm{2.5}{}^{\circ}\times \mathrm{2.5}{}^{\circ}$ horizontal resolution. The electron density is plotted on pressure surface 2, which is near the F_{2}region peak (∼300 km of altitude). Figure 6a corresponds to the standard TIEGCM with the FFT filter, while Fig. 6b is the result using the ring average technique. Generally, the electron densities in the two simulations in Fig. 6a and b are similar below 60^{∘} N, with an evident electron density enhancement seen in the afternoon sector and negative storm effects in the morning at 10:50 UT during the storm. The dense ionospheric plasma in the afternoon sector is transported in the antisunward direction into the polar cap region by the dusk cell of the convection pattern. Consequently, prominent polar tongue of ionization (TOI) features can be seen as a narrow density plume on the dayside, which stretches from 65^{∘} N at noon to latitudes greater than 80^{∘} N inside the polar cap. Those TOI features agree well with the polar Global Position System (GPS) total electron content (TEC) observations (e.g., Foster et al., 2005; Thomas et al., 2013). It is evident that, in Fig. 6a, the TOI cannot go through the polar cap region and generates an artificial “hole” structure above 80^{∘} N. This nonphysical depletion is associated with the loss of electron density induced by the removal of high frequency in the FFT filter, as also indicated in the advection experiment in Fig. 4f. Consequently, the plasma within the TOI accumulates around the hole and a “ringlike” structure appears at about 70^{∘} N. In contrast, for the ring average technique, the electron densities in Fig. 6b can successfully flow through the polar cap and arrive at the nightside, which is consistent with Fig. 4c. Thus, using the ring average, the artificial structures no longer exist in the polar cap region, indicating the advantage of the ring average technique in handling numerical instability by causing fewer artificial structures and preserving the real mesoscale structures.
Figure 7 shows the comparison of the polar maps of electron densities between simulations with different spatial resolutions using the TIEGCM bolstered by the ring average technique. The simulation results after using the ring average technique are generally similar among different simulations, with finer structures at higher spatial resolutions. Besides the ionospheric parameters, we have also tested the performance of the ring average in thermospheric simulations (not shown here), which indicates that the thermospheric variables generally converge between different spatial resolutions. The thermospheric temperature, $\mathrm{O}/{\mathrm{N}}_{\mathrm{2}}$, and thermospheric density simulated by two kinds of filters do not show distinct differences compared with the ionospheric simulations due to the relatively smoother variations of neutral parameters. Only slight deviation exists locally on a smaller scale in the polar thermosphere. The results from Figs. 6 and 7 demonstrate that the ring average technique can be applied in the finitedifference method, which is usually considered to be less stable than the finitevolume scheme. The ring average method can successfully maintain numerical stability, even with the structures of large spatial gradients, and conserve true mesoscale structures. Meanwhile, the ring average technique shows advantages of inexpensive computational cost and easy implementation, as indicated by Table 3. By using the ring average, the time step has been greatly relaxed in the ideal advection experiment and the highresolution TIEGCM, which would maintain the computational cost at an acceptable level. Furthermore, ring averaging can be applied as a postfilter after each simulation step and would not require a modification of the underlying code; this would make the technique easily applied.
Benefiting from the ring average technique, the newly developed highresolution TIEGCM has been applied to explore the mesoscale variations in the I–T system during space weather events. For instance, based on the $\mathrm{0.625}{}^{\circ}\times \mathrm{0.625}{}^{\circ}$ highresolution TIEGCM simulations and satellite observations, Dang et al. (2019) reported the occurrence of double TOIs and carried out a comprehensive study on the dynamic evolution and formation mechanism of double TOIs. Lu et al. (2020) used the highresolution model to study ionospheric disturbances such as traveling ionospheric disturbances and enhanced storm density during geomagnetic disturbances. The highresolution TIEGCM has also been utilized to simulate the subauroral polarization stream (Lin et al., 2019), neutral wind variabilities (Wu et al., 2019), and the responses of the I–T system to solar eclipses (Dang et al., 2018a, b; Lei et al., 2018; Wang et al., 2019). These works highlight the enhanced capability of the highresolution TIEGCM in resolving ionospheric and thermospheric mesoscale structures enabled by the ring average technique.
Simulating the mesoscale structures also requires a more realistic input from the upper boundary, corresponding to the electric field and auroral precipitation from the magnetosphere, and the bottom boundary, corresponding to the upward propagation of tides and waves from the lower atmosphere, of the I–T system. In the TIEGCM, these inputs are directly adopted from two empirical models, the Weimer model and the Global Scale Wave Model (GSWM), which might not necessarily represent the complexity of the actual physical processes from the boundaries. To obtain a more physical upper boundary condition, the CMIT has been developed (Wang et al., 2004; Wiltberger et al., 2004), which couples the LFM global magnetosphere model with the I–T model TIEGCM. The LFM provides the TIEGCM with highlatitude electric fields and auroral electron precipitation, and the TIEGCM feeds ionospheric heightintegrated conductance back to the LFM. The standard resolution of the ionosphere and thermosphere in CMIT is $\mathrm{2.5}{}^{\circ}\times \mathrm{2.5}{}^{\circ}$, which is the same as the standard TIEGCM. By implementing the highresolution TIEGCM in CMIT, the thermosphere–ionosphere part in CMIT has a horizontal resolution of 1.25^{∘} in both latitude and longitude, which is comparable to the magnetospheric resolution of 100 km mapped to the ionospheric reference altitude. Figure 8 shows an example CMIT simulation of the ionosphere and thermosphere at 17:30 UT during the 17 March 2013 geomagnetic storm. The TEC in Fig. 8a shows more dynamic and finer TOI variations driven by the magnetospheric convection during the storm time. Meanwhile, the thermospheric temperature in Fig. 8b also exhibits distinct mesoscale structures associated with changes in the neutral wind circulation and ion collisional heating. The results illustrate that, with the implementation of the ring average technique, the highresolution CMIT show advantages in resolving the dynamic evolution of mesoscale structures in the coupled magnetosphere–ionosphere–thermosphere system.
Furthermore, the ring average technique has also been applied in the WACCMX, which can provide a relatively more realistic bottom boundary for the I–T simulation. The WACCMX is a wholeatmosphere chemistry–climate general circulation model spanning the range of altitude from the Earth's surface to the upper thermosphere to simulate the entire atmosphere and ionosphere (Liu et al., 2018). The ionosphere and electrodynamo parts in WACCMX are the same as in the TIEGCM. The ring average scheme has been successfully implemented in the O^{+} transport module of the WACCMX to get a higher spatial resolution for the ionosphere. Figure 9 shows the simulation results for the 17 March 2013 geomagnetic storm from WACCMX. For this simulation, the horizontal resolution is $\mathrm{1.25}{}^{\circ}\times \mathrm{0.9}{}^{\circ}$ in longitude and latitude directions, respectively, and the vertical resolution in the upper atmosphere is $\mathrm{1}/\mathrm{4}$ of the scale height. Detailed analyses and exploration of the CMIT and WACCMX results are beyond the scope of this study and will be studied in the future. Ongoing efforts also include improving the resolution of the vertical direction and the electrodynamo of the TIEGCM as well as applying the ring average technique to highresolution data assimilation and space weather prediction.
In summary, a postprocessing technique with an averaging–reconstruction (ring average) algorithm is developed to solve the problem of clustering of azimuthal cells in a spherical coordinate with the finitedifference method. The ring average technique is applied based on a reduced effective polar grid by first averaging quantities within azimuthal effective “chunks” and then reconstructing them within each chunk. The ring average technique shows the advantages of inexpensive computational cost, easy implementation, time step relaxation, and maintenance of the mesoscale structures without introducing artifacts, which allows for the development of highresolution GCMs to resolve mesoscale structures. We have developed a new version of the TIEGCM, which has a horizontal resolution of $\mathrm{0.625}{}^{\circ}\times \mathrm{0.625}{}^{\circ}$ in a geographic longitude–latitude grid, by implementing the ring average technique as a postprocessing step. The nonphysical “hole” and “ring” structures, which are induced by an FFT filter in the previous TIEGCM version, no longer exist in the highresolution TIEGCM associated with the ring average technique. The simulation results illustrate that the highresolution TIEGCM is capable of resolving the mesoscale structures in the I–T system during a geomagnetic storm event. Moreover, the ring average scheme has also been implemented in CMIT and WACCMX to enable highspatialresolution selfconsistent simulations of the whole geospace from the ground to the magnetosphere.
The source codes of TIEGCM 2.1 are provided through the GitHub repository (https://github.com/dangtustc/TIEGCM2.1, last access: February 2021) and Zenodo archive (https://doi.org/10.5281/zenodo.4506230, Dang, 2021b).
The simulation data used in this study are available from the Zenodo archive (https://doi.org/10.5281/zenodo.4506273, Dang, 2021a).
TD developed the model code, performed the simulations, and wrote the paper. BZ, JL, and KAS proposed the original idea and revised the paper. WW and AB helped to develop the code and edit the paper. HLL and KP contributed to coupling the highresolution TIEGCM to WACCMX and CMIT, respectively.
The authors declare that they have no conflict of interest.
We are grateful for support from the ISSI/ISSIBJ workshop “MultiScale Magnetosphere–Ionosphere–Thermosphere Interaction”. We would also like to acknowledge the Supercomputing Center of the University of Science and Technology of China for the numerical simulations in this study.
This work was supported by the Btype Strategic Priority Program of the Chinese Academy of Sciences (XDB41000000), the National Natural Science Foundation of China (41831070, 41974181), the preresearch project on Civil Aerospace Technologies no. D020105 funded by China's National Space Administration, and the Open Research Project of Large Research Infrastructures of CAS – “Study on the interaction between the low and midlatitude atmosphere and ionosphere based on the Chinese Meridian Project”. Tong Dang was supported by the National Natural Science Foundation of China (41904138), the National Postdoctoral Program for Innovative Talents (BX20180286), the China Postdoctoral Science Foundation (2018M642525), and the Fundamental Research Funds for the Central Universities. Binzheng Zhang was supported by the RGC General Research Fund (17300719, 17308520) and the Excellent Young Scientists Fund (Hong Kong and Macau) of the National Natural Science Foundation of China (41922060). Hanli Liu's work was supported by the National Center for Atmospheric Research, which is a major facility sponsored by the National Science Foundation under cooperative agreement no. 1852977. The National Center for Atmospheric Research is sponsored by the National Science Foundation.
This paper was edited by Josef Koller and reviewed by two anonymous referees.
Basu, S., Basu, S., Sojka, J. J., Schunk, R. W., and MacKenzie, E.: Macroscale modeling and mesoscale observations of plasma density structures in the polar cap, Geophys. Res. Lett., 22, 881–884, 1995. a
Bouaoudia, S. and Marcus, P.: Fast and accurate spectral treatment of coordinate singularities, J. Comput. Phys., 96, 217–223, 1991. a
Codrescu, M. V., FullerRowell, T. J., and Foster, J. C.: On the importance of Efield variability for Joule heating in the highlatitude thermosphere, Geophys. Res. Lett., 22, 2393–2396, 1995. a
Colella, P. and Woodward, P. R.: The Piecewise Parabolic Method (PPM) for gasdynamical simulations, J. Comput. Phys., 54, 174–201, 1984. a
Courant, R., Friedrichs, K., and Lewy, H.: Über die partiellen Differenzengleichungen der mathematischen Physik, Math. Annal., 100, 32–74, 1928. a
Crowley, G., Knipp, D. J., Drake, K. A., Lei, J., Sutton, E., and Lühr, H.: Thermospheric density enhancements in the dayside cusp region during strong BY conditions, Geophys. Res. Lett., 37, L07110, https://doi.org/10.1029/2009GL042143, 2010. a
Dang, T.: Supplementary Material for the paper “Azimuthal averaging–reconstruction filtering techniques for finitedifference general circulation models in spherical geometry”, Zenodo, https://doi.org/10.5281/zenodo.4506273, 2021a. a
Dang, T.: Highresolution TIEGCM code, Zenodo, https://doi.org/10.5281/zenodo.4506230, 2021b. a
Dang, T., Lei, J., Wang, W., Burns, A., Zhang, B., and Zhang, S.R.: Suppression of the polar tongue of ionization during the 21 August 2017 solar eclipse, Geophys. Res. Lett., 45, 2918–2925, 2018a. a
Dang, T., Lei, J., Wang, W., Zhang, B., Burns, A., Le, H., Wu, Q., Ruan, H., Dou, X., and Wan, W.: Global responses of the coupled thermosphere and ionosphere system to the August 2017 Great American Solar Eclipse, J. Geophys. Res.Space, 123, 7040–7050, 2018b. a
Dang, T., Lei, J., Wang, W., Wang, B., Zhang, B., Liu, J., Burns, A., and Nishimura, Y.: Formation of double tongues of ionization during the 17 March 2013 geomagnetic storm, J. Geophys. Res.Space, 124, 10619–10630, https://doi.org/10.1029/2019JA027268, 2019. a
Foster, J. C., Coster, A. J., Erickson, P. J., Holt, J. M., Lind, F. D., Rideout, W., McCready, M., van Eyken, A., Barnes, R. J., Greenwald, R. A., and Rich, F. J.: Multiradar observations of the polar tongue of ionization, J. Geophys. Res.Space, 110, A09S31, https://doi.org/10.1029/2004JA010928, 2005. a, b
Fukagata, K. and Kasagi, N.: Highly energyconservative finite difference method for the cylindrical coordinate system, J. Comput. Phys., 181, 478–498, 2002. a
Fuller‐Rowell, T. J., Rees, D., Quegan, S., Moffett, R. J., Codrescu, M. V., and Millward, G. H.: A coupled thermosphere‐ionosphere model (CTIM), in: Handbook of ionospheric models, edited by: Schunk, R. W., State Univ., Logan, Utah, 217–238, 1996. a
Hagan, M. E. and Forbes, J. M.: Migrating and nonmigrating diurnal tides in the middle and upper atmosphere excited by tropospheric latent heat release, J. Geophys. Res., 107, 4754, https://doi.org/10.1029/2001JD001236, 2002. a
Hagan, M. E. and Forbes, J. M.: Migrating and nonmigrating semidiurnal tides in the upper atmosphere excited by tropospheric latent heat release, J. Geophys. Res., 108, 1062, https://doi.org/10.1029/2002JA009466, 2003. a
Heelis, R. A., Lowell, J. K., and Spiro, R. W.: A model of the highlatitude ionospheric convection pattern, J. Geophys. Res.Space, 87, 6339–6345, 1982. a
Lei, J., Dang, T., Wang, W., Burns, A., Zhang, B., and Le, H.: Long‐lasting response of the global thermosphere and ionosphere to the 21 August 2017 solar eclipse, J. Geophys. Res.Space Physics, 123, 4309–4316, 2018. a
Lin, D., Wang, W., Scales, W. A., Pham, K., Liu, J., Zhang, B., Merkin, V., Shi, X., Kunduri, B., and Maimaiti, M.: SAPS in the 17 March 2013 Storm Event: Initial Results From the Coupled MagnetosphereIonosphereThermosphere Model, J. Geophys. Res.Space, 124, 6212–6225, 2019. a
Liu, H.L., Bardeen, C. G., Foster, B. T., Lauritzen, P., Liu, J., Lu, G., Marsh, D. R., Maute, A., McInerney, J. M., Pedatella, N. M., Qian, L., Richmond, A. ., Roble, R. G., Solomon, S. C., Vitt, F. M., and Wang, W.: Development and Validation of the Whole Atmosphere Community Climate Model With Thermosphere and Ionosphere Extension (WACCMX 2.0), J. Adv. Model. Earth Syst., 10, 381–402, 2018. a
Lotko, W. and Zhang, B.: Alfvénic Heating in the Cusp IonosphereThermosphere, J. Geophys. Res.Space, 123, 10368–10383, 2018. a
Lu, G., Zakharenkova, I., Cherniak, I., and Dang, T.: Large‐scale ionospheric disturbances during the 17 March 2015 storm: A modeldata comparative study, J. Geophys. Res.Space, 125, e2019JA027726, https://doi.org/10.1029/2019JA027726, 2020. a
Lühr, H., Rother, M., Köhler, W., Ritter, P., and Grunwaldt, L.: Thermospheric upwelling in the cusp region: Evidence from CHAMP observations, Geophys. Res. Lett., 31, L06805, https://doi.org/10.1029/2003GL019314, 2004. a
Lyon, J. G., Fedder, J. A., and Mobarry, C. M.: The Lyon–Fedder–Mobarry (LFM) global MHD magnetospheric simulation code, J. Atmos. Sol.Terr. Phy., 66, 1333–1350, 2004. a
Makela, J. J. and Otsuka, Y.: Overview of Nighttime Ionospheric Instabilities at Low and MidLatitudes: Coupling Aspects Resulting in Structuring at the Mesoscale, Space Sci. Rev., 168, 419–440, 2012. a
Matsuo, T. and Richmond, A. D.: Effects of highlatitude ionospheric electric field variability on global thermospheric Joule heating and mechanical energy transfer rate, J. Geophys. Res.Space, 113, A07309, https://doi.org/10.1029/2007JA012993, 2008. a
Prusa, J. M.: Computation at a coordinate singularity, J. Comput. Phys., 361, 331–352, 2018. a
Purser, J.: Degradation of numerical differencing caused by Fourier filtering at high latitudes, Mon. Weather Rev., 116, 1057–1066, 1988. a
Qian, L., Burns, A. G., Emery, B. A., Foster, B., Lu, G., Maute, A., Richmond, A. D., Roble, R. G., Solomon, S. C., and Wang, W.: The NCAR TIEGCM: A community model of the coupled thermosphere/ionosphere system, American Geophysical Union, Washington, D.C., 73–83, 2014. a
Ren, Z., Wan, W., and Liu, L.: GCITEMIGGCAS: A new global coupled ionosphere–thermosphereelectrodynamics model, J. Atmos. Sol.Terr. Phy., 71, 2064–2076, 2009. a
Richmond, A. D., Ridley, E. C., and Roble, R. G.: A thermosphere/ionosphere general circulation model with coupled electrodynamics, Geophys. Res. Lett., 19, 601–604, 1992. a, b
Ridley, A. J., Deng, Y., and Tóth, G.: The global ionosphere–thermosphere model, J. Atmos. Sol.Terr. Phy., 68, 839–864, 2006. a
Roble, R. G., Ridley, E. C., Richmond, A. D., and Dickinson, R. E.: A coupled thermosphere/ionosphere general circulation model, Geophys. Res. Lett., 15, 1325–1328, 1988. a
Skamarock, W. C., Klemp, J. B., Dudhia, J., Gill, D. D., Barker, D. M., Duda, M. G., Huang, X.Y., Wang, W., and Powers, J. G.: A description of the advanced research WRF version 3, NCAR Technical Note NCAR/TN475+STR, Nat. Cent. Atmos. Res., Boulder, CO, 113 pp., 2008. a
Sun, L., Xu, J., Wang, W., Yue, X., Yuan, W., Ning, B., Zhang, D., and M., D.: Mesoscale field‐aligned irregularity structures (FAIs) of airglow associated with medium‐scale traveling ionospheric disturbances (MSTIDs), J. Geophys. Res.Space, 120, 9839–9858, 2015. a
Takacs, L., Sawyer, W., Suarez, M. J., and FoxRabinowitz, M. S.: Filtering Techniques on a Stretched Grid General Circulation Model, NASA/TM1999104606, Tech. Rept. Ser., vol. 16, Global Modeling Data Assimilation, Goddard SFC, Greenbelt, MD, 1999. a
Thomas, E. G., Baker, J. B. H., Ruohoniemi, J. M., Clausen, L. B. N., Coster, A. J., Foster, J. C., and Erickson, P. J.: Direct observations of the role of convection electric field in the formation of a polar tongue of ionization from storm enhanced density, J. Geophys. Res.Space, 118, 1180–1189, 2013. a
Wang, W.: A thermosphere‐ionsopehre nested grid (ting) model, PhD thesis, University of Michigan, Ann Arbor, 1998. a
Wang, W., Wiltberger, M., Burns, A. G., Solomon, S. C., Killeen, T. L., Maruyama, N., and Lyon, J. G.: Initial results from the coupled magnetosphere–ionosphere–thermosphere model: thermosphere–ionosphere responses, J. Atmos. Sol.Terr. Phy., 66, 1425–1441, 2004. a
Wang, W., Dang, T., Lei, J., Zhang, S.R., Zhang, B., and Burns, A.: Physical Processes Driving the Response of the F2 Region Ionosphere to the 21 August 2017 Solar Eclipse at Millstone Hill, J. Geophys. Res.Space, 124, 2978–2991, 2019. a
Weimer, D. R.: Improved ionospheric electrodynamic models and application to calculating Joule heating rates, J. Geophys. Res.Space, 110, A05306, https://doi.org/10.1029/2004JA010884, 2005. a
Williamson, D. L. and Browning, G. L.: Comparison of Grids and Difference Approximations for Numerical Weather Prediction Over a Sphere, J. Appl. Meteorol., 12, 264–274, 1973. a
Williamson, D. L., Drake, J. B., Hack, J. J., Jakob, R., and Swarztrauber, P. N.: A standard test set for numerical approximations to the shallow water equations in spherical geometry, J. Comput. Phys., 102, 211–224, 1992. a
Wiltberger, M., Wang, W., Burns, A. G., Solomon, S. C., Lyon, J. G., and Goodrich, C. C.: Initial results from the coupled magnetosphere ionosphere thermosphere model: magnetospheric and ionospheric responses, J. Atmos. Sol.Terr. Phy., 66, 1411–1423, 2004. a
Wu, Q., Knipp, D., Liu, J., Wang, W., Varney, R., Gillies, R., Erickson, P., Greffen, M., Reimer, A., Häggström, I., Jee, G., and Kwak, Y.: HIWIND Observation of Summer Season Polar Cap Thermospheric Winds, J. Geophys. Res.Space, 124, 9270–9277, 2019. a
Zhang, B., Sorathia, K. A., Lyon, J. G., Merkin, V. G., and Wiltberger, M.: Conservative averagingreconstruction techniques (Ring Average) for 3D finitevolume MHD solvers with axis singularity, J. Comput. Phys., 376, 276–294, 2019. a, b, c
Zhang, Q. H., Zhang, B., Lockwood, M., Hu, H., Moen, J., Ruohoniemi, J. M., Thomas, E. G., Zhang, S.R., Yang, H., Liu, R., McWilliams, K. A., and Baker, J. B. H.: Direct Observations of the Evolution of Polar Cap Ionization Patches, Science, 339, 1597–1600, 2013. a
Zhu, Q., Deng, Y., Richmond, A., and Maute, A.: SmallScale and Mesoscale Variabilities in the Electric Field and Particle Precipitation and Their Impacts on Joule Heating, J. Geophys. Res.Space., 123, 9862–9872, 2018. a