The Global Numerical Model of the Earth’s Upper Atmosphere

Abstract

The global numerical first principle 3D model of the upper atmosphere (UAM) for the heights 60–100,000 km is presented. The physical continuity, motion, heat balance and electric potential equations for the neutral, ion and electron gases and their numerical solution method are described. The numerical grids, spatial and time integration steps are given together with the boundary and initial conditions and inputs. Testing and obtained geophysical results are given for many observed situations at various levels of solar, geomagnetic and seismic activity.

1. Introduction

The upper atmosphere is a part of the gas envelope of the Earth. It is located from the height h ≈ 60 km to several R_E of the geocentric distance, where R_E = 6371 km is the Earth’s radius. It is characterized by a sharp transition from the predominance of the neutral particles to the charged particles.

The upper atmosphere is divided into several height regions depending on its gas composition and dominating physical process: the thermosphere (from ~80–90 km to ~400–800 km) and the exosphere (above ~400–800 km) in relation to the neutral particles; the ionosphere, the plasmasphere and the magnetosphere in relation to the charged particles.

The ionosphere is located at the height range from ~50 km (the upper part of the middle atmosphere—mesosphere) to ~1000 km. The plasmasphere is located above the ionosphere up to the plasmapause—the geomagnetic force line with L~(2.5 – 7), where L is the geocentric distance to the geomagnetic line top expressed in units of R_E. The magnetosphere is the region below the magnetopause, where the solar wind pressure balances the geomagnetic field one. The magnetopause height is ~10R_E geocentric distance on the day side and ~100R_E on the night one.

The upper atmosphere state is part of Space Weather. It experiences regular annual, seasonal and diurnal variations as well as disturbances caused by the solar, geomagnetic and lithosphere activities, both globally and locally. Along with the meteorological weather, the Space Weather greatly affects human activity. Such disturbances as geomagnetic storms and substorms, auroras lead to disruptions of the radio communication in the HF range up to its blackout, faulty operation of the navigation satellite systems and electronic onboard equipment of the aircrafts. They generate intense geomagnetic-induced currents in the long conducting lines (transmission facilities, telecommunication, cables, railways and oil and gas pipelines), which create failures of the automatic and relay protection systems and, as a consequence, they cause emergency shutdowns of power supply systems, etc. Therefore, monitoring and forecasting of the upper atmosphere state is an extremely important task.

Monitoring of the near-Earth space is conducted by measurements of different physical plasma parameters (temperatures, concentrations and velocities of neutral and charged components, electric and magnetic fields, etc.) at different heights and areas of the globe, both ground-based and on the board of airplanes, rockets and satellites. Despite the growing level of the technical perfectness, it is impossible to provide stable and global monitoring, and “white areas” still remain over the Earth. These areas are those where measurements cannot be conducted or experience various difficulties (especially over the oceans and near the poles). In such cases, the method of the numerical simulation becomes valuable, and the calculations of the desired parameters fill the “white gaps”. The numerical models use the basic fundamental laws of nature to describe quantitatively the near-Earth environment and/or interpret the measurements.

The Upper Atmosphere Model (UAM), described in this chapter, is a global, 3D, self-consistent, numerical model covering h = 60 − 100,000 km. It solves the system of equations describing the basic laws of mass, energy, momentum and electric charge conservation and calculates variations of the neutral (O₂, N₂, O, H) and charged (electrons and ions O₂⁺, NO⁺, O⁺, H⁺) particles’ concentrations; temperatures of neutral gas, ions and electrons; the velocities of the particles as well as the electric field potential for any geo- and heliophysical conditions. The UAM takes into account the non-coincidence of the Earth’s geographic and geomagnetic axes.

The UAM was developed previously in Kaliningrad under the supervision of Prof. A.A. Namgaladze as the Global Self-consistent Model of the Thermosphere, Ionosphere and Protonosphere (GSM TIP) [1, 2] and further improved in Murmansk. The modern version called as UAM [3, 4] differs from the GSM TIP by implementation of algorithms for the integration with variable latitude steps and by including empirical models of the ionosphere and thermosphere to use their data for initial and boundary conditions and for the UAM testing.

Most of the modern numerical models (NCAR TIE-GCM [5], CTIM [6], CTIPe [7], GITM [8], SWMF [9, 10]) either cover the near-Earth space in very limited ranges in heights and latitudes, or perform a combination of several models without physical coupling between the parts. The UAM covers the near-Earth space as a coupled system and is still superior to all existing models by spatial coverage and resolution. This makes the UAM suitable to investigate a variety of physical processes, both globally and locally.

2. Basic equations

The fundamental conservation laws are used in the UAM: the quasi hydrodynamic equations of the continuity (1), momentum (2), heat balance (3) and electric current continuity (4).

$$\partial n_{\alpha }/\partial t+\nabla \left(n_{\alpha }{\overrightarrow{v}}_{\alpha }\right)=Q_{\alpha }-L_{\alpha },$$E1

$$\left(\partial {\overrightarrow{v}}_{\alpha }/\partial t+\left[\overrightarrow{\Omega }\times \left[\overrightarrow{\Omega }\times \overrightarrow{r}\right]\right]+2\left[\overrightarrow{\Omega }\times {\overrightarrow{v}}_{\alpha }\right]\right)=\overrightarrow{F_{\alpha }},$$E2

$$\begin{array}{c}\overrightarrow{F_{\alpha }}=-\nabla p_{\alpha }+n_{\alpha }{m}_{\alpha }\overrightarrow{g}-{\sum }_i{\mu }_{\mathit{\alpha }\mathit{i}}{V}_{\mathit{\alpha }\mathit{i}}{n}_{\alpha }\left({\overrightarrow{v}}_{\alpha }-{\overrightarrow{v}}_i\right)+e{n}_{\alpha }\left(\overrightarrow{E}+\left[{\overrightarrow{v}}_{\alpha }\times \overrightarrow{B}\right]\right)\\ {\rho }_{\alpha }{c}_{\mathit{V}\mathit{\alpha }}\partial T_{\alpha }/\partial t+p_{\alpha }\nabla {\overrightarrow{v}}_{\alpha }=\nabla \left({\lambda }_{\alpha }\nabla T_{\alpha }\right)+P_{\mathit{Q}\mathit{\alpha }}-P_{\mathit{L}\mathit{\alpha }}+P_{\mathit{T}\mathit{\alpha }},\end{array}$$E3

$$div\overrightarrow{j}=0.$$E4

where n_α is the concentration of the α-th gas, ${\overrightarrow{v}}_{\alpha }$ is its velocity vector, Q_α and L_α are its production and loss rates, respectively, ρ_α is the mass density of the α-th gas, $\overrightarrow{\mathrm{\Omega }}$ is the angular velocity of the Earth’s rotation, $\overrightarrow{r}$ is the radius-vector directed from the center of the Earth, $\overrightarrow{F_{\alpha }}$ is the force vector acting on the unit gas volume, p_α = n_αkT_α is the partial pressure of the α-th gas, T_α is its temperature, k is Boltzmann’s constant, $\overrightarrow{g}$ is the vector of the gravity acceleration, m_α is the mass of the α-th gas particle, μ_αi is the reduced mass of colliding particles, v_αi is the collision frequency, $\left({\overrightarrow{v}}_{\alpha }-{\overrightarrow{v}}_i\right)$ is the relative velocity of the α-th and i-th gases, e is the elementary charge, $\overrightarrow{E}$ and $\overrightarrow{B}$ are the electric and magnetic fields, respectively, c_Vα is the specific heat per unit volume of the α-th gas, λ_α is the thermal conductivity, P_Qα is the rate of heat of the α-th gas, P_Lα is its rate of cooling, P_Tα is the rate of heat exchange between the α-th gas and other particles and $\overrightarrow{j}$ is the electric current density.

Due to the changes in gas composition and the predominance of different physical processes at different heights, the UAM is divided into several computational blocks: the neutral atmosphere and lower ionosphere block; the F2-layer and plasmasphere block and the magnetosphere block. A special block calculates the electric field potential by solving the continuity equation for $\overrightarrow{j}$. Each block covers a particular height range, uses its own particular equations, coordinate system and calculates its own set of parameters and exchanges with other blocks at each time step of numerical integration.

The block of the neutral atmosphere and lower ionosphere uses the spherical geomagnetic coordinate system and calculates the 3D variations of n_n for O, O₂ and N₂ gases, n_e, T_n and $\overrightarrow{v}$, where ​​v → ​​ is the neutral wind velocity vector, within the height range from the model’s lower boundary to 520 km, as well as T_i and T_e up to 175 km.

Concentrations n(N₂) are obtained from the barometric law, n(O) and n(O₂) from the continuity Eq. (1), which takes view as:

$$\partial n_n/\partial t+\nabla \left[n_n\left(\overrightarrow{v}+{\overrightarrow{v}}_{\mathit{dn}}\right)\right]=Q_n-L_n,$$E5

where $\overrightarrow{v}$ is the neutral wind velocity, ${\overrightarrow{v}}_{\mathit{dn}}$ is the diffusion velocity vector, which has only a vertical component equal to the sum of molecular and turbulent diffusion velocities, Q_n and L_n are the production and loss rates, respectively, taking into account the photodissociation of the molecular oxygen by the solar radiation and recombination of O and O₂ in the triple collisions and radiative recombination.

The meridional and zonal components of $\overrightarrow{v}$ are obtained from the projection of Eq. (2) on the horizontal axes of the geographical coordinate system:

$$\rho {\left(\partial \overrightarrow{v}/\partial t+\left(\overrightarrow{v},\nabla \right)\overrightarrow{v}+2\left[\overrightarrow{\Omega }\times \overrightarrow{r}\right]\right)}_{\mathit{hor}}=-{\left(\nabla p\right)}_{\mathit{hor}}-{\overrightarrow{R}}_{\mathit{ni}\mathit{hor}}+\eta {\left({\nabla }^2\overrightarrow{\nu }\right)}_{\mathit{hor}},$$E6

where ρ is the mean neutral density, $p=\sum _n{n}_{n}k{T}_n$ is the mean partial pressure, ${\overrightarrow{R}}_{\mathit{ni}\mathit{hor}}={\sum }_n{\sum }_i{\mu }_{\mathit{ni}}{\nu }_{\mathit{ni}}{n}_n{\left(\overrightarrow{v}-{\overrightarrow{v}}_i\right)}_{\mathit{hor}}$ is the horizontal projection of the neutral-ion friction force, η is the coefficient of viscosity, n and i lowercases stand for the neutral particles and ions, respectively, ν_ni is the frequency of collisions between the neutrals and ions, $\left(\overrightarrow{v}-{\overrightarrow{v}}_i\right)$ is the relative velocity of the neutrals and ions.

Neutral density ρ is calculated from the projection of Eq. (2) on the radii vector $\overrightarrow{r}$ directed from the center of the Earth to the particular point in space, taking into account only the gravity force and vertical component of ∇p, that is, fulfilling the hydrostatic equilibrium:

$$\mathit{\rho }\mathit{g}=-\partial p/\partial r,\rho ={\sum }_n{\rho }_n={\sum }_n{n}_n{m}_n,$$E7

where g is the vertical projection of the sum of the gravitational and centrifugal accelerations.

By summing up for n the continuity Eq. (5) and assuming that the total number of particles remains constant, that is, $\sum _n{Q}_n=\sum _n{L}_n$,we obtain:

$$\partial \rho /\partial t+\nabla \left(\rho \overrightarrow{v}\right)=0.$$E8

With the calculated ρ and horizontal components of $\overrightarrow{v}$, we calculate the vertical wind velocity from Eq. (8).

Finally, T_n is calculated from the heat balance Eq. (3) written as:

$$\rho c_V\partial T_n/\partial t+\left(\overrightarrow{v},\nabla \right)T_n+p\nabla \overrightarrow{v}=\nabla \left({\lambda }_n\nabla T_n\right)+P_{\mathit{nQ}}^{\mathit{UV}}+P_{\mathit{nQ}}^J+P_{\mathit{nQ}}^C-P_{\mathit{nL}},$$E9

where c_V and λ_n are the specific heat per unit volume and thermal conductivity, respectively, $P_{\mathit{nQ}}^{\mathit{UV}}$, $P_{\mathit{nQ}}^J$ and $P_{\mathit{nQ}}^C$ are the rates of the heating by the ultraviolet (UV) and extra-ultraviolet (EUV) solar radiation, Joule heating due to collisions with ions and heating by particles precipitating from the magnetosphere, respectively; P_nL is the rate of the cooling due to radiation.

For concentrations of O₂⁺, N₂⁺ and NO⁺ ions in the D, E and F1 layers, the photochemical, local heating and heat exchange processes dominate over the transport processes, since the lifetime of the molecular ions is many times smaller than the transport characteristic time due to high collision frequencies of the charged particles with neutrals and each other. Thus, Eq. (1) for the molecular ions in the lower ionosphere block is written as:

$$\partial n_i/\partial t=\partial n\left({\mathit{XY}}^{+}\right)/\partial t=Q\left({\mathit{XY}}^{+}\right)-L\left({\mathit{XY}}^{+}\right),$$E10

where n(XY⁺) = n_e is the total concentration of the molecular ions, Q(XY⁺) and L(XY⁺) are the production and loss rates due to the ionization by direct UV solar radiation, scatter radiation and ionization by the precipitating electrons, as well as the result of the ion-molecular reactions and dissociative recombination.

The momentum equation for the molecular ions is written as:

$$n_i{m}_i\overrightarrow{g}-\nabla \left(n_{i}k{T}_i\right)-{\sum }_n{\mu }_{\mathit{\text{in}}}{\nu }_{\mathit{\text{in}}}{n}_i\left({\overrightarrow{v}}_i-\overrightarrow{v}\right)+e{n}_i\left(\overrightarrow{E}+\left[{\overrightarrow{v}}_i\times \overrightarrow{B}\right]\right)=0.$$E11

In the momentum equation for electrons, we neglect all terms containing m_e due to its smallness:

$$\nabla \left(n_{e}k{T}_e\right)+e{n}_e\left(\overrightarrow{E}+\left[{\overrightarrow{v}}_e\times \overrightarrow{B}\right]\right)=0.$$E12

This gives the electric field of the ambipolar diffusion:

$${\overrightarrow{E}}_{\mid \mid }=\nabla \left(n_{e}k{T}_e\right)/e{n}_e.$$E13

The temperatures T_i and T_e are calculated in this block taking into account that up to h = 175 km the local heating dominates, and the processes of the heating transport are neglected in the heat balance Eq. (3):

$$3/2n_{i}k\partial T_i/\partial t=P_{\mathit{iQ}}^J+P_{\mathit{iT}}^e+P_{\mathit{iT}}^n,$$E14

$$3/2n_{e}k\partial T_e/\partial t=P_{\mathit{eQ}}^P+P_{\mathit{eQ}}^C+P_{\mathit{eT}}^i+P_{\mathit{eT}}^n,$$E15

where $P_{\mathit{iQ}}^J$ is the rate of Joule’s heating, $P_{\mathit{iT}}^e=-P_{\mathit{iT}}^n$ are the rates of the ions’ heat exchange with electrons and neutrals, respectively, $P_{\mathit{eT}}^i$ and $P_{\mathit{eT}}^n$ are the rates of the electrons’ heat exchange with ions and neutrals, respectively, $P_{\mathit{eQ}}^P$ and $P_{\mathit{eQ}}^C$ are the rates of the electrons’ heating by the photoelectrons and electrons precipitating from the magnetosphere, respectively.

In the block of the F2-layer and plasmasphere n(O⁺), n(H⁺), $\overrightarrow{v}$_i, T_i and T_e are calculated for h = 175 − 100,000 km. At these heights collision frequencies ν_in of charged particles with neutrals are much smaller than ions’ gyrofrequencies Ω_i, that is, collisions do not interfere the cyclotron rotation and drift of charged particles. Ions and electrons are fully magnetized, that is, they are tied to the geomagnetic field lines and can move across the lines only under the action of an extraneous force. Because the geomagnetic field has such great influence on the behavior of ions and electrons, we use the magnetic dipole coordinate system in this block. The system of modeling Eqs. (1)–(3), thus, is integrated along the field lines taking into account the electromagnetic plasma drift perpendicular to the magnetic field.

In the continuity equation, along with the production rates Q_i and loss rates L_i of ions O⁺ and H⁺ by photoionization, corpuscular ionization, chemical reactions between O⁺, O₂ and N₂, as well as the charge exchange between O⁺ and H and between O and H⁺, the transport of ionized components (the divergence of charged particles’ flow) is taken into account:

$$D{n}_i/\mathit{Dt}+{\nabla }^{\mathit{par}}\left(n_i{v}_i^{\mathit{par}}\right)=Q_i-L_i-n_i{\nabla }^{\mathit{per}}{v}^{\mathit{per}}.$$E16a

Here, superscripts par and per refer to directions parallel and perpendicular to the geomagnetic field lines, respectively. The operator

$$D/\mathit{Dt}=\partial /\partial t+\left(v_i^{\mathit{per}},\nabla \right),$$E16b

gives the Lagrangian time derivatives along the trajectory of the charged particle’s electromagnetic drift perpendicular to the geomagnetic field line with the velocity:

$$v_i^{\mathit{per}}=v_e^{\mathit{per}}=\left[\overrightarrow{E}\times \overrightarrow{B}\right]/B^2.$$E16c

The electrons’ velocity along the field lines is calculated as:

$$v_e^{\mathit{par}}=\sum _i{n}_i{v}_i^{\mathit{par}}/n_e.$$E17

The atomic ions’ motion equation for the i-th atomic ion is given as:

$$\begin{array}{c}2m_i{n}_i{\left(\overrightarrow{\Omega }\times {\overrightarrow{v}}_i\right)}^{par}=m_i{n}_i{g}^{par}-{\nabla }^{par}\left(n_{i}k{T}_i\right)-n_i/n_e{\nabla }^{par}\left(n_{e}k{T}_e\right)\\ -{\displaystyle \sum }_n{\mu }_{in}{\nu }_{in}{n}_i\left(v_i^{par}-v^{par}\right)-{\displaystyle \sum }_j{\mu }_{ij}{\nu }_{ij}{n}_i\left(v_j^{par}-v^{par}\right),\end{array}$$E18

where subscripts i, j, e and n refer to O⁺, H⁺, electrons and neutrals, respectively and g^par is the projection of the sum of the gravity and centrifugal accelerations along the geomagnetic field line.

The temperatures T_i and T_e are calculated from the heat balance equations written as:

$$3/2n_{i}k\left(D{T}_i/\mathit{Dt}+v_i^{\mathit{par}}{\nabla }^{\mathit{par}}{T}_i\right)+n_{i}k{T}_i\nabla v_i^{\mathit{par}}-{\nabla }^{\mathit{par}}\left({\lambda }_i{\nabla }^{\mathit{par}}{T}_i\right)=P_{\mathit{iQ}}^J+P_{\mathit{iT}}^e+P_{\mathit{iT}}^j+P_{\mathit{iT}}^n,$$E19

$$3/2n_{e}k\left(D{T}_e/\mathit{Dt}+v_e^{\mathit{par}}{\nabla }^{\mathit{par}}{T}_e\right)+n_{e}k{T}_e\nabla v_e^{\mathit{par}}-{\nabla }^{\mathit{par}}\left({\lambda }_e{\nabla }^{\mathit{par}}{T}_e\right)=P_{\mathit{eT}}^i+P_{\mathit{eT}}^j+P_{\mathit{eT}}^n+P_{\mathit{eQ}}^P+P_{\mathit{eQ}}^C,$$E20

where $P_{\mathit{iQ}}^J$ is the rate of Joule’s heating for the ion gas; $P_{\mathit{iT}}^j$, $P_{\mathit{iT}}^e$ and $P_{\mathit{iT}}^n$ are the rates of the heat exchange between ions, between ions and electrons and between ions and neutrals, respectively; $P_{\mathit{eT}}^i$, $P_{\mathit{eT}}^j$ and $P_{\mathit{eT}}^n$ are the rates of the electrons’ heat exchange with the i-th and j-th ions as well as the neutral gas, respectively; $P_{\mathit{eQ}}^P$ and $P_{\mathit{eQ}}^C$ are the rates of the electrons’ heating by the photoelectrons and electrons precipitating from the magnetosphere, respectively.

The magnetospheric block of the UAM calculates the plasma layer n_i from the continuity Eq. (21), ${\overrightarrow{v}}_i$ from the motion Eq. (22), ion pressure p_i from Eq. (23) and the field-aligned current density ${\overrightarrow{j}}_{\mid \mid }$ (FACs) from Eq. (24) above 175 km:

$$\partial n_i/\partial t+\overrightarrow{\nabla }\left(n_i{\overrightarrow{v}}_i\right)=0,$$E21

$$e{n}_i\left(\overrightarrow{E}+\left[{\overrightarrow{v}}_i\times \overrightarrow{B}\right]\right)=\overrightarrow{\nabla }{p}_i,$$E22

$$d\left(p_i{V}^{\gamma }\right)/\mathit{dt}=0,$$E23

$${\overrightarrow{j}}_{\mid \mid }=\frac{{\overrightarrow{e}}_z\left[\overrightarrow{\nabla }V\times \overrightarrow{\nabla }{p}_i\right]}{B},$$E24

where V is a half-volume of the geomagnetic field tube, γ = 5/3 is the adiabatic coefficient for the plasma layer ions and ${\overrightarrow{e}}_z$ is the unit vector along $\overrightarrow{B}$. V is calculated as:

$$V=B{\int }_0^{z_{\mathit{max}}}{B}^{-1}\mathit{dz},$$E25

where z is a distance along $\overrightarrow{B}$, z_max is the distance to the top of the geomagnetic field line.

We assume that the pressure of the plasma layer is isotropic and constant along $\overrightarrow{B}$. The pressure of the plasma layer electrons is neglible in comparison with p_i, that is, the electrons are considered to be cold.

In the block of the electric field potential ϕ, it is calculated taking into account electric fields of magnetospheric, thermospheric (dynamo), and lower atmospheric (due to the charges transfer into the ionosphere from bellow) origins. The equation of the electric current continuity is written as:

$$\nabla \overrightarrow{j}=\nabla \left({\overrightarrow{j}}_i+{\overrightarrow{j}}_m+{\overrightarrow{j}}_s\right)=0,$$E26a

where ${\overrightarrow{j}}_m$ and ${\overrightarrow{j}}_s$ are the densities of the magnetosphere and lower ionosphere electric currents, and ${\overrightarrow{j}}_i$ is the ionosphere electric current density given by Ohm’s law for the plasma:

$${\overrightarrow{j}}_i={\sigma }_{\mathit{par}}{\overrightarrow{E}}_{\mathit{par}}+{\sigma }_P{\overrightarrow{E}}_{\mathit{per}}+{\sigma }_H\left[\overrightarrow{B}\times {\overrightarrow{E}}_{\mathit{per}}\right].$$E26b

Here, σ_P and σ_H are the Pedersen and Hall conductivities, respectively, and σ_par is the conductivity along $\overrightarrow{B}$.

$${\sigma }_P=e^2{n}_e\left[\frac{{\nu }_{\mathit{\text{in}}}}{m_i\left({\mathrm{\Omega }}_i^2+{\nu }_{\mathit{\text{in}}}^2\right)}+\frac{{\nu }_{\mathit{en}}}{m_e\left({\mathrm{\Omega }}_e^2+{\nu }_{\mathit{en}}^2\right)}\right],$$E26c

$${\sigma }_H=e^2{n}_e\left[\frac{{\mathrm{\Omega }}_e}{m_e\left({\mathrm{\Omega }}_e^2+{\nu }_{\mathit{en}}^2\right)}-\frac{{\mathrm{\Omega }}_i}{m_i\left({\mathrm{\Omega }}_i^2+{\nu }_{\mathit{\text{in}}}^2\right)}\right],$$E26d

$${\sigma }_{\mathit{par}}=e^2{n}_e\left[\frac{1}{m_i{\nu }_{\mathit{\text{in}}}}+\frac{1}{m_e{\nu }_{\mathit{en}}}\right],$$E26e

where Ω_i and Ω_e are the ion and electron cyclotron frequencies, respectively.

$${\overrightarrow{E}}_{\mathit{par}}=\overrightarrow{B}\left(\overrightarrow{B},\overrightarrow{E}\right)/B^2,$$E26f

$${\overrightarrow{E}}_{\mathit{per}}=\overrightarrow{B}\times \left[\overrightarrow{B}\times \overrightarrow{E}\right]/B^2,$$E26g

where $\overrightarrow{E}$ is the electric field vector defined as the sum of the electrostatic field with the potential ϕ and the induced dynamo field:

$$\overrightarrow{E}=-\nabla \varphi +\left[\overrightarrow{v}\times \overrightarrow{B}\right].$$E26h

Thus, the equation for the electric current continuity is given as:

$$-\nabla \left(\widehat{\sigma }\overrightarrow{E}+{\overrightarrow{j}}_m+{\overrightarrow{j}}_s\right)=\nabla \left(\widehat{\sigma }\left(\nabla \varphi -\left[\overrightarrow{v}\times \overrightarrow{B}\right]\right)-{\overrightarrow{j}}_m-{\overrightarrow{j}}_s\right)=0,$$E27a

where $\widehat{\sigma }$ is the tensor of the ionospheric electric conductivity for the coordinate system with axes along $\overrightarrow{B}$, $\overrightarrow{E}$ and $\left[\overrightarrow{E}\times \overrightarrow{B}\right]$:

$$\widehat{\sigma }=\left(\begin{array}{ccc}{\sigma }_P& 0& {\sigma }_H\\ 0& {\sigma }_{\mathit{par}}& 0\\ -{\sigma }_H& 0& {\sigma }_P\end{array}\right).$$E27b

The integration of Eq. (27a) is conducted along the electric conducting layer, from the lower boundary of the UAM up to 175 km, where we neglect the dependency of the electric field intensity on h. Above 175 km, we assume that the plasma is magnetized and the geomagnetic field lines are electrically equipotential.

3. Equations solutions: initial and boundary conditions

The integration of the model equations is carried out using a finite-difference numerical method. The near-Earth environment is considered as a discrete three-dimensional grid, and each computational block uses its own coordinate system. The derivative operators are replaced by the difference ratios, and the solution is obtained in the nodes of the numerical grid. Steps in the numerical grid are chosen depending on the task and the characteristic scale of the investigated process. Usually, the integration step in longitude is chosen as a constant value between 5° and 15°. The integration steps in latitude are variable and depend on the latitude. Near the geomagnetic equator the step is chosen in the range of 2.5–5.0° and 1° near the poles, because a tighter numerical grid is required for a precise calculation in the aurora zones and across the polar caps. Steps of the numerical integration along the height are also variable. In the blocks where a spherical coordinate system is used, the step is 3 km at the lower boundary of the UAM and increases exponentially with altitude.

The solutions of Eqs. (1)–(4), which are containing derivations of coordinates, require boundary conditions: the distribution of the desired parameters at the boundaries of the numerical blocks. At the lower boundary T_n and n_n are obtained from the empirical model NRLMSISE-00 (further simply MSIS [11]), $\overrightarrow{v}$ is calculated from Eqs. (6) and (8) neglecting the viscosity and ion-neutral friction, in the geostrophic approximation, when Eq. (6) contains only the Coriolis acceleration and the pressure gradient of the neutral gas.

For the upper boundary conditions in the block of the neutral atmosphere (at 520 km) the diffusion equilibrium for n_n is used, $\overrightarrow{v}$ and T_n are considered to be independent of the height:

$$\partial n_n/\partial r+m_{n}g/\mathit{kT}=0,$$E28

$$\partial v/\partial r=0,$$E29

$$\partial T_n/\partial r=0.$$E30

In the F2-layer and plasmasphere block the boundary conditions are defined as follows. At the altitude of 520 km, the concentration of the neutral hydrogen is set according to the empirical model of the neutral atmosphere [12]. For atomic ions n_i at the bases of the geomagnetic field lines at 175 km, the production rates are equal to the loss rates:

$$Q_i=L_i.$$E31

Eqs. (14) and (15) are used for the T_i and T_e at 175 km.

The model equations are integrated along the geomagnetic field lines in the areas with closed lines of force (±75° in geomagnetic latitude). In the regions with open geomagnetic field lines (poleward of 75° geomagnetic latitude) n_i, as well as the heat fluxes of the ions and electrons are set equal to zero at the upper boundary of the UAM. Thus, the model reproduces the condition of the polar wind: a supersonic outflow of plasma from the F2-layer and upper ionosphere of the polar caps along the open lines of the geomagnetic field.

Solutions of those equations, which are containing time derivatives, require initial conditions: parameter distributions at the start of the numerical calculation. For quiet geophysical conditions, stationary solutions are used as initial conditions after the multiple runs of the model. The initial conditions for perturbed periods are the solutions obtained for the previous quiet days. It is also possible to use data from empirical models (NRLMSISE-00 [11], HWM-93 [13], IRI-2001 [14] or IRI-2007 [15]) as initial conditions.

4. Model inputs

The inner state of the simulated space and external forcings acting on it are characterized by the model inputs set up by the user: (1) date and time to set an initial placement of the numerical grid nodes relative to the Sun; (2) spectra of the solar UV and EUV radiation; (3) fluxes of the high energetic electrons precipitating from the magnetosphere; (4) the FACs connecting the ionosphere with the magnetosphere and/or (5) the distribution of ϕ at the boundaries of the polar cap; (6) the local ${\overrightarrow{j}}_s$ flowing through the lower boundary from below; (7) indices of the geomagnetic activity and (8) components of the interplanetary magnetic field (IMF) and solar wind.

The solar UV and EUV spectra define the coefficients of O₂ dissociation and O₂⁺, N₂⁺, NO⁺ and O⁺ production rates due to the photoionization of the corresponding neutral components. The UV and EUV spectra dependence on the solar activity is set up according to Ref. [16]. The intensity of the night scatter radiation intensity is 5 kR for the emission with wave length λ = 121.6 nm and 5 R for the rest emission lines (102.6, 58.4 and 30.4 nm).

The precipitating electrons’ fluxes are set up at the upper boundary of the thermosphere, at 520 km, and their intensity I is written as:

$$I\left(\mathrm{\Phi },\mathrm{\Lambda },E\right)=I_m\left(E\right)\mathit{exp}\left[-{\left(\mathrm{\Phi }-{\mathrm{\Phi }}_m\left(E\right)\right)}^2/\mathrm{\Delta \Phi }{\left(E\right)}^2-{\left(\mathrm{\Lambda }-{\mathrm{\Lambda }}_m\left(E\right)\right)}^2/\mathrm{\Delta \Lambda }{\left(E\right)}^2\right],$$E32

$${\mathrm{\Phi }}_m=\left({\mathrm{\Phi }}_{\mathit{md}}+{\mathrm{\Phi }}_{\mathit{mn}}\right)/2+cos\mathrm{\Lambda }\left({\mathrm{\Phi }}_{\mathit{md}}+{\mathrm{\Phi }}_{\mathit{mn}}\right)/2,$$E33

where Φ and Λ are the geomagnetic latitude and longitude, respectively (Λ = 0 corresponds to the midday magnetic meridian); E is the energy of the precipitating electrons; I_m(E) is the maximum intensity of the precipitating electron flux; ∆Φ and ∆Λ are the half-widths of the maximum precipitations in latitude and longitude; Φ_md and Φ_mn are the magnetic latitudes of the maximum precipitations at the midday and midnight meridians, respectively. Specific values for the precipitating parameters in Eqs. (32) and (33) are taken from the empirical models [17, 18].

The magnetospheric sources ${\overrightarrow{j}}_m$ in Eq. (27a) specify the distribution of FACs. The FACs of the Region 1 (R1) flow from the magnetosphere into the ionosphere on the dawn side and out of the ionosphere on the dusk side, at latitudes higher ±75°. The FACs of the Region 2 (R2) flow in the opposite direction, at areas equatorward of the R1 currents. The distribution of current densities depends on the geophysical conditions and is setup in the UAM in several different ways depending on the task, either as distribution of the FACs in the R1, 2 and the cusp region according to the model [19], or as the distribution of electric potential at the boundary of the polar cap [20] with the FACs in the cusp region and the R2.

The so-called seismogenic electric currents are the vertical electric currents switched on to simulate the ionosphere effects of various mesoscale phenomena in the lower ionosphere, such as earthquakes, thunderstorms, etc. Used as a model input, the vertical ${\overrightarrow{j}}_s$ are added locally to Eq. (27a) at the lower boundary of the UAM.

The geomagnetic activity is used in the UAM by setting up the planetary geomagnetic indices Kp and Ap, indicating global geomagnetic field disturbance, and aurora indices AE, AL and AU, where AU and AL indicate, respectively, the largest increase and lowest decrease of the Northern component of the geomagnetic field in comparison to the background (quiet) value, AE is the sum of AL and AU and characterizes the largest scale of the magnetic field during the substorm in the high-latitude regions.

5. The UAM versions

The UAM provides the possibility of integrating various empirical models and data of the upper atmosphere. The comparison of the self-consistent UAM version and the UAM versions with different combinations of the empirical models allows testing both the UAM and the empirical models.

In the UAM-MSIS version, n_n and T_n are calculated directly from the MSIS [11]. The thermospheric circulation is calculated by the numerical solution of Eqs. (6) and (8) where ∇p from the MSIS is used. Finally, the MSIS is used to set up T_n and n_n at the lower boundary and initial conditions as well.

In the UAM-HWM version the distribution of the horizontal thermospheric wind is calculated using the empirical model HWM-93 [13]. The vertical component of the wind velocity is calculated by the numerical solution of the continuity equation for ρ (Eq. (8)). The HWM-93 is used for the set of the initial conditions and for comparison of the theoretical model winds with observations.

The magnetospheric block of the UAM simulates the transport processes in the plasma sheet by solving the system of the equations for the plasma sheet ions (see item 2). In Ref. [21], the initial values are taken as p_i = 0, 4 nPa and n_i = 0, 4 cm⁻³, correspondingly. The program produces more or less realistic p_i distribution and R2 FACs. The problem is that the obtained solution is not stable, and it falls apart after approximately 1 h.

There are several ways to set up FACs spatial-temporal distributions in the UAM, such as empirical data from the magnetic field measurements from the Dynamics Explorer 2 [22] and the Magsat satellites [23], the FACs empirical models by Papitashvili [24] in [25, 26], by Lukianova [27] in Ref. [28] and MFACE [29] in Ref. [30]. All these versions with various FACs take into account the dependence of FACs on the interplanetary magnetic field (IMF). Such methods of setting the FACs distribution allow using any other empirical data of FACs.

In the UAM version [31], the positions of the auroral oval boundaries, the values of electron flux intensities and the latitudes and longitudes of the intensity maxima were set from precipitation patterns observed by DMSP. The spectra of the precipitating electrons are assumed to be Maxwellian in this case.

The UAM-P version [32, 33], created in Potsdam, differs from the UAM by the electric field block simulation. This block uses magnetic dipole coordinates instead of spherical geographical ones within h = 80 − 526 km. It is assumed that $\overrightarrow{E}$ does not change along $\overrightarrow{B}$ inside the ionospheric current layer. After the integration, it is also represented as a 2D-distribution. It allows to exclude the conductivity parallel to the magnetic field and to keep the vertical electric field inside the current-carrying layer. As a result, the lower latitudinal and equatorial electric field distributions as well as the current system of these areas look more correctly.

The Canadian Ionosphere and Atmosphere Model (Canadian IAM or C-IAM) is comprised of the extended Canadian Middle Atmosphere (CMAM) and the UAM, currently coupled in a one-way manner [34]. This version was used to investigate the physical mechanisms responsible for forming the four-peak longitudinal structure of the 135.6 nm ionospheric emission observed by the IMAGE satellite over the tropics at 20:00 local time from March 20 to April 20, 2002. The study showed that main mechanism is driven by the diurnal eastward propagating tide with zonal number 3.

6. Results

During its development and improvement, the UAM was used to perform a number of numerical experiments aimed at testing and comparison of calculation results with measurements and other models. The UAM successfully showed its ability to reproduce the general behavior of ionospheric and thermospheric parameters such as at low and high geomagnetic and solar activity conditions. A good agreement of the numerical simulations’ results was achieved in comparison with observations by incoherent scatter radars located at various latitudes and longitudes (ISRs) [31], digital ionosonde CADI in the Voeykovo Main Geophysical Observatory [35], several chains of the ionospheric tomographical receivers [36, 37, 38], satellites, including CHAMP and GPS [3, 4, 26, 30, 39, 40, 41, 42], as well as with empirical models such as various types of IRI, MSIS, HWM [42, 43, 44, 45, 46, 47], etc. and other theoretical models [42].

Numerical modeling of the upper atmosphere behavior during substorms. These relatively short disturbances of the geomagnetic field have duration from about 0.5 to 3 h. They are generated at the magnetopause and in the magnetosphere tail via the magnetic field lines reconnection processes and connected with the polar upper atmosphere in the auroral zones via FACs including the current wedges. The substorm auroral currents are reflected by the auroral magnetic activity indexes AL, AU and AE. The UAM takes into account these indexes. This allows the modeling of the upper atmospheric behavior during substorms via the UAM simulations. The results were presented in [4, 48, 49, 50, 51, 52, 53] including the cusp and auroral oval behavior, energetic magnetospheric electron precipitations, electric fields, current wedge and internal atmospheric gravity waves generation. The main role of the thermospheric heating due to the soft electron precipitation was shown for the thermosphere substorm effects.

Numerical modeling of the upper atmosphere behavior during geomagnetic storms. These geomagnetic disturbances are global and have duration of about several days. Usually they include several substorms on the background of the geomagnetic field depression created by the ring current development. The geomagnetic activity indexes Dst and Kp characterize the geomagnetic storms. These indexes are included in the UAM inputs as well. This allows to model the upper atmospheric behavior during magnetic storms via the UAM simulations. The calculation results were presented in [25, 26, 30, 36, 37, 38, 54, 55, 56, 57, 58, 59] including the main ionospheric trough dynamics due to the ionosphere-magnetosphere convection and non-coincidence of the geographical and geomagnetic axes of the Earth. The physical mechanisms of the negative and positive F2-layer ionospheric storms (electron concentration decreases and increases, correspondingly) formation were described. The main role of the thermospheric composition (atomic to molecular neutral gas concentration ratio) and winds disturbances in the magnetic storm ionospheric effects was demonstrated in these UAM calculations.

In addition to the UAM testing and comparison with the observations for the different levels of solar and geomagnetic activities, the model has been widely used to study the ionosphere response on the local sources in the lower atmosphere, such as disturbances associated with the processes of the earthquakes’ preparation processes [26, 60, 61, 62, 63]. Numerical UAM calculations showed that the electric fields of 5–10 mV/m effects on the F2-layer plasma by the electromagnetic drift in the crossed geomagnetic field and the electric field of the seismic origin. The vertical electric current, flowing through the lower boundary of the ionosphere with the density of ~20 nA/m², is required for the generation of the seismogenic electric field of ~10 mV/m [60, 61, 62, 63, 64]. The important role of the aerosols over the tectonic faults was underlined in this process due to the very low recombination rate for the charged aerosols.

Thus, the UAM was tested and used in many helio- and geophysical situations. Nevertheless, the amount of the UAM simulations remains to be insufficient despite the IT progress. This is related with the specifics of the geophysics as science at all. The near-Earth space environment varies due to the solar, seismic and human activities. This does not allow performing the repeated fixed experiments as in usual physical laboratories to obtain correctly the standard statistical error estimates. Moreover, the observations themselves are very limited. None of them has 3D spatial and time resolutions satisfying to the requirements of the modern technical means of the Space weather practical usage. This was well known long ago [65] and such models as UAM are aimed solving this important problem.

7. Conclusions

Further development of the UAM means a huge amount of further numerical experiments to its mathematical and physical quality. These experiments have to take into account all modern achievements of the numerical mathematics and computer science. The numerical grids, steps, various iterations, etc. should be tested to find their optimal combinations for the best stability and accuracy. The user’s manuals should be constructed, including the UAM website. The UAM prognostic features have to be improved by modeling many case studies for various helio- and geophysical situations especially for geomagnetic and seismogenic disturbances. Comparisons with ground-based and satellite observations, empirical and other theoretical models have to be made continuously. The frame approach should be widely used by including separate observational, empirical and theoretical blocks into the UAM, such as the real geomagnetic field, polar wind, plasma sheet, electric fields, lower atmosphere, aerosols, tides, etc. An international cooperation is absolutely necessary for these future scientific UAM perspectives.

Acknowledgments

We thank many people who have worked with the UAM, developed and improved it, such as A. N. Namgaladze, M. A. Volkov, E. N. Doronina, Yu. V. Romanovskaya, I. V. Korableva, Yu. A. Shapovalova, M. G. Botova, M. I. Rybakov, I. V. Artamonov, E. V. Vasilieva, V. A. Medvedeva, V. A. Shlykov, L. A. Chernyuk and many others.