**1. Introduction**

The upper atmosphere is a part of the gas envelope of the Earth. It is located from the height *h* ≈ 60km 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.

Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. © 2018 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

© 2016 The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons

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.

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 investi-

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

→

<sup>→</sup>]] + 2[*Ω*

→ *<sup>α</sup>* − *v* →

→*α*

is the force vector acting on the unit gas volume, *pα* = *nαkTα* is the partial pressure of the α-th

*mα* is the mass of the α-th gas particle, *μαi* is the reduced mass of colliding particles, *Vαi* is the

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

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

coordinate system and calculates its own set of parameters and exchanges with other blocks

The block of the neutral atmosphere and lower ionosphere uses the spherical geomagnetic coor-

*ρα cV* <sup>∂</sup>*Tα*/∂*<sup>t</sup>* <sup>+</sup> *<sup>p</sup>α*∇*<sup>v</sup>*

<sup>→</sup> × *v* →

*<sup>α</sup>* = ∇(*λα*∇*Tα*) + *PQ* − *PL* + *PT*

*<sup>i</sup>*) + *e nα*(*E*

<sup>→</sup> is the radius-vector directed from the center of the Earth, →*F<sup>α</sup>*

) is the relative velocity of the α-th and *i*-th gases, *e* is the elementary

<sup>→</sup> are the electric and magnetic fields, respectively, *cVα* is the specific heat per

<sup>→</sup>. Each block covers a particular height range, uses its own particular equations,

for O, O2

and N2

gases, *ne*

, *Tn* and *v*

<sup>→</sup>, where

<sup>→</sup> <sup>+</sup> [*<sup>v</sup>* → *<sup>α</sup>* × *B* → ])

<sup>→</sup> = 0. (4)

is its velocity vector, *Qα* and *Lα* are its pro-

<sup>→</sup> is the vector of the gravity acceleration,

*<sup>α</sup>*) = *Q<sup>α</sup>* − *Lα*, (1)

The Global Numerical Model of the Earth's Upper Atmosphere

http://dx.doi.org/10.5772/intechopen.71139

5

*<sup>α</sup>*]) <sup>=</sup> <sup>→</sup>*Fα*, (2)

, (3)

<sup>→</sup> is the angular

gate a variety of physical processes, both globally and locally.

∂*nα*/∂*t* + ∇(*n<sup>α</sup> v*

div*j*

where *nα* is the concentration of the α-th gas, *<sup>v</sup>*

→ *<sup>α</sup>* − *v* → *i*

at each time step of numerical integration.

dinate system and calculates the 3D variations of *nn*

gas, *Tα* is its temperature, *k* is Boltzmann's constant, *<sup>g</sup>*

<sup>→</sup> is the electric current density.

→

<sup>→</sup>*F<sup>α</sup>* <sup>=</sup> <sup>−</sup>∇*p<sup>α</sup>* <sup>+</sup> *<sup>n</sup><sup>α</sup> <sup>m</sup>α<sup>g</sup>*

*<sup>α</sup>*/∂*t* + [*Ω*

<sup>→</sup> × [*Ω* <sup>→</sup> × *r*

→

duction and loss rates, respectively, *ρα* is the mass density of the α-th gas, Ω

<sup>→</sup> − ∑*<sup>i</sup> μ<sup>i</sup> V<sup>i</sup> nα*(*v*

**2. Basic equations**

(∂*v*

velocity of the Earth's rotation, *r*

collision frequency, (*v*

<sup>→</sup> and *<sup>B</sup>*

charge, *E*

particles and *j*

equation for *j*

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 (O2 , N2 , O, H) and charged (electrons and ions O2 + , NO<sup>+</sup> , O<sup>+</sup> , H<sup>+</sup> ) 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 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

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

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

ticles' 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.

, O, H) and charged (electrons and ions O2

+ , NO<sup>+</sup> , O<sup>+</sup> , H<sup>+</sup> ) par-

forecasting of the upper atmosphere state is an extremely important task.

night one.

4 Numerical Simulations in Engineering and Science

the measurements.

variations of the neutral (O2

, N2

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 \boldsymbol{n}\_a/\partial \boldsymbol{t} + \nabla \left( \boldsymbol{n}\_a \vec{\boldsymbol{v}}\_a \right) = \boldsymbol{Q}\_a - \boldsymbol{L}\_{\boldsymbol{n}'} \tag{1}
$$

$$\left(\partial\vec{\boldsymbol{v}}\_{a'}\big|\partial\boldsymbol{t}+\left[\vec{\boldsymbol{\Omega}}\times\left[\vec{\boldsymbol{\Omega}}\times\vec{\boldsymbol{r}}\right]\right]+\boldsymbol{\mathcal{D}}\left[\vec{\boldsymbol{\Omega}}\times\vec{\boldsymbol{v}}\_{a}\right]\right)=\overrightarrow{F}\_{a'}\tag{2}$$

$$\left(\partial\vec{v}\_{a}/\partial t + \left[\vec{\mathcal{Q}}\times\left[\vec{\mathcal{Q}}\times\vec{r}\right]\right] + 2\left[\vec{\mathcal{Q}}\times\vec{v}\_{a}\right]\right) = \overrightarrow{F}\_{a'} \tag{2}$$

$$\begin{split} \vec{F}\_{a} &= -\nabla p\_{a} + \mathfrak{n}\_{a}\mathfrak{m}\_{a}\vec{\mathcal{S}} - \sum\_{i} \mu\_{ai}\,\nu\_{ai}\,\mathfrak{n}\_{a}(\vec{\mathcal{v}}\_{a} - \vec{\mathcal{v}}\_{i}) + c\,\mathfrak{n}\_{a}(\vec{\mathcal{E}} + \left[\vec{\mathcal{v}}\_{a}\times\vec{\mathcal{B}}\right]) \\ \rho\_{a}\,c\_{\sqrt{a}}\,\partial T\_{a'}/\partial t + p\_{a}\,\nabla\vec{\mathcal{v}}\_{a} &= \nabla(\lambda\_{a}\nabla T\_{a}) + P\_{Qa} - P\_{La} + P\_{Td} \end{split} \tag{3}$$

$$\text{div}\vec{j} \tag{4}$$

where *nα* is the concentration of the α-th gas, *<sup>v</sup>* →*α* is its velocity vector, *Qα* and *Lα* are its production and loss rates, respectively, *ρα* is the mass density of the α-th gas, Ω <sup>→</sup> is the angular velocity of the Earth's rotation, *r* <sup>→</sup> is the radius-vector directed from the center of the Earth, →*F<sup>α</sup>* 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, *<sup>g</sup>* <sup>→</sup> 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, (*v* → *<sup>α</sup>* − *v* → *i* ) is the relative velocity of the α-th and *i*-th gases, *e* is the elementary charge, *E* <sup>→</sup> and *<sup>B</sup>* <sup>→</sup> are the electric and magnetic fields, respectively, *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 *j* <sup>→</sup> 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 *j* <sup>→</sup>. 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 *nn* for O, O2 and N2 gases, *ne* , *Tn* and *v* <sup>→</sup>, where

*v* <sup>→</sup> is the neutral wind velocity vector, within the height range from the model's lower boundary to 520 km, as well as *Ti* and *Te* up to 175 km.

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

$$\left[\partial n\_{n'} \partial t + \nabla \left[n\_{n} \left(\vec{\boldsymbol{\upsilon}} + \vec{\boldsymbol{\upsilon}}\_{dn}\right)\right] \right] = Q\_n - L\_{n'} \tag{5}$$

where *cV* and *λ<sup>n</sup>*

where *n*(*XY*<sup>+</sup>

*ni mi*

The temperatures *Ti*

heat balance Eq. (3):

where *PiQ*

For concentrations of O2

∂*ni*

) = *ne*

tions and dissociative recombination.

*g*

∇(*ne k Te*

*E*

and *Te*

3/2 *ni k*∂*Ti*

3/2 *ne k*∂*Te*

electrons and neutrals, respectively, *PeT*

with ions and neutrals, respectively, *PeQ*

*<sup>J</sup>* is the rate of Joule's heating, *PiT*

<sup>→</sup> − ∇(*ni k Ti*

This gives the electric field of the ambipolar diffusion:

+ , N2 +

*PnQ <sup>J</sup>* and *PnQ* are the specific heat per unit volume and thermal conductivity, respectively, *PnQ*

ions in the D, E and F1 layers, the photochemical, local

The Global Numerical Model of the Earth's Upper Atmosphere

http://dx.doi.org/10.5772/intechopen.71139

) and *L*(*XY*<sup>+</sup>

]) = 0. (11)

]) = 0. (12)

. (13)

*<sup>n</sup>* , (14)

*<sup>n</sup>* , (15)

are the rates of the ions' heat exchange with

*<sup>n</sup>* are the rates of the electrons' heat exchange

*<sup>C</sup>* are the rates of the electrons' heating by the

due to its smallness:

) are the

/∂*t* = ∂*n*(*XY*<sup>+</sup>)/∂*t* = *Q*(*XY*<sup>+</sup>) − *L*(*XY*<sup>+</sup>), (10)

*<sup>C</sup>* 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

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

is the total concentration of the molecular ions, *Q*(*XY*<sup>+</sup>

) − ∑*<sup>n</sup> μin νin ni*(*v*

) + *e ne*(*E*

∣∣ = ∇(*ne k Te*

the local heating dominates, and the processes of the heating transport are neglected in the

/∂*t* = *PiQ*

/∂*t* = *PeQ*

*<sup>e</sup>* <sup>=</sup> <sup>−</sup>*PiT n*

*<sup>i</sup>* and *PeT*

*<sup>P</sup>* and *PeQ*

photoelectrons and electrons precipitating from the magnetosphere, respectively.

In the momentum equation for electrons, we neglect all terms containing *me*

→

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 reac-

> → *<sup>i</sup>* − *v* →

<sup>→</sup> <sup>+</sup> [*<sup>v</sup>* → *<sup>e</sup>* × *B* →

)/*e ne*

*<sup>J</sup>* + *PiT <sup>e</sup>* + *PiT*

*<sup>C</sup>* + *PeT*

*<sup>P</sup>* + *PeQ*

are calculated in this block taking into account that up to *h* = 175 km

*<sup>i</sup>* + *PeT*

) + *e ni*(*E*

<sup>→</sup> <sup>+</sup> [*<sup>v</sup>* → *<sup>i</sup>* × *B* →

the magnetosphere, respectively; *PnL* is the rate of the cooling due to radiation.

and NO<sup>+</sup>

the molecular ions in the lower ionosphere block is written as:

The momentum equation for the molecular ions is written as:

*UV*,

7

where *v* <sup>→</sup> is the neutral wind velocity, <sup>→</sup>*<sup>v</sup> dn* is the diffusion velocity vector, which has only a vertical component equal to the sum of molecular and turbulent diffusion velocities, *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 O2 in the triple collisions and radiative recombination.

The meridional and zonal components of *v* <sup>→</sup> are obtained from the projection of Eq. (2) on the horizontal axes of the geographical coordinate system:

$$\rho \left( \partial \vec{v} / \partial t + (\vec{v} \,\nabla) \vec{v} + 2 \left[ \vec{\Omega} \times \vec{r} \right] \right)\_{\text{hor}} = - (\nabla p)\_{\text{hor}} - \vec{R}\_{\text{nithor}} + \eta \left( \nabla^2 \vec{v} \right)\_{\text{hor}'} \tag{6}$$

where *ρ* is the mean neutral density, *<sup>p</sup>* <sup>=</sup> <sup>∑</sup>*<sup>n</sup> nn <sup>k</sup> Tn* is the mean partial pressure, *R* → *nihor* = ∑*<sup>n</sup>* ∑*<sup>i</sup> μni <sup>ν</sup>ni nn* (*<sup>v</sup>* <sup>→</sup> − *v* →*i* )*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, (*v* <sup>→</sup> − *v* →*i* ) is the relative velocity of the neutrals and ions.

Neutral density *ρ* is calculated from the projection of Eq. (2) on the radii vector *r* <sup>→</sup> 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:

$$
\rho g = -\partial p/\partial r, \rho = \sum\_{n} \rho\_{n} = \sum\_{n} n\_{n} m\_{n} \tag{7}
$$

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, ∑*<sup>n</sup> Qn* <sup>=</sup> <sup>∑</sup>*<sup>n</sup> Ln* ,we obtain:

$$
\partial \rho / \partial t + \nabla (\rho \vec{v}) = 0. \tag{8}
$$

With the calculated *ρ* and horizontal components of *v* <sup>→</sup>, we calculate the vertical wind velocity from Eq. (8).

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

$$\rho \ c\_{\vee} \partial T\_{\underline{u}} / \partial t + \left(\vec{v}, \nabla\right) T\_{\underline{u}} + p \nabla \vec{v} = \nabla \left(\lambda\_{\underline{u}} \nabla T\_{\underline{u}}\right) + P\_{\underline{u}Q}^{\perp \underline{v}} + P\_{\underline{u}Q}^{\parallel} + P\_{\underline{u}Q}^{\mathbb{C}} - P\_{\underline{u}\mathcal{L}} \tag{9}$$

where *cV* and *λ<sup>n</sup>* are the specific heat per unit volume and thermal conductivity, respectively, *PnQ UV*, *PnQ <sup>J</sup>* and *PnQ <sup>C</sup>* 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; *PnL* is the rate of the cooling due to radiation.

For concentrations of O2 + , N2 + and NO<sup>+</sup> 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 \mathfrak{u}\_{\neq} \partial t = \partial \mathfrak{u}(XY^{\circ}) / \partial t = \mathbb{Q}(XY^{\circ}) - L(XY^{\circ}),
\tag{10}
$$

where *n*(*XY*<sup>+</sup> ) = *ne* is the total concentration of the molecular ions, *Q*(*XY*<sup>+</sup> ) and *L*(*XY*<sup>+</sup> ) 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:

*v*

where *v*

*<sup>ν</sup>ni nn* (*<sup>v</sup>* <sup>→</sup> − *v* →*i* )*hor*

to 520 km, as well as *Ti*

Concentrations *n*(N2

ity Eq. (1), which takes view as:

6 Numerical Simulations in Engineering and Science

and radiative recombination.

*ρ* (∂*v*

neutrals and ions.

from Eq. (8).

Finally, *Tn*

and *Te*

∂*nn* /∂*t* + ∇[*nn*(*v*

The meridional and zonal components of *v*

horizontal axes of the geographical coordinate system:

<sup>→</sup>, ∇)*v*

frequency of collisions between the neutrals and ions, (*v*

<sup>→</sup>/∂*t* + (*v*

where *ρ* is the mean neutral density, *<sup>p</sup>* <sup>=</sup> <sup>∑</sup>*<sup>n</sup> nn <sup>k</sup> Tn*

remains constant, that is, ∑*<sup>n</sup> Qn* <sup>=</sup> <sup>∑</sup>*<sup>n</sup> Ln*

*ρ cV*∂*Tn* /∂*t* + (*v*

∂*ρ*/∂*t* + ∇(*ρv*

With the calculated *ρ* and horizontal components of *v*

up to 175 km.

molecular oxygen by the solar radiation and recombination of O and O2

<sup>→</sup> + 2[*Ω* <sup>→</sup> × *r*

Neutral density *ρ* is calculated from the projection of Eq. (2) on the radii vector *r*

force and vertical component of ∇*p*, that is, fulfilling the hydrostatic equilibrium:

<sup>→</sup> is the neutral wind velocity vector, within the height range from the model's lower boundary

) are obtained from the barometric law, *n*(O) and *n*(O2

<sup>→</sup> + *v* →

the production and loss rates, respectively, taking into account the photodissociation of the

cal component equal to the sum of molecular and turbulent diffusion velocities, *Qn*

<sup>→</sup> is the neutral wind velocity, <sup>→</sup>*<sup>v</sup> dn* is the diffusion velocity vector, which has only a verti-

<sup>→</sup>])*hor* = −(∇*p*)*hor* − *R*

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

the center of the Earth to the particular point in space, taking into account only the gravity

*g* = −∂*p*/∂*r*, *ρ* = ∑*<sup>n</sup> ρ<sup>n</sup>* = ∑*<sup>n</sup> nn mn* , (7)

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

<sup>→</sup> = ∇(*λ<sup>n</sup>*

∇*Tn*) + *PnQ*

*UV* + *PnQ*

*<sup>J</sup>* + *PnQ*

,we obtain:

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

<sup>→</sup>, ∇) *Tn* + *p*∇*v*

*dn*)] = *Qn* − *Ln*

) from the continu-

and *Ln*

in the triple collisions

*hor*, (6)

*nihor* = ∑*<sup>n</sup>* ∑*<sup>i</sup>*

<sup>→</sup> directed from

*μni*

→

) is the relative velocity of the

are

, (5)

<sup>→</sup> are obtained from the projection of Eq. (2) on the

*ni hor* + *η* (∇<sup>2</sup>

is the mean partial pressure, *R*

*ν* →)

<sup>→</sup>) = 0. (8)

<sup>→</sup>, we calculate the vertical wind velocity

*<sup>C</sup>* − *PnL*, (9)

→

<sup>→</sup> − *v* →*i*

$$m\_i m\_j \vec{\mathbf{g}} - \nabla \{ n\_i k \, T\_i \} - \sum\_n \mu\_{ia} \nu\_{ia} n\_i \left( \vec{\boldsymbol{v}}\_i - \vec{\boldsymbol{v}} \right) + e \ln \left( \vec{\boldsymbol{E}} + \left[ \vec{\boldsymbol{v}}\_i \times \vec{\boldsymbol{B}} \right] \right) \tag{11}$$

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

$$\nabla(n\_{\epsilon}k\,T\_{\epsilon}) + e\, n\_{\epsilon}\left(\vec{E} + \left[\vec{v}\_{\epsilon} \times \vec{B}\right]\right) = 0. \tag{12}$$

This gives the electric field of the ambipolar diffusion:

$$
\vec{E}\_{\parallel} = \nabla \langle \mathbf{n}\_{\epsilon} k \, T\_{\epsilon} \rangle / e \, \mathbf{n}\_{\epsilon}. \tag{13}
$$

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):

$$\Im \mathcal{D} \, n\_{\,\,l} k \partial \, T\_{\,\,l} \prime \partial t \,\, = \, P\_{\,\,l}^{l} + P\_{\,\,l}^{\epsilon} + P\_{\,\,l}^{u} \,\tag{14}$$

$$\Im \slash \mathcal{D} \, n\_e k \partial T\_e / \partial t \,\tag{15} = P\_{eQ}^p + P\_{eQ}^c + P\_{eT}^l + P\_{eT}^u \tag{15}$$

where *PiQ <sup>J</sup>* is the rate of Joule's heating, *PiT <sup>e</sup>* <sup>=</sup> <sup>−</sup>*PiT n* are the rates of the ions' heat exchange with electrons and neutrals, respectively, *PeT <sup>i</sup>* and *PeT <sup>n</sup>* are the rates of the electrons' heat exchange with ions and neutrals, respectively, *PeQ <sup>P</sup>* and *PeQ <sup>C</sup>* 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<sup>+</sup> ), *n*(H<sup>+</sup> ), *v* → 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 *Ω<sup>i</sup>* , 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*<sup>i</sup> and loss rates *Li* of ions O<sup>+</sup> and H<sup>+</sup> by photoionization, corpuscular ionization, chemical reactions between O<sup>+</sup> , O2 and N2 , as well as the charge exchange between O<sup>+</sup> and H and between O and H<sup>+</sup> , the transport of ionized components (the divergence of charged particles' flow) is taken into account:

 *D ni* /*Dt* + ∇*par*(*ni vi par*) = *Qi* − *Li* − *ni* ∇*per vper*. (16a)

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

$$D/Dt = \left\| \langle \partial t + \langle v \rangle^{\circ r}, \nabla \rangle\_{\prime} \right\| \tag{16b}$$

where *PiQ*

tively; *PeT*

Eq. (21), *v* →*i*

density *j* → ∣∣

*<sup>i</sup>* , *PeT*

*<sup>j</sup>* and *PeT*

well as the neutral gas, respectively; *PeQ*

∂*ni*

*e ni*(*E*

*j*

the plasma layer ions and *e*

where *z* is a distance along *B*

∇*j*

*j*

*σ<sup>P</sup>* = *e* <sup>2</sup> *ne*

considered to be cold.

written as:

where *j* → *<sup>m</sup>* and *<sup>j</sup>* → *s*

tivity along →*<sup>B</sup>* .

and *j* → *i*

*d*(*pi V<sup>γ</sup>*

*V* = *B*

→*z*

*<sup>J</sup>* is the rate of Joule's heating for the ion gas; *PiT*

*j* , *PiT e* and *PiT n*

*<sup>n</sup>* are the rates of the electrons' heat exchange with the *i*-th and *j*-th ions as

exchange between ions, between ions and electrons and between ions and neutrals, respec-

*<sup>P</sup>* and *PeQ*

/∂*t* + ∇ → (*ni v* →

<sup>→</sup> <sup>+</sup> [*<sup>v</sup>* → *<sup>i</sup>* × *B* → ]) = ∇ →

where *V* is a half-volume of the geomagnetic field tube, *γ* = 5/3 is the adiabatic coefficient for

∫

0 *zmax*

We assume that the pressure of the plasma layer is isotropic and constant along →*<sup>B</sup>* . The pressure of the plasma layer electrons is neglible in comparison with *pi* , that is, the electrons are

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

→ ∣∣ <sup>=</sup> *<sup>e</sup>* → *<sup>z</sup>*[∇ → *V* × ∇ → *<sup>p</sup>* \_\_\_\_\_\_\_\_\_\_*i*]

is the unit vector along *B*

<sup>→</sup> = ∇(*j* → *<sup>i</sup>* + *j* → *<sup>m</sup>* + *j* →

[

→

*<sup>i</sup>* = *σpar E* →

is the ionosphere electric current density given by Ohm's law for the plasma:

*par* + *σ<sup>P</sup> E* →

Here, *σP* and *σ<sup>H</sup>* are the Pedersen and Hall conductivities, respectively, and *σpar* is the conduc-

*<sup>ν</sup>* \_\_\_\_\_\_\_\_\_ *in mi* (Ω*<sup>i</sup>* <sup>2</sup> + *νin* 2

photoelectrons and electrons precipitating from the magnetosphere, respectively.

The magnetospheric block of the UAM calculates the plasma layer *ni*

from the motion Eq. (22), ion pressure *pi*

(FACs) from Eq. (24) above 175 km:

are the rates of the heat

9

from the continuity

*<sup>C</sup>* are the rates of the electrons' heating by the

The Global Numerical Model of the Earth's Upper Atmosphere

http://dx.doi.org/10.5772/intechopen.71139

from Eq. (23) and the field-aligned current

*<sup>i</sup>*) = 0, (21)

)/*dt* = 0, (23)

*<sup>B</sup>* , (24)

*B*−1 *dz*, (25)

*<sup>s</sup>*) = 0, (26a)

*per*]. (26b)

, (26c)

<sup>→</sup>. *V* is calculated as:

<sup>→</sup>, *zmax* is the distance to the top of the geomagnetic field line.

are the densities of the magnetosphere and lower ionosphere electric currents,

<sup>→</sup> × *E* →

) <sup>+</sup> *<sup>ν</sup>* \_\_\_\_\_\_\_\_\_ *en me* (Ω*<sup>e</sup>* <sup>2</sup> + *νen* <sup>2</sup> )]

*per* + *σH*[*B*

*pi* , (22)

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

$$\left. \upsilon\_{l}^{\nu \upsilon} \right| = \left. \upsilon\_{\iota}^{\nu \upsilon} \right| = \left[ \vec{E} \times \vec{B} \right] / B^{2}. \tag{16c}$$

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

$$
\upsilon\_e^{\nu w} = \sum\_{l} \mathbf{n}\_l \upsilon\_l^{\nu w} / \mathbf{n}\_e. \tag{17}
$$

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

$$\begin{split} \left(\mathbf{2}\ m\_{i}\mathbf{n}\_{i}\left(\stackrel{\cdot}{\mathbf{\widehat{D}}}\times\stackrel{\cdot}{\mathbf{\widehat{v}}}\right)^{\mu\nu} = \ m\_{i}\boldsymbol{n}\_{i}\mathbf{g}^{\mu\nu} - \nabla^{\mu\nu}\{\mathbf{n}\_{i}\boldsymbol{k}\cdot\boldsymbol{T}\_{i}\} - \boldsymbol{n}\_{i}\{\boldsymbol{n}\_{e}\,\nabla^{\mu\nu}\{\mathbf{n}\_{e}\,\boldsymbol{k}\,\boldsymbol{T}\_{e}\} \\ - \sum\_{\boldsymbol{n}}\mu\_{i\boldsymbol{n}}\boldsymbol{\nu}\_{i\boldsymbol{n}}\,\boldsymbol{n}\_{i}\mathbf{n}\_{i}\{\mathbf{\nu}^{\mu\nu}\_{i} - \boldsymbol{\nu}^{\mu\nu}\} - \sum\_{\boldsymbol{l}}\mu\_{\boldsymbol{q}}\boldsymbol{\nu}\_{\boldsymbol{q}}\,\boldsymbol{n}\_{\boldsymbol{q}}\{\mathbf{\nu}^{\mu\nu}\_{\boldsymbol{q}} - \boldsymbol{\nu}^{\mu\nu}\}. \end{split} \tag{18}$$

where subscripts *i*, *j*, *e* and *n* refer to O<sup>+</sup> , H<sup>+</sup> , 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:

$$\Im \mathcal{J} \mathcal{D} \mathfrak{u}\_{\!} \not\! \mathcal{D} \left( \mathcal{D} \not\!\! \mathcal{T}\_{\!\!\!/ } \mathcal{D} \!\!\! t + \mathcal{v}\_{\!\!\!\! \! \mathcal{V}} \mathcal{V}^{\mu \nu} \!\!\!\!\!\!\!\!\!\!\!\/) + \mathfrak{u}\_{\!\!\!\!\!} \mathcal{k} \!\!\!\!\!\!\!\!\!\/) \mathfrak{v}\_{\!\!\!\!\!\!\!\/) \mathcal{B} \!\!\!\!\!\!\/) \mathcal{B} \!\!\!\!\!\!\!\/) \mathcal{B} \!\!\!\!\!\!\!\/) \mathcal{B} \!\!\!\!\!\!\!\/) \mathcal{B} \!\!\!\!\!\!\/) \mathcal{B} \!\!\!\!\!\!\/) \mathcal{B} \!\!\!\!\!\/) \mathcal{B} \!\!\!\!\/) \mathcal{B} \!\!\!\!\/) \mathcal{B} \!\!\!\!\/) \mathcal{B} \!\!\!\/) \mathcal{B} \!\!\!\/) \mathcal{B} \!\!\!\/) \mathcal{B} \!\!\!\/) \mathcal{B} \!\!\!\/) \mathcal{B} \!\!\!\/) \mathcal{B} \!\!\!\/) \mathcal{B} \!\!\/) \mathcal{B} \!\!\/) \mathcal{B} \!\!\/) \mathcal{B} \!\!\/) \mathcal{B} \!\!\/) \mathcal{B} \!\!\/) \mathcal{B} \!\!\/) \mathcal{B} \!\!\/) \mathcal{B} \!\!\/) \mathcal{B} \!\!\/) \mathcal{B} \!\!\/) \mathcal{B} \!\!\/) \mathcal{B} \!\!\/) \mathcal{B} \!\!\/) \mathcal{B} \!\!\/$$

$$\mathcal{Z}/\mathcal{Z}\ n\_{\boldsymbol{\epsilon}}k\Big[\mathcal{D}\ T\_{\boldsymbol{\epsilon}}/\mathcal{D}t+\boldsymbol{\upsilon}\_{\boldsymbol{\epsilon}}^{\boldsymbol{u}\boldsymbol{v}}\,\nabla^{\boldsymbol{u}\boldsymbol{v}}\,T\_{\boldsymbol{\epsilon}}\Big]+\boldsymbol{\eta}\_{\boldsymbol{\epsilon}}k\,T\_{\boldsymbol{\epsilon}}\nabla\boldsymbol{\upsilon}\_{\boldsymbol{\epsilon}}^{\boldsymbol{u}\boldsymbol{v}}-\nabla^{\boldsymbol{u}\boldsymbol{v}}\Big(\mathcal{A}\_{\boldsymbol{\epsilon}}\,\nabla^{\boldsymbol{u}\boldsymbol{v}}\,T\_{\boldsymbol{\epsilon}}\Big)=\mathcal{P}\_{\boldsymbol{\epsilon}\mathcal{T}}^{\boldsymbol{i}}+\mathcal{P}\_{\boldsymbol{\epsilon}\mathcal{T}}^{\boldsymbol{i}}+\mathcal{P}\_{\boldsymbol{\epsilon}\mathcal{T}}^{\boldsymbol{u}}+\mathcal{P}\_{\boldsymbol{\epsilon}\mathcal{Q}}^{\boldsymbol{o}}+\mathcal{P}\_{\boldsymbol{\epsilon}\mathcal{Q}'}^{\boldsymbol{o}}\tag{20}$$

where *PiQ <sup>J</sup>* is the rate of Joule's heating for the ion gas; *PiT j* , *PiT e* and *PiT n* are the rates of the heat exchange between ions, between ions and electrons and between ions and neutrals, respectively; *PeT <sup>i</sup>* , *PeT <sup>j</sup>* and *PeT <sup>n</sup>* are the rates of the electrons' heat exchange with the *i*-th and *j*-th ions as well as the neutral gas, respectively; *PeQ <sup>P</sup>* and *PeQ <sup>C</sup>* 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<sup>+</sup>

are much smaller than ions' gyrofrequencies *Ω<sup>i</sup>*

8 Numerical Simulations in Engineering and Science

plasma drift perpendicular to the magnetic field.

well as the charge exchange between O<sup>+</sup>

netic field lines, respectively. The operator

*vi*

*ve*

2 *mi ni* (*Ω*

where subscripts *i*, *j*, *e* and *n* refer to O<sup>+</sup>

and *Te*

/*Dt* + *vi*

/*Dt* + *ve*

*D*/*Dt* = ∂/∂*t* + (*vi*

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

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

− ∑*<sup>n</sup> μin νin ni*

*par* ∇*par Ti*

*par* ∇*par Te*

<sup>→</sup> × *v* → *i* ) *par*

, H<sup>+</sup>

) + *ni k Ti*

) + *ne k Te*

projection of the sum of the gravity and centrifugal accelerations along the geomagnetic field

∇ *vi*

∇*ve*

*D ni*

H<sup>+</sup>

line.

The temperatures *Ti*

3/2 *ni k*(*D Ti*

3/2 *ne k*(*D Te*

In the continuity equation, along with the production rates *Q*<sup>i</sup>

by photoionization, corpuscular ionization, chemical reactions between O<sup>+</sup>

components (the divergence of charged particles' flow) is taken into account:

/*Dt* + ∇*par*(*ni vi*

magnetic drift perpendicular to the geomagnetic field line with the velocity:

*per* = *ve*

), *n*(H<sup>+</sup>

and H and between O and H<sup>+</sup>

*h* = 175 − 100,000 km. At these heights collision frequencies *νin* of charged particles with neutrals

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

Here, superscripts *par* and *per* refer to directions parallel and perpendicular to the geomag-

gives the Lagrangian time derivatives along the trajectory of the charged particle's electro-

*per* = [*E*

*par* = ∑*<sup>i</sup> ni vi*

<sup>→</sup> × *B* <sup>→</sup>]/*B*<sup>2</sup>

*par*/*ne*

= *mi ni gpar* − ∇*par*(*ni k Ti*

are calculated from the heat balance equations written as:

*par* − ∇*par*(*λ<sup>i</sup>* ∇*par Ti*

*par* − ∇*par*(*λ<sup>e</sup>* ∇*par Te*) = *PeT*

(*vi*

), *v* → i , *Ti*

and *Te*

, that is, collisions do not interfere the cyclotron

and loss rates *Li*

*par*) = *Qi* − *Li* − *ni* ∇*per vper*. (16a)

*per*, ∇), (16b)

. (16c)

. (17)

/*ne* ∇*par*(*ne k Te*

)

*par* − *vpar*), (18)

*<sup>n</sup>* , (19)

*<sup>C</sup>* , (20)

) − *ni*

, electrons and neutrals, respectively and *gpar* is the

) = *PiQ*

*<sup>i</sup>* + *PeT*

*<sup>J</sup>* + *PiT <sup>e</sup>* + *PiT <sup>j</sup>* + *PiT*

*<sup>j</sup>* + *PeT*

*<sup>n</sup>* + *PeQ*

*<sup>P</sup>* + *PeQ*

*par* − *vpar*) − ∑*<sup>j</sup> μij νij ni*(*vj*

are calculated for

of ions O<sup>+</sup>

and N2

, O2

, the transport of ionized

and

, as

The magnetospheric block of the UAM calculates the plasma layer *ni* from the continuity Eq. (21), *v* →*i* from the motion Eq. (22), ion pressure *pi* from Eq. (23) and the field-aligned current density *j* → ∣∣ (FACs) from Eq. (24) above 175 km:

$$
\partial \mathfrak{u}\_{\flat} / \partial \mathfrak{t} + \vec{\nabla} \left( \mathfrak{u}\_{\flat} \vec{\boldsymbol{\upsilon}}\_{i} \right) = \mathfrak{0},\tag{21}
$$

$$e\, n\_i \left(\vec{E} + \left[\vec{v}\_i \times \vec{B}\right]\right) = \vec{\nabla} \, p\_{i\cdot\cdot} \tag{22}$$

$$d(p\_\perp V^\gamma)/dt = 0,\tag{23}$$

$$
\vec{f}\_{\parallel} = \frac{\vec{e}\_{\perp} \left[ \mathcal{V} \, V \times \mathcal{V} \, p\_{\parallel} \right]}{B} \, \tag{24}
$$

where *V* is a half-volume of the geomagnetic field tube, *γ* = 5/3 is the adiabatic coefficient for the plasma layer ions and *e* →*z* is the unit vector along *B* <sup>→</sup>. *V* is calculated as:

$$V = B \int\_0^{\frac{\pi}{\omega}} B^{-1} dz,\tag{25}$$

where *z* is a distance along *B* <sup>→</sup>, *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 →*<sup>B</sup>* . The pressure of the plasma layer electrons is neglible in comparison with *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:

$$
\vec{\nabla}\_l^\dagger = \nabla \left( \vec{j}\_l + \vec{j}\_m + \vec{j}\_s \right) = 0,\tag{26a}
$$

where *j* → *<sup>m</sup>* and *<sup>j</sup>* → *s* are the densities of the magnetosphere and lower ionosphere electric currents, and *j* → *i* is the ionosphere electric current density given by Ohm's law for the plasma:

$$
\vec{f}\_{\perp} = \sigma\_{\mu\nu} \vec{E}\_{\mu\nu} + \sigma\_p \vec{E}\_{\mu\nu} + \sigma\_H \left[ \vec{B} \times \vec{E}\_{\mu\nu} \right]. \tag{26b}
$$

Here, *σP* and *σ<sup>H</sup>* are the Pedersen and Hall conductivities, respectively, and *σpar* is the conductivity along →*<sup>B</sup>* .

$$\sigma\_p = e^2 n\_e \left[ \frac{\nu\_{\rm in}}{m \langle \Omega\_i^2 + \nu\_{\rm in}^2 \rangle} + \frac{\nu\_m}{m \langle \Omega\_i^2 + \nu\_m^2 \rangle} \right] \tag{26c}$$

$$
\sigma\_H = e^2 \,\eta\_e \left[ \frac{\Omega\_\epsilon}{m\_\epsilon (\Omega\_\epsilon^2 + \nu\_m^2)} - \frac{\Omega\_i}{m\_i (\Omega\_i^2 + \nu\_m^2)} \right] \,\tag{26d}
$$

$$
\sigma\_{\mu\nu} = \, ^\mu \pi\_\epsilon \left[ \frac{1}{m\_i \, \nu\_{in}} + \frac{1}{m\_i \, \nu\_m} \right] \,\tag{26e}
$$

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

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

ion-neutral friction, in the geostrophic approximation, when Eq. (6) contains only the Coriolis

For the upper boundary conditions in the block of the neutral atmosphere (at 520 km) the

∂*nn* /∂*r* + *mn g*/*kT* = 0, (28)

∂*v*/∂*r* = 0, (29)

∂*Tn* /∂*r* = 0. (30)

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

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

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],

<sup>→</sup> and *Tn*

are obtained from the empirical model NRLMSISE-00

The Global Numerical Model of the Earth's Upper Atmosphere

http://dx.doi.org/10.5772/intechopen.71139

11

are considered to be independent of the height:

at the bases of the geomagnetic field

. (31)

, as well as the heat fluxes of the ions and

<sup>→</sup> is calculated from Eqs. (6) and (8) neglecting the viscosity and

lower boundary of the UAM and increases exponentially with altitude.

and *nn*

acceleration and the pressure gradient of the neutral gas.

model of the neutral atmosphere [12]. For atomic ions *ni*

*Qi* = *Li*

field lines (poleward of 75° geomagnetic latitude) *ni*

IRI-2001 [14] or IRI-2007 [15]) as initial conditions.

Eqs. (14) and (15) are used for the *Ti*

lines at 175 km, the production rates are equal to the loss rates:

and *Te*

blocks. At the lower boundary *Tn*

diffusion equilibrium for *nn* is used, *<sup>v</sup>*

(further simply MSIS [11]), *v*

where Ω*<sup>i</sup>* and Ω*<sup>e</sup>* are the ion and electron cyclotron frequencies, respectively.

$$
\vec{E}\_{\mu\nu} = \vec{B}(\vec{B}, \vec{E}) / B^2 \,. \tag{26f}
$$

$$
\vec{E}\_{per} = \vec{B} \times \left[\vec{B} \times \vec{E}\right] / B^2 \,\tag{26g}
$$

where *E* <sup>→</sup> is the electric field vector defined as the sum of the electrostatic field with the potential *ϕ* and the induced dynamo field:

$$
\vec{E} = -\nabla\varphi + \left|\vec{v} \times \vec{B}\right|.\tag{26h}
$$

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

$$-\nabla \left(\vec{\sigma}\vec{E} + \vec{j}\_m + \vec{j}\_s\right) = \nabla \left(\partial \left(\nabla \varphi - \left[\vec{\upsilon} \times \vec{B}\right]\right) - \vec{j}\_m - \vec{j}\_s\right) = 0,\tag{27a}$$

where *σ*̂ is the tensor of the ionospheric electric conductivity for the coordinate system with axes along *B* <sup>→</sup>, *<sup>E</sup>* <sup>→</sup> and [*<sup>E</sup>* <sup>→</sup> × *B* →]:

$$
\boldsymbol{\hat{o}} = \begin{pmatrix} \sigma\_p & 0 & \sigma\_H \\ 0 & \sigma\_{\mu\nu} & 0 \\ -\sigma\_H & 0 & \sigma\_p \end{pmatrix} . \tag{27b}
$$

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]), *v* <sup>→</sup> 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 *nn* is used, *<sup>v</sup>* <sup>→</sup> and *Tn* are considered to be independent of the height:

$$
\partial m\_{\mu}/\partial r + m\_{\mu}g/kT = 0,\tag{28}
$$

$$
\partial v/\partial r = 0,\tag{29}
$$

$$
\partial T\_{\mu}/\partial r = 0.\tag{30}
$$

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:

$$\mathbf{Q}\_{\parallel} = \mathbf{L}\_{\mathbf{r}}.\tag{31}$$

Eqs. (14) and (15) are used for the *Ti* and *Te* at 175 km.

*<sup>σ</sup><sup>H</sup>* <sup>=</sup> *<sup>e</sup>* <sup>2</sup> *ne*[

10 Numerical Simulations in Engineering and Science

*E*

*E*

tial *ϕ* and the induced dynamo field:

*E*

<sup>→</sup> × *B* →]:

*σ*̂ =

field lines are electrically equipotential.

−∇(*σ*̂*E*

<sup>→</sup>, *<sup>E</sup>* <sup>→</sup> and [*<sup>E</sup>*

where Ω*<sup>i</sup>*

where *E*

axes along *B*

and Ω*<sup>e</sup>*

*σpar* = *e* <sup>2</sup> *ne*[

\_\_\_\_\_\_\_\_\_ <sup>Ω</sup>*<sup>e</sup> me* (Ω*<sup>e</sup>* <sup>2</sup> + *νen*

> \_\_\_\_ 1 *mi νin*

are the ion and electron cyclotron frequencies, respectively.

<sup>→</sup> × [*B* <sup>→</sup> × *E* <sup>→</sup>]/*B*<sup>2</sup>

<sup>→</sup> = −∇*ϕ* + [*v*

*<sup>s</sup>*) = ∇(*σ*̂

⎛ ⎜ ⎝

<sup>→</sup> is the electric field vector defined as the sum of the electrostatic field with the poten-

(∇*ϕ* − [*v*

*σ<sup>P</sup>* 0 *σ<sup>H</sup>* 0 *σpar* 0 −*σ<sup>H</sup>* 0 *σ<sup>P</sup>*

where *σ*̂ is the tensor of the ionospheric electric conductivity for the coordinate system with

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

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

<sup>→</sup> × *B*

<sup>→</sup> × *B* <sup>→</sup>]) − *j* → *<sup>m</sup>* − *j* →

> ⎞ ⎟ ⎠

→ *par* = *B* <sup>→</sup>(*B* <sup>→</sup>, *E* <sup>→</sup>)/*B*<sup>2</sup>

→ *per* = *B*

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

**3. Equations solutions: initial and boundary conditions**

<sup>→</sup> + *j* → *<sup>m</sup>* + *j* → <sup>2</sup> ) <sup>−</sup> \_\_\_\_\_\_\_\_\_ <sup>Ω</sup>*<sup>i</sup> mi* (Ω*<sup>i</sup>* <sup>2</sup> + *νin* 2 )]

+ \_\_\_\_\_ <sup>1</sup>

, (26d)

*me <sup>ν</sup>en*], (26e)

, (26f)

, (26g)

<sup>→</sup>]. (26h)

. (27b)

*<sup>s</sup>*) = 0, (27a)

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.
