## 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.

### Keywords

- global numerical model
- UAM
- equations
- solution method
- grid
- integration steps
- input
- output
- upper atmosphere
- thermosphere
- ionosphere
- plasmasphere
- magnetosphere
- space weather
- solar and geomagnetic 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 *RE* of the geocentric distance, where *RE* = 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 *RE*. The magnetosphere is the region below the magnetopause, where the solar wind pressure balances the geomagnetic field one. The magnetopause height is ~10*RE* geocentric distance on the day side and ~100*RE* 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_{2}, N_{2}, O, H) and charged (electrons and ions O_{2}^{+}, 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).

where *nα* is the concentration of the α-th gas, *Qα* and *Lα* are its production and loss rates, respectively, *ρα* is the mass density of the α-th gas, *pα* = *nαkTα* is the partial pressure of the α-th gas, *Tα* is its temperature, *k* is Boltzmann’s constant, *mα* is the mass of the α-th gas particle, *μαi* is the reduced mass of colliding particles, *vαi* is the collision frequency, *i*-th gases, *e* is the elementary charge, *cVα* is the specific heat per unit volume of the α-th gas, *λα* is the thermal conductivity, *PQα* is the rate of heat of the α-th gas, *PLα* is its rate of cooling, *PTα* is the rate of heat exchange between the α-th gas and other particles and

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

The block of the neutral atmosphere and lower ionosphere uses the spherical geomagnetic coordinate system and calculates the 3D variations of *nn* for O, O_{2} and N_{2} gases, *ne*, *Tn* and *Ti* and *Te* up to 175 km.

Concentrations *n*(N_{2}) are obtained from the barometric law, *n*(O) and *n*(O_{2}) from the continuity Eq. (1), which takes view as:

where *Qn* and *Ln* 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_{2} in the triple collisions and radiative recombination.

The meridional and zonal components of

where *ρ* is the mean neutral density, *η* 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,

Neutral density *ρ* is calculated from the projection of Eq. (2) on the radii vector *p*, that is, fulfilling the hydrostatic equilibrium:

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,

With the calculated *ρ* and horizontal components of

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

where *cV* and *λn* are the specific heat per unit volume and thermal conductivity, respectively, *PnL* is the rate of the cooling due to radiation.

For concentrations of O_{2}^{+}, N_{2}^{+} 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:

where *n*(*XY*^{+}) = *ne* 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:

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

This gives the electric field of the ambipolar diffusion:

The temperatures *Ti* and *Te* 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):

where

In the block of the F2-layer and plasmasphere *n*(O^{+}), *n*(H^{+}), _{i}, *Ti* and *Te* 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 *Li* of ions O^{+} and H^{+} by photoionization, corpuscular ionization, chemical reactions between O^{+}, O_{2} and N_{2}, 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:

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

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

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

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

where subscripts *i*, *j*, *e* and *n* refer to O^{+}, H^{+}, electrons and neutrals, respectively and *gpar* is the projection of the sum of the gravity and centrifugal accelerations along the geomagnetic field line.

The temperatures *Ti* and *Te* are calculated from the heat balance equations written as:

where *i*-th and *j*-th ions as well as the neutral gas, respectively;

The magnetospheric block of the UAM calculates the plasma layer *ni* from the continuity Eq. (21), *pi* from Eq. (23) and the field-aligned current density

where *V* is a half-volume of the geomagnetic field tube, *γ* = 5/3 is the adiabatic coefficient for the plasma layer ions and *V* is calculated as:

where *z* is a distance along *zmax* 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 *pi*, 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:

where

Here, *σP* and *σH* are the Pedersen and Hall conductivities, respectively, and *σpar* is the conductivity along

where Ω*i* and Ω*e* are the ion and electron cyclotron frequencies, respectively.

where *ϕ* and the induced dynamo field:

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

where

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 *Tn* and *nn* are obtained from the empirical model NRLMSISE-00 (further simply MSIS [11]),

For the upper boundary conditions in the block of the neutral atmosphere (at 520 km) the diffusion equilibrium for *nn* is used, *Tn* are considered to be independent of the height:

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 *ni* at the bases of the geomagnetic field lines at 175 km, the production rates are equal to the loss rates:

Eqs. (14) and (15) are used for the *Ti* and *Te* 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) *ni*, 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

The solar UV and EUV spectra define the coefficients of O_{2} dissociation and O_{2}^{+}, N_{2}^{+}, 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:

where *Φ* and *Λ* are the geomagnetic latitude and longitude, respectively (*Λ* = 0 corresponds to the midday magnetic meridian); *E* is the energy of the precipitating electrons; *Im*(*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

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

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, *nn* and *Tn* 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 *Tn* and *nn* 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 *pi* = 0, 4 nPa and *ni* = 0, 4 cm^{−3}, correspondingly. The program produces more or less realistic *pi* 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

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^{2}, 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.