Computational study of temperature and density perturbations on atmospheric dynamics

This study examines the perturbation effect of temperature and density of moist air on atmospheric variables at 9°1’48″N, 38°44’24″E and 6.324 km above the Earth’s surface. The atmosphere is a compressible neutral moist air flowing on a rotating Earth as a model and it’s basic atmospheric parameters such as gas constant, transport coefficients, mixing ratio and specific heat capacities are considered to be temperature dependent and the Earth’s gravity changes with latitude and altitude. To describe the dynamics, we carried out a numerical computation using finite difference method on an unstaggered grid. Our results revealed that the response of all the variables have a plane wave pattern, in which specific heat capacities (SHCs), resultant wind speed (RWS) and water vapor mixing ratio (MR) increase with time at each latitude but vertical wind speed (VWS), specific enthalpy (SE) and pressure decrease with time at each latitude. The increase of MR with time is the result of water vapor flux into the air parcel and the decrease of SE and increase of RWS with time is the result of thermal energy to mechanical energy transformation. The decrease of VWS with time is the effect of the viscous force due to temperature dependence of viscosity.


Introduction
The atmosphere is a complex fluid system that produces chaotic (unpredictable) motions and this complexity result from numerous interactions between different physical processes operating at various locations [1][2][3][4]. Several papers have been published that examine the possible source of atmospheric perturbations on the Earth's atmosphere [5,6]. The restoring forces on the atmosphere are commonly identified as changes in air density, pressure and temperature fields of the air and are induced by compressibility, gravity and rotation of the globe [7][8][9][10].
Atmospheric waves provide efficient communication between various areas of the atmosphere [11,12]. Because of the wide range of applications, research of air waves employing 'perturbation techniques' has been conducted for a long time.The time-dependent disturbances in a uniform basic current caused by differential heating on a flat rotating earth surface was investigated by Stern [13] in 1954. The results showed that the mean motion of an atmospheric fluid would be specified in terms of an equivalent mountain that depends on the Coriolis parameter as well as the surface temperature and eddy conductivity of the heated region.
Fifteen years later, the numerical integration of the shallow water equation on a rotating Earth was done by Akira [14] to indicate the influence of orography on large-scale atmospheric dynamics. Andreev [15] investigated the role of atmospheric non-equilibrium in the generation of wave perturbations caused by solar radiation flux. He suggested that such radiation interaction with the atmospheric gas can be responsible to contribute significant effect on spectrum of wave perturbations. Moreover, the theoretical studies of atmospheric perturbations related to inhomogeneity of gravity in the two-dimensional case were also studied by Ingel et al [16].
The horizontal equations of motion for an inviscid homogeneous fluid under the effect of pressure disturbances and waves have recently been used by [17] to investigate the nonlinear processes of solitary waves and cyclone genesis pushed by a moving pressure disturbance in the atmosphere. The finding demonstrated that the nonlinear evolution equation for wave amplitude meets the Korteweg-de Vries equation with a forcing term.
For extended spatial and temporal propagation, the research outlined above takes into account the constant transport coefficients, specific heat capacities, gas constant and gravity of the Earth. These properties of moist air, on the other hand, vary noticeably when temperature and density differences exist in the surface atmosphere even for a short time [18][19][20][21][22].
The purpose of this work is to explain how the atmosphere initially in equilibrium responds if it is subjected to local density and temperature perturbations, how do these locally generated disturbances propagate through the atmosphere spatially as a function of time and, how do the various atmospheric variables characterizing the atmosphere, like pressure, wind velocity components, enthalpy, internal energy, heat capacity and mixing ratio, respond and change in space with time. In particular, we investigate the effects of these perturbations on the dynamics of the tropospheric portion of the moist atmosphere over Ethiopia. We accomplished this by carrying out a rigorous numerical computation of the various conservation equations of the atmosphere so as to generate numerical values of atmospheric variables at each point in space as a function of time. The moist air has temperature dependent transport coefficients and mixing ratio dependent specific heat capacities and gas constant. Furthermore, the effects of altitude and latitude dependence of gravitational acceleration on the dynamics of the atmosphere are being studied.
The contents of this paper are organized as follows. Section 2, describes the model, including a brief explanation of the equations. A numerical method for solving the conservation equations is presented in section 3. Section 4 contains the outcomes of our research as well as a discussion of the findings. Section 5 ends with a conclusion.

Model of atmospheric system
The atmosphere is modeled as a compressible neutral ideal moist air moving in the rotating Earth's reference frame under the influence of the pressure gradient, gravity, Coriolis, centrifugal and viscous forces. The moist air is subjected to different scales of perturbation of density, pressure and temperature. We considered the middle troposphere located in the range of latitude 9°N to 15°N, longitude 35°E to 45°E and altitude 4 km to 7 km above the surface of the Earth. A perturbation of density and temperature is applied at a point of 9°1′48″N, 38°44′24″E and 6.324 km above the Earth's surface. The model describes unsaturated moist atmosphere in dry seasons and in particular night time anywhere in the glob where the air is less humid, that is, humidity is much less than one. This means cloud and rain droplet formation can be neglected and thus no internal heat sources due to phase transitions are considered. Since both solar and terrestrial infrared radiations are absorbed mostly by the lower planetary boundary layer as it contains most of the green house gases, the radiation heat transfer in the middle and upper troposphere can be neglected.
Troposphere of the atmosphere is thus considered as unsaturated moist air consisting of water vapor and dry air, which behave as ideal gases [23]. The moist atmospheres gas constant (R m ) and specific heat capacities at constant volume (c v ) and pressure (c p ) vary with the mixing ratio (ω) as shown below.
Here, R d , c vd and c pd are the gas constant and specific heat capacity at constant volume and pressure at sea level of dry air, respectively, whose numerical values are given in appendix A.1.
Since the specific enthalpy h m of moist-air is an additive function (Dalton's Law), it is equal to the weighted average of the partial specific enthalpies of the dry air (h d ) and water vapor (h v ), resulting in The specific enthalpies h d and h v can be expressed in terms of specific heat capacity and temperature [24] as follows: Here, h do and h vo are the specific enthalpies of dry air and water vapor at the reference temperature T o = 273.15K and according to [25] their value is h do = 0 and h vo = L v (T o ) = 2501 kJkg −1 is the latent heat of vaporization and T m is the temperature of moist air in kelvin that varies with space and time. In the atmospheric range of temperatures 283.15 K to 323.15 K, the specific heat capacities c pd and c pv of dry air and water vapor, respectively, have constant values, which are given in table A.1. Using these numerical values in equation (3) and plugging the resulting expressions into equation (2), we get the following expression of the specific enthalpy of the moist air. But, according to thermodynamics of moist atmosphere, change in specific enthalpy is given by: where c p is given by equation (1). Equating the enthalpies in equations (5) with (4), we get a quadratic equation for ω whose physically acceptable root is :

Linear perturbation
We use the concept of an air parcel to describe a small mass element in the atmosphere consisting of dry air and water vapor. The actual state variables of the moist atmosphere after it is subjected to density and temperature perturbations can be written as the sum of equilibrium and perturbed state quantities [23,28,29]: where, ρ mo , T mo and p mo are the equilibrium state mass density, absolute temperature and pressure of the moist air, respectively, before perturbation. r ¢ m , ¢ T m and ¢ p m are the magnitudes of the mass density, temperature and pressure perturbations of the moist air, respectively. As one can see, the equilibrium state variables only spatially but perturbed state variables vary in both space and time. The variables λ, fandz denote longitude, latitude and altitude space coordinates, respectively.
The expressions of the total mass density ρ m and pressure p m of the moist air are determined by Dalton's law [29]: where ρ d and p d are the partial mass density and pressure dry air, respectively, and ρ v and p v are the partial mass density and pressure of water vapor, respectively. Applying this concept for the moist air in equilibrium and using the ideal gas equation for each gas, we get: 3), respectively. p mo and ρ mo depend only on spatial variables and don't change with time. They show the equilibrium pressure and density profile of the atmosphere as a function of altitude, latitude and longitude. ω o denotes the mixing ratio of moist air at the equilibrium state evaluated at temperature T mo and ε = R d /R v . Moreover, T mo is obtained by solving the latitudinally varying tropospheric temperature lapse rate as described by [30]: Here, a, b, and c are constants and their values are given in table A.1. Using equation (9) and the ideal gas equations for dry air and water vapor, the total pressure of the moist air becomes: Hence, the perturbed pressure of the moist air is the difference of the total pressure in equation (12) and equilibrium pressure in equation (10): The atmosphere obeys various conservation laws. As a result, we shall develop the equations of conservation of mass, momentum and energy using a spherical coordinate system, taking into account the effect of the perturbation on the dynamics of moist air.

The continuity equation formoist air
Upon using equation (8), the continuity equation describing the conservation of the total mass of the moist air in the air parcel can be described in spherical coordinate system as follows:  v in spherical coordinate, R e is the radius of the Earth and z is the altitude. Moreover, u, v and w are the wind velocities in the zonal, meridional and radial directions of the Earth, respectively.

The Navier-Stokes equations for moist air
The motion of an air parcel subjected to pressure gradient, Coriolis, centrifugal, gravitational and viscous atmospheric forces is governed by Newtons second law [28,29], which is called Navier-Stokes equation in fluid dynamics. Using the velocity vector  v in the Earth's rotating reference frame, the momentum balance is given by Navier-Stokes equation: where  W is the angular velocity of the Earth,  r is the radius vector from the center of Earth and  g is the gravitational acceleration of the Earth that vary with both latitude (g f ) and altitude (g r ) in the atmosphere [31] and their corresponding expressions are given by: where G and J n are the universal gravitational and Jeffery's constants, respectively, and M is the mass of the Earth. The terms P 2 , P 3 and P 4 are Legendre polynomials of degrees 2, 3 and 4, respectively, that can be generated with the help of recurrence formulas in a gravity model [31].
The last term in equation (15) is the viscous force on the moist air in the air parcel due to its motion in the surrounding atmosphere. It is defined as the divergence of the viscous stress tensor of rank two (J) as: vis For Newtonian fluids [32], the viscous stress tensor (J) is expressed as where σ represents a unit dyadic which is the gradient of the position vector and the exponent T indicates the the transpose of the tensor. Due to temperature dependence of the coefficient of viscosity, as shown in equation (7), the viscous force can be expressed as: f o is the viscous force exerted on the air parcel when the dynamic viscosity is constant and  m f is the viscous force that arises due to temperature dependence of the dynamic viscosity, which most models neglect. Their corresponding expressions are given by: Plugging equation (20) into (19) and using the result in equation (15), we get the momentum or Navier-Stokes equation. Upon decomposing the momentum equation into the the eastward (u), northward (v) and upward (w) components of the velocity vector  v , we get a coupled set of differential equations for the velocity components u, v and w in spherical coordinates which are solved numerically in this work.

The first law of thermodynamics formoist air
As discussed in describing the model of our system in section (2), there is no radiative heat transfer and internal heat source due to phase change. The conservation of internal energy is given by the first law of thermodynamics [29].
where e is the local specific internal energy of the moist air which as a function of temperature and mixing ratio is given by: m v m The vector  q represents sensible heat flux due to thermal conduction and according to Fourier's law of heat conduction [29], it is given by where k T is the temperature dependent thermal conductivity given in equation (7). Thus, in equation (21), the first term on the right side represents the rate of change of internal energy due to heat flux, the second term on the right side represents the rate at which work is done on the moist air parcel by the pressure gradient force and the last term is the dissipation function that represents the rate at which work is done on the moist air parcel by the viscous force. This is the rate at which mechanical energy is converted to internal energy.
Since thermal conductivity depends on temperature, which varies in both space and time, after some mathematical procedure the heat term can be dissected into two terms: where Q o is the component of heat flux term when thermal conductivity k T is constant and Q k arises as a result of variation of thermal conductivity with temperature, which is neglected by most numerical models. The expressions of Q o and Q k in the spherical coordinate system are given by: The expression of the dissipation function (ä) is defined by a double dot product of the tensors J and   v .

( )   Î = v J•• 26
Using a tensor identity, equation (26) can be expressed as: Using equations (7) and (18) and performing the divergence operation on the viscous stress tensor taking into account the temperature dependence of μ T , the dissipation function can be split into two components.
where ä o denotes the dissipation function for a constant coefficient of viscosity and ä μ is the dissipation function that arises due to temperature dependence of the coefficient of viscosity, which is also neglected by most numerical computations. Their mathematical expressions are given by: f are given by equation (20). Upon substituting the expressions in equations (24), (25) and (29) into (21) and applying the material derivative, we get the expression of the time variation of the internal energy of the moist air parcel in spherical coordinate system. In the above energy equation, Q k and ä μ are the additional terms that arise as a result of temperature dependence of thermal conductivity and coefficient of viscosity, respectively, which is a new contribution of this work.

Numerical analysis
In this work, the perturbation is a local increase in temperature accompanied by a decrease in density, but the model can be applied to any type of local or global perturbation. Local heating of the atmosphere generates a heat wave that results in expansion of the gas and thus decrease in density which generates a local or internal gravity wave. Therefore, taking this into account, the moist air is subjected to a small perturbation of temperature ¢ = T T 0.25 m mo and density r r ¢ = -0.025 m mo at initial time t=0 at a point described in the model system, section 2.
We developed a FORTRAN90 program and carried out a rigorous numerical computation using the finite difference method (FDM) on an unstaggered grid for all atmospheric variables. This method is efficient, powerful, and flexible for calculating the discrete points in both space and time [33]. To fulfill the stability condition, the spatial grid size along latitude, longitude, and altitude are Δf = 0.05°, Δλ = 0.05°, and Δz = 20 meters, respectively. The temporal time step for this specific study were 3 minutes and we set the wind speeds to (u = 10, v = 10, w = 0.1) ms −1 in the zonal, meridional and radial directions as initial values.
At t=0, the temperature T mo is determined by solving equation (11), the pressure p mo and the mass density ρ mo are determined by substituting equations (A-2), (A-8) and (6) into (10). The numerical value of all the constant parameters used in this study were presented in table A.1.
Analysis and explanation of atmospheric dynamics requires numerical values of the various atmospheric variables at each point in space at each time. The purpose of the numerical computation is therefore to explain how the atmosphere initially in equilibrium responds if it is subjected to density and temperature local perturbations by generating numerical values of the state variables at each point in space as a function of time. To get these numerical values, we carried out a rigorous numerical computation of the coupled set of Navier-Stokes equation (15) for the velocity components u, v, w, the continuity equation (14) for the mass density ρ m and the first law of thermodynamics (30) for the internal energy e, from which one can get the temperature T m using equation (22) at each time step.

Results and discussion
In this section, we present the plots of the data generated by the numerical computation and the corresponding discussions. The temporal and spatial pattern of the = + HWS u v 2 2 , VWS = w and = + + RWS u v w 2 2 2 are shown in figures 1, 2 and 3, respectively. As one can see in figures 1 and 3, the HWS and RWS have a similar pattern so that they have a constant value at each latitude but both behaves like a plane wave that travels with time over a planner surface that is inclined from the time axis by a fixed angle. The increase of RWS with time implies that there is a mechanical energy gain at the expense of internal energy. At time t = 0, thermal energy (temperature perturbation) was injected into moist air and with time that excess thermal energy results in a collective air movement due to pressure work term in the internal energy conservation equation, that is, . This means, the temperature and thus internal and enthalpy decrease with time as the thermal energy injected into the moist air as temperature perturbation results in the expansion of the gas and thus work is done by the pressure on the environment at the expense of internal energy. The decrease of SE and atmospheric pressure of the moist air as a function of time at each latitude is the consequence of this temperature decrease.
From figure 2, we observe that VWS decreases as a function of time at each latitude and has the same value at each latitude at each time step. The decrease of VWS as a function of time may be associated with the dragging effect of the radial component of gravity, g r given in equation (16) and the viscous force that arises from temperature dependence of dynamic viscosity,  m f given in equation (20). This force depends on the gradient of temperature, which decreases with time at each latitude. Since the vertical temperature gradient is much larger than the horizontal one it has a negligible effect on the HWS. The observation that HWS, VWS and RWS have a constant value at each latitude at each time step implies that the atmosphere is moving collectively so that each point of the moist air moves with the same speed at each time.  The response of SE to the perturbations is determined by solving equation (4) numerically and the plot of the resulting data generated by the numerical computation is depicted in figure 4. As the SE decreases as a function of time the RWS increases as a function of time at each latitude. This indicates that the thermal and mechanical specific energies are undergoing transformation due to the inclusion of viscous dissipation (favors transformation of mechanical energy to thermal energy) and the pressure work (favors transformation of thermal energy to mechanical energy) terms in the internal energy conservation equation, which serves as a means of energy exchange between thermal and mechanical energies. The moist air parcel also loses energy by conduction, which has two terms The component of heat flux that arises due to temperature dependence of thermal conductivity depends on temperature gradient and is thus one of the causes of the rapid decay of specific enthalpy and internal energy.
The effect of the perturbations on MR ω is determined by solving equation (6) numerically at each time step and the plot of the resulting data generated by the numerical computation is depicted in figure 5.
The MR has same value at each latitude as a function of time but it behaves like a plane wave that travels with time over a planner surface that is inclined from the time axis by a fixed angle. When we analyze its frequency by Fourier transform, the highest peak of power occurs at a temporal frequency of 0.89956 hertz. The increase of MR with time shows that there is a water vapor flux into the air parcel that is driven by gradient of water vapor concentration, gradient of pressure and gradient of temperature, as the general sensible heat flux and mass flux are functions of the gradients of these parameters. However, the sum of all the diffusive mass fluxes is zero and the sum of all the heat fluxes is also zero and thus don't appear in the conservation equations of the total mass and total internal energy conservation equations.
The effect of the perturbations on SHC c v and c p are determined by solving equation (1) numerically at each time step and the plot of the resulting data generated by the numerical computation are depicted in figures 6 and 7.
From the figures we observe that both c v and c p have a wave pattern similar to MR so that they have the same value at each latitude but increase with time. This is because c v and c p are functions of ω and thus as ω increases with time, c v and c p also increase with time in a similar wave pattern.  Finally, the total pressure is calculated using equation (12) after the perturbations and the plot of the data generated by the numerical computation is shown in figure 8. From figure 8, it is noticed that the perturbation results in a pressure decrease as a function of time at each latitude and a pressure increase as a function of latitude at each time. This pressure decrease is the consequence of temperature decrease with time at each latitude due to transformation of thermal energy to mechanical energy. Moreover, as the pressure profile as a function of time and latitude depicted in the above figure shows, the pressure increases with latitude at each time step. This result   agrees well with the global pressure pattern of the Earth's atmosphere, where the equator (f = 0°) is a low pressure region and at f = 30°south or north is a high pressure region, that is, pressure increases with latitude.

Conclusions
In this paper, a finite difference method has been used for computing the effect of temperature and density perturbation on compressible neutral ideal moist air that moves on the Earth's rotating reference frame. The gas constant and specific heat capacities depend on water vapor mixing ratio and the transport coefficients are temperature dependent. Moreover, the Earth's gravity changes with latitude and altitude.
The applied perturbations materialize as waves of the various fields of atmospheric variables traveling through the moist atmosphere in various wave patterns. The waves that are associated with VWS and SE have similar wave pattern so that they have the same value at each latitude at each time step but decrease as a function of time at each latitude. The waves that are associated with atmospheric pressure have values that decrease as a function of time at each latitude but increase as a function of latitude at each time step. The traveling wave patterns associated with MR, SHC, HWS and RWS have similar wave patterns so that they have a constant value as a function of latitude at each time step but increase as a function of time at each latitude.
The increase of MR as a function of time at each latitude implies the inward water vapor flux, which results in increase in the SHCs as a function of time at each latitude. The observation that HWS, VWS and RWS have a constant value as a function of latitude at each time step implies that the atmosphere is moving collectively so that each point of the moist air moves with the same speed at each time. The increase of HWS and RWS and the decrease of SE as a function of time at each latitude implies that there is a strong coupling between wind speeds and SE and thus there is a transformation of energy between thermal and kinetic energies throughout the system that give rise to temperature changes. Moreover, the decrease of VWS as a function of time at each latitude is the effect of the viscous force that arises due to temperature dependence of dynamic viscosity.
The complex mode of motion of the atmosphere is the result of competition between the various forces acting on air parcel, that is, gravity, pressure gradient, Coriolis, centrifugal and viscous forces. In addition to these, the complex processes heat conduction, viscous dissipation and pressure work energy transfer mechanisms also mediate these coupled interactions and energy transfer processes. The presence of water vapor in the atmosphere also affects the dynamics of the atmosphere.
In conclusion, the main cause of the different patterns of waves in the atmospheric layer that we considered is linearity and dissipation processes. This effect leads to the formation of secondary wave sources in the atmosphere and affects climate change. Thus, this study will serve as a tool to bring the new scholar up to date on the current development in the numerical weather prediction model with considering the effect of perturbation at a regional scale.