© Author(s) 2008. This work is distributed under the Creative Commons Attribution 3.0 License. Advances in

Abstract. A mathematical model of the global neutral wind system of the Earth's atmosphere, developed earlier in the Polar Geophysical Institute (PGI), is utilized to simulate the large-scale global circulation of the middle atmosphere for January conditions. The utilized model enables to calculate not only the horizontal components but also the vertical component of the neutral wind velocity by means of a numerical solution of a generalized Navier-Stokes equation for compressible gas, so the hydrostatic equation is not applied. Global distributions of the horizontal and vertical wind, calculated for January conditions, are compared with simulation results, obtained earlier for conditions corresponding to summer in the northern hemisphere. It was found that the global distributions of the neutral wind, calculated both for winter and for summer periods in the northern hemisphere, in particular, the large-scale circumpolar vortices, are consistent with the planetary circulation of the atmosphere, obtained from observations.


Introduction
Not long ago, a new mathematical model of the global neutral wind system in the Earth's atmosphere has been developed in the PGI (Mingalev and Mingalev, 2005;Mingalev et al., 2007).This model produces three-dimensional global distributions of the zonal, meridional, and vertical components of the neutral wind and neutral gas density in the troposphere, stratosphere, mesosphere, and lower thermosphere.The characteristic feature of the model is that the vertical component of the neutral wind velocity, as well as horizontal components of the neutral wind, are obtained by means of a numerical solution of the appropriate momentum equation Correspondence to: V. S. Mingalev (mingalev@pgi.kolasc.net.ru) for a viscous gas without any simplification.The hydrostatic equation, usually utilized in the global circulation models, is not used in the new model.
This model has been applied in order to simulate the global circulation of the middle atmosphere for conditions corresponding to summer in the northern hemisphere (Mingalev and Mingalev, 2005;Mingalev et al., 2007).In the present study, this model is utilized to simulate the global distributions of the atmospheric parameters in the middle atmosphere for conditions corresponding to winter in the northern hemisphere.The simulation results are compared with results, obtained earlier, and with observational data.

Numerical model
The new model, developed in the PGI and utilized in this study, differs from existing global circulation models of the atmosphere, on principle.Firstly, as pointed out previously, the model does not include the pressure coordinate equations of atmospheric dynamic meteorology, in particular, the hydrostatic equation.Instead, the vertical component of the neutral wind velocity is obtained by means of a numerical solution of the appropriate momentum equation, with whatever simplifications of this equation being absent.Secondly, the model does not include the internal energy equation for the neutral gas.Instead, the global temperature field is assumed to be a given distribution.It is known that the atmospheric temperature distributions, calculated by using the existing global circulation models containing the internal energy equation, as a rule, differ from observed distributions of the atmospheric temperature.These differences are conditioned by uncertainty and complexity in various chemicalradiational heating and cooling rates.Therefore, there is no reason to expect an exact correspondence between the calculated and measured temperatures of the neutral gas.On the other hand, over the last years empirical models of the Published by Copernicus Publications on behalf of the European Geosciences Union.
global atmospheric temperature field have been successfully developed which have a satisfactory accuracy.In the present study, we take the global temperature distribution from the NRLMSISE-00 empirical model (Picone et al., 2002) and consider it to be an input parameter of the utilized model.
The mathematical model, utilized in the present study, is based on the numerical solution of the system of equations containing the dynamical equation and continuity equation for the neutral gas.For solving the system of equations, the finite-difference method is applied.The dynamical equation for the neutral gas in vectorial form can be written as where ρ is the neutral gas density, V is the neutral wind velocity, F is the acceleration comprising the gravity acceleration, Coriolis acceleration, acceleration of translation, and acceleration due to elastic collisions with the ion gas, and P is the total stress tensor.The latter tensor can be decomposed as follows: where p is the pressure, I is the unit tensor, and τ is the extra stress tensor whose components are given by the rheological equation of state or the law of viscous friction.A spherical coordinate system rotatable together with the Earth is utilized in model calculations.Therefore, from the dynamical equation, Eq. (1), momentum equations for the zonal, meridional, and vertical components of the neutral gas velocity may be derived.These equations include not only the pressure gradients but also partial derivatives of components of the extra stress tensor, τ .The latter tensor is composed of a Newtonian part, τ 0 , and a complementary part, τ 1 , namely, The former tensor, τ 0 , is given by the well-known Newton's law of viscous friction, where µ is the coefficient of molecular viscosity, whose dependence on the temperature is assumed to obey the Sutherland's law, and εis the tensor defined as where D 0 is the strain rate tensor and T r( ) denotes the trace of a tensor.The complementary stress tensor, τ 1 , is supposed to be conditioned by a small-scale turbulence having the scales equal and less than the steps of the finite-difference approximations.It is assumed that this tensor represents the effect of the turbulence on the mean flow and is given by an expression, analogous to the Newton's law of viscous friction, Eq. ( 4), with the scalar coefficient of viscosity, µ, being replaced by three distinct coefficients describing the eddy viscosities in the directions of the basis vectors of the utilized spherical coordinate system.For computing the eddy viscosities, the turbulence theory of Obukhov (1988) is applied.Thus, the momentum equations for the zonal, meridional, and vertical components of the neutral gas velocity acquire ultimately a form of a generalized Navier-Stokes equation for compressible gas on scales which are more than the steps of the finite-difference approximations, with the effect of the turbulence on the mean flow being taken into account by using an empirical subgrid-scale parameterization.The steps of the finite-difference approximations in the latitude and longitude directions are identical and equal to 1 degree.The height step is non-uniform and does not exceed the value of 1 km.
The simulation domain is the layer surrounding the Earth globally and stretching from the ground up to the altitude of 126 km at the equator.Upper boundary conditions provide the conservation law of mass in the simulation domain.The Earth's surface is supposed to coincide approximately with an oblate spheroid whose radius at the equator is more than that at the pole.More complete details of the utilized model may be found in the studies of Mingalev and Mingalev (2005) and Mingalev et al. (2007).

Simulation results and discussion
For the present study, the simulations were performed for winter period in the northern hemisphere (16 January) under moderate solar activity (F 10.7 =101) and low geomagnetic activity (Kp=1).The atmospheric temperature is one of major input parameters of the mathematical model, utilized in the present study.It is known that the atmospheric temperature varies in the course of time.A diurnal variation of the atmospheric temperature is described by the NRLMSISE-00 empirical model, applied in the present study.However, this variation is rather smooth, so the temperature may be approximately considered as a quasi-stationary function.To calculate distributions of atmospheric parameters, corresponding to a fixed moment, with the help of the utilized mathematical model, the method of establishment was used.Boundary conditions and inputs to the model were assumed to be time-independent and correspondent to a fixed moment of the universal time (UT).The stationary distributions of the atmospheric parameters were obtained for inputs to model corresponding to 10:30 UT.
Let us consider the distributions of the atmospheric parameters, calculated using the horizontally non-uniform temperature, obtained from the NRLMSISE-00 empirical model.The simulation results are partly shown in Figs.1-3.It is easy to see that the horizontal and vertical components of the wind velocity are changeable functions of latitude, longitude, and altitude.A maximal absolute value of the horizontal component of the wind velocity at a fixed height is larger at higher altitudes.The horizontal wind velocity can have various directions which may be opposite at the near points, displaced for a distance of a few steps of the finite-difference approximation.It is seen that horizontal domains exist where the steep gradients in the horizontal velocity field take place.
Incidentally, the presence of horizontal domains, where the steep gradients in the horizontal velocity field take place, is a novel feature of simulation results obtained with the help of the new mathematical model, developed in the PGI and utilized in the present study.The fact is that the meridional and zonal components of the wind velocity, calculated with the help of other mathematical models at levels of the middle atmosphere and lower thermosphere, possess usually the smooth variations with latitude and longitude.
From results of simulation, it can be seen that horizontal domains exist in which the vertical neutral wind component has opposite directions.The horizontal domains, where the vertical neutral wind component is upward, have a great length and large width, as a rule.Unlike, the horizontal domains, where the vertical neutral wind component is downward, have a large length and little width.Usually, the latter domains have a configuration like a long narrow band and coincide, as a rule, with the regions, where the steep gradients in the horizontal velocity field take place.Maximal absolute values of the upward vertical wind component are less than the maximal module of the downward vertical wind component.The peculiarities of configurations of horizontal domains, where the vertical neutral wind component is upward and downward, are novel features of simulation results, which were not produced by currently used general circulation models.
The vertical wind velocity can achieve values of a few m/s at levels of the lower thermosphere in the horizontal domains having a configuration like a limited narrow band.Comparing the vertical wind velocity distributions, calculated in the present study, with analogous results, obtained with the help of other existing general circulation models, one discovers a great discrepancy in absolute values of velocities.Commonly used general circulation models produce the vertical component of the wind velocity having the values of several centimeters per second at levels of the lower thermosphere.Such values are approximately a factor of 10 less than those obtained in the present study.Therefore, the presented results may appear as unusual and questionable.It is probable that the elucidation of the discussed discrepancy will require additional computational studies.However, the authors can express a supposition concerning this discrepancy.It is known that commonly used general circulation models are based on the application of the hydrostatic approach, which implies the utilization of a hydrostatic equation.Let us consider a problem of deriving a hydrostatic equation.This equation can be derived from the momentum equation for the vertical component of the neutral gas veloc-ity under specific conditions.The main condition consists in statics of air in the vertical direction.In other words, the vertical motion of atmospheric gas is assumed to be absent, whereas the neutral gas is not supposed to be at rest in the horizontal directions.Now, let us consider an inverse problem that consists in obtaining the neutral gas velocity with the help of the system of momentum equations derived using a hydrostatic approach.It is difficult to expect that the obtained flow of the neutral gas will conspicuously differ from a motion statical in the vertical direction.It is probable that the vertical component of the gas velocity will be close to that which was utilized for deriving a hydrostatic equation, that is, close to zero.It is the supposition of the authors that general circulation models, based on the utilization of the hydrostatic approach, cannot produce the vertical component of the gas velocity having visible values, on principle.However, the computational examination of the authors' supposition is beyond the scope of the present paper.
It is of interest to compare the simulation results with the experimental data.Unfortunately, the authors have failed to find in the literature observational data on global distributions of the vertical neutral wind component at levels of the mesosphere and lower thermosphere.Nevertheless, verticalwind measurements have been made in the mesosphere and lower thermosphere at some high-latitude stations.Peteherych et al. (1985) have observed vertical winds by optical Doppler measurements of the 557.7 nm emission in the auroral zone using Fabry-Perot spectrometer, located at Saskatoon, Canada (52 • N).Both upward and downward winds were observed, of 15 m/s magnitude, at altitudes of approximately 110 and 135 km.Price and Jacka (1991) found vertical winds of up to 30 m/s using Fabry-Perot spectrometer, located at Mawson station, Antarctica (67 • S), at altitudes less than 110 km.Ishii (2005) has observed vertical winds of up to 40 m/s using Fabry-Perot interferometer, installed at Poker Flat, Alaska (65 • N), at altitudes between 110 and 140 km in December 1999.Widdel (1987) has observed vertical velocities of air in the mesosphere using foil cloud experiments at Andenes, Norway (69 • N) during January and February 1984.He measured vertical winds of up to 10 m/s between 75 and 90 km.Measurements using the EISCAT UHF incoherent scatter radar show vertical winds in the range 6 m/s downwards to 2 m/s upwards in the height range 80-85 km over northern Scandinavia in November 1985and April 1986(Hoppe and Hansen, 1988).It can be seen that values of observed vertical winds may be much more than maximal values of the vertical wind component calculated in the present study.
It is know that the global atmospheric circulation can contain sometimes so called circumpolar vortices that are the largest scale inhomogeneities in the global neutral wind system.Their extent can be very large, sometimes reaching the latitudes close to the equator.It is well known from numerous observations that circumpolar vortices are formed at heights of the stratosphere and mesosphere in the periods close to summer and winter solstices, when there is no rebuilding of the atmosphere.The circumpolar anticyclone arises in the northern hemisphere under summer conditions, while the circumpolar cyclone arises in the southern hemisphere under winter conditions.On the contrary, the circumpolar cyclone arises in the northern hemisphere under winter conditions, while the circumpolar anticyclone arises in the southern hemisphere under summer conditions.Let us compare these experimental data with the simulation results obtained.
From the results shown in Figs.1-3, we can see that, for winter period in the northern hemisphere, at levels of the middle atmosphere, the motion of the neutral gas in the northern hemisphere is primarily eastward, so a circumpolar cyclone is formed.It can be noticed that the center of the northern cyclone is displaced from the pole for a distance corresponding to approximately 10 • of latitude.Simultaneously, the motion of the neutral gas is primarily westward in the southern hemisphere at levels of the middle atmosphere, so a circumpolar anticyclone is formed for summer period of the southern hemisphere.It can be seen that the circumpolar vortices of the northern and southern hemispheres, simulated in the present study at levels of the middle atmosphere, correspond qualitatively to the global circulation, obtained from observations.It is obvious that the calculated circulation of the atmosphere is completely conditioned by the horizontal non-uniformity of the temperature in the rotatable atmosphere.
It may be recalled that simulation results analogous to those, presented in Figs.1-3, have been obtained in the studies by Mingalev and Mingalev (2005) and Mingalev et al. (2007) for summer period in the northern hemisphere (10 June and 16 July, respectively).It turned out that simulation results obtained for 10 June and 16 July are rather similar, with primary directions of the motions of the neutral gas in the northern and southern hemispheres at levels of the middle atmosphere being opposite to those, calculated in the present study for January conditions.It has been found in the studies just mentioned that the global distributions of the neutral wind, calculated for summer period in the northern hemisphere, in particular, the large-scale circumpolar vortices, are consistent with the planetary circulation of the atmosphere, obtained from observations.
Incidentally, the qualitative agreement between the circumpolar vortices of the northern and southern hemispheres, simulated for winter and summer periods of the northern hemisphere, and the ones, obtained experimentally, manifests the adequacy of the distributions of the atmospheric temperature, calculated with the help of the NRLMSISE-00 empirical model (Picone et al., 2002) and utilized in our simulations.
The mathematical model of the global neutral wind system, utilized in the present study, has some intrinsic limitations.The internal energy equation for the neutral gas is not included in the model system of transport equations.Instead, the global temperature field is assumed to be a given distribution, i.e. the steady input parameter of the model, which is utilized for obtaining the steady-state distributions of the atmospheric parameters, using the method of establishment.Therefore, a transfer of energy is not considered in the model calculations.The steady state of the temperature is supposed to be conditioned by a balance between various chemicalradiational heating and cooling processes, which are complex and uncertain.The point of view is taken that this balance is retained during some period of time and provides the stationary state of the final global distributions of the zonal, meridional, and vertical components of the neutral gas velocity, density, and temperature.
Moreover, as a consequence of a stationary state of the atmospheric temperature, the mathematical model, utilized in the present study, cannot deal with a temporal behavior of velocity field, in particular, with a propagation of atmospheric waves.
In the mathematical model of the global neutral wind system utilized in the present study, the Earth's surface is assumed to coincide approximately with an oblate spheroid whose radius at the equator is more than that at the pole.It is of interest to study how the non-sphericity of the Earth influences on the formation of the planetary circulation of the lower and middle atmosphere.But a computational study of this problem is beyond the scope of the present paper and will be considered in future.

Conclusions
The mathematical model of the global neutral wind system of the Earth's atmosphere, developed earlier in the PGI, has been applied for simulation of the three-dimensional global distributions of the neutral wind components and neutral gas density for January conditions.A specific feature of the applied model is that the vertical component of the neutral wind velocity, as well as horizontal components of the neutral wind, are obtained by means of a numerical solution of the appropriate momentum equation for a viscous gas without any simplifications.Moreover, the internal energy equation for the neutral gas is not included in the system of transport equations of the model.Instead, the global temperature field is supposed to be a given distribution, i.e. the input parameter of the model, which is calculated with the help of the NRLMSISE-00 empirical model.
The results of simulation indicated that the horizontal domains exist in which the vertical neutral wind component is upward, with the lengths and widths of these domains being large.Simultaneously, the horizontal domains exist where the vertical neutral wind component is downward.The latter domains have configurations like a long narrow band and coincide, as a rule, with the regions in which the steep gradients in the horizontal velocity field take place.The vertical component of the wind velocity can achieve values of a few m/s at levels of the lower thermosphere in the limited horizontal domains.
The results of simulation indicated that, at levels of the stratosphere and mesosphere, the circumpolar cyclone is formed in the northern hemisphere and the circumpolar anticyclone is formed in the southern hemisphere.Simulation results of the present paper differ essentially from the results, obtained earlier for June and July conditions (Mingalev and Mingalev, 2005;Mingalev et al., 2007).The primary directions of the neutral wind in the northern and southern hemispheres at levels of the middle atmosphere, calculated in the present study for January conditions, are opposite to those, obtained earlier for June and July conditions.Both simulation results, calculated in the present study, and those, obtained earlier, are consistent with the planetary circulation of the middle atmosphere, obtained from observations.Incidentally, this fact manifests the adequacy of the applied mathematical model of the global neutral wind system.
It can be noticed that the utilized mathematical model of the global wind system of the Earth's atmosphere has the potential to be applied for numerical simulations of global circulations of atmospheres of other planetary bodies.A description of one such application of the utilized mathematical model to simulation of Titan's atmosphere circulation may be found in the study by Mingalev et al. (2006).

Fig. 1 .
Fig. 1.The distributions of the given neutral gas temperature (top panel), vector of the calculated horizontal component of the neutral wind velocity (middle panel), and calculated vertical component of the neutral wind velocity (bottom panel) as functions of longitude and latitude at the altitude of 50 km.The temperature is given in K and wind velocities are given in m/s, with positive direction of the vertical velocity being upward.

Fig. 2 .
Fig. 2. The distributions of the given neutral gas temperature (top panel), vector of the calculated horizontal component of the neutral wind velocity (middle panel), and calculated vertical component of the neutral wind velocity (bottom panel) as functions of longitude and latitude at the altitude of 70 km.The temperature is given in K and wind velocities are given in m/s, with positive direction of the vertical velocity being upward.

Fig. 3 .
Fig. 3.The distributions of the given neutral gas temperature (top panel), calculated zonal component of the neutral wind velocity (middle panel), and calculated meridional component of the neutral wind velocity (bottom panel) as functions of latitude and altitude at 23 • E longitude (noon meridian).The temperature is given in K and wind velocities are given in m/s, with positive directions being eastward and poleward, respectively.