2.1 Dynamics

2.1.1 Governing equations

The dynamics module solves a simplified Navier-Stokes equation set:

     

$${\frac{\partial \rho_0 \overline{U}_i}{\partial t}} \overbrace{- \displaystyle \frac{\partial \overline{U}_j \rho_0 \overline{U}_i}{\partial x_j}}^{\rm transport}_{} \overbrace{- \overline{\frac{U_j'' \rho_0 \partial U_i'' }{\partial x_j}}}^{\rm turbulence}_{} \overbrace{- \frac{\partial \overline{p'}}{\partial x_i}}^{\rm press.\ grad.}_{} \overbrace{- 2 \epsilon_{ijk} \Omega_j \rho_0 \overline{U}_k}^{\rm Coriolis\ force}_{} \overbrace{- \delta_{i3} \rho_0 \left(\frac{\theta'}{\theta_0} - \frac{C_v}{C_p} \frac{p'}{p_0} \right) g}^{\rm buoyancy}_{}$$ = (1)

$${\frac{\partial \rho_0 \overline{\theta}}{\partial t}} \nabla \overline{\bf U} \rho_{0}^{} \overline{\theta} \nabla \overline{{\bf U}'' \rho_0 \theta''} \overbrace{+ \rho_0 S_{\theta}}^{\rm sources}_{}$$ = (2)

$${\frac{\partial \rho_0 \overline{q}}{\partial t}} \nabla \overline{\bf U} \rho_{0}^{} \overline{q} \nabla \overline{{\bf U''} \rho_0 q''} \rho_{0}^{}$$ = (3)

$${\frac{\partial \overline{C_m}}{\partial t}} \nabla \overline{\bf U} \overline{C_m} \nabla \overline{{\bf U}'' C_m''}$$ = (4)

$$\nabla \left( \frac{\partial \rho_0 \overline{U}_i}{\partial t} \right)$$ = 0. (5)

and the ideal gas law:

 

$$\rho {\frac{p}{R_A \cdot T}}$$ (6)

where eq. 1-4 are the conservation of momentum, heat, water vapor, and chemical species, respectively. Eq. 5 is called the anelastic approximation and eq. 6 is the ideal gas law. The definitions of the symbols are listed in Table 1.

 

 

$$\rho \overline{w''\theta''}$$ A₀ the synoptic value of A . A₀ if a function of the height above sea level, Q_H

$$\overline{A} \downarrow$$ the average of entity A within the volume of a grid cell, down-welling IR-radiation,

$$\overline{A} \uparrow$$ A' up-welling IR-radiation,

$$\overline{A}$$ A'' R_A the ideal gas constant for air,

a_H2O,CO2 absorption coefficient due to water vapor and carbon dioxide. S_A source term of entity A ,

a_O3 absorption coefficient due to total ozone, T temperature of air,

C_m concentration of the chemical species m , T_T the temperature at the top of the modeling domain

$$\cdot$$ C_p T_f temperature of the big leaf,

D_R indirect solar irradiance (sect. 2.2.1), T_s temperature of the soil surface,

E turbulent kinetic energy, U_i component i of wind U,

$$\bf$$ E_f the energy balance on the big leaf, vectorial wind speed,

$$\bf$$ E_s the energy balance at the soil surface, vectorial wind speed at the cell faces,

f_zL() a function which calculates z/L in dependence of T_s,T_f, and z/L , u the normal path length (in vertical direction).

I direct irradiance, z/L the Monin-Obukhov stability parameter (Monin and Obukhov; 1954),

$$\epsilon$$ I₀ earth-sun distance corrected solar constant, turbulence dissipation rate,

$$\epsilon$$ K eddy diffusivity, the broadband emissivity, depending on path length and temperature,

$$\Omega_{j}^{}$$ L_v specific evaporisation energy of water, the components of the earth rotation vector,

$$\rho_{0}^{}$$ m optical path length (along the sun rays), middle air density at a certain height,

$$\sigma$$ p air pressure at ground, the Stefan-Boltzmann constant,

$$\overline{p} \sigma_{f}^{}$$ p' the fraction of the soil covered by vegetation,

$$\tau_{D}^{}$$ p_s standard air pressure ( = 10⁵Pa ), correction term to account for the influence of aerosols,

$$\tau_{R}^{}$$ q water vapor mixing ratio, mean optical depth of the atmosphere along the light ray.

$$\rho \overline{w''q''} \theta \theta$$ Q_E

Q_G soil heat flux.

A is used here as a general place holder for any physical quantity.

2.1.2 Grid structure and solving procedure

  

Figure 1: Metphomod uses a rectangular coordinate system. Most of the physical parameters are stored in the center of the according grid cells. Only turbulence parameters and fluxes are stored at the faces. Topography is considered, by having two categories of grid cells: normal grid cells, and underground grid cells. All fluxes at the boundaries of underground grid cells set to zero.

Metphomod uses a Cartesian grid with rectangular geometry. The grid structure and organization is explained in Figure 1. Topography is represented by introducing a below-ground grid cell category, and by the fact that the fluxes at the lateral boundaries of these cells are zero. Most state variables are stored in the center of the grid cells and flux variables at the cell boundaries. Many models that are able to consider complex topography are now available (Buty; 1988; Flassak; 1988; Grell et al.; 1993; Liu et al.; 1997; Schluenzen et al.; 1996; Thunis; 1995; Vogel et al.; 1987; Walko et al.; 1993). Most of them use terrain-following coordinates, by applying a transformed set of governing equations, that was derived by applying tensor analysis. This method has several advantages: It makes it possible to represent fine topographical features without the need for introducing many model layers. It, however, has the disadvantage of introducing complicated mathematical terms into the system, which are not easy to understand and to test within the program. The magnitude of numerical errors produced in this way is difficult to estimate. Most of these terrain following models have problems in complex terrain with sharp angles (steep mountain walls; deep valleys) and require the underlying topography to be smoothed.

To keep things simple, Cartesian coordinates were used in the Metphomod -model. The grid cells can be stretched in the vertical direction. (A similar approach was chosen by Bartzis et al.; 1993.)

The solving procedures for the governing equations were adapted from a method proposed by Rhie and Chow (1983). Both, a non-hydrostatic and a hydrostatic solver were implemented. The non-hydrostatic solver includes the following steps ( $\bf$u represents the wind speed vector at the cell faces; $\bf$u ^*, and $\bf$u ^** are temporary wind speed vectors occurring within the numerical scheme.):

1.

Interpolate wind speeds from the center of the cells to their faces (u).

2.

Apply the buoyancy term to the upper and lower boundaries of the grid cells:

$$\Delta \left(\frac{\theta'}{\theta_0} - \frac{C_v}{C_p} \frac{p'}{p_0} \right)g$$ (7)

wherein u ^*₃ is the vertical wind component at the vertical cell faces, C_v and C_p are the heat capacity of air at constant volume at constant pressure, respectively. The buoyancy has to be interpolated from the center to the faces.

3.

Using the anelastic assumption (eq. 5), a diagnostic expression for pressure can easily be derived:

$$\nabla^{2}_{} {\textstyle\frac{1}{\Delta t}} \nabla \cdot \rho_{0}^{} \bf$$ (8)

Metphomod solves this equation with the preconditioned conjugated gradient method (Press et al.; 1988).

4.

Apply the pressure gradient to the wind speed at the faces:

$$\bf \bf \Delta {\textstyle\frac{1}{\rho_0}} \nabla$$ (9)

The wind speeds at the faces now describe a divergence free wind field.

5.

Calculate the change of ( $\Delta$$\bf$u = $\bf$u ^** - $\bf$u ) and interpolate it to the center.

6.

Calculate the Coriolis forces.

7.

Calculate advection using $\bf$u ^** to calculate the fluxes (section 2.1.3).

8.

Calculate turbulence (section 2.1.4).

9.

Calculate the soil atmosphere fluxes (section 2.3).

10.

Calculate IR-radiation divergence (section 2.2.2).

11.

Calculate emissions (section 2.4).

12.

Calculate deposition (section 2.4).

13.

Calculate gas phase chemistry (section 2.5).

Steps 1-13 then are repeated iteratively.

As an alternative, the model can be run in hydrostatic mode. The solving procedure then is as follows:

1.

Identical to non-hydrostatic mode.

2.

Calculate the pressure field by vertically integrating the hydrostatic pressure equation:

$${\frac{\partial p'}{\partial z}} \rho_{0}^{} \left(\frac{\theta'}{\theta_0} - \frac{C_v}{C_p} \frac{p'}{p_0} \right)g$$ (10)

3.

Apply the pressure gradient to the wind speeds at the horizontal cell boundaries:

$$\Delta {\textstyle\frac{1}{\rho_0}} \nabla$$ (11)

4.

Calculate the vertical wind speed at the top and bottom cell boundaries, using the anelastic assumption (eq. 5).

Steps 5-13 are identical to the non-hydrostatic mode. Both methods were tested and validated. The tests presented in this publication were carried out with the non-hydrostatic method.

2.1.3 Advection

  The advection equations are solved by using the multidimensional form of Smolarkievicz's (1984) mpdata scheme. It is based upon the upwind formulation, where only the adjacent upwind grid cells influence the advection to any specific model cell. Numerical diffusion is reduced by applying an anti-diffusion correction step (Smolarkievicz; 1984). The mpdata-scheme is mass conservative, monotone and has low numerical diffusion.

2.1.4 Turbulence closure scheme

  Metphomod considers three turbulent regions to account for vertical mixing:

1.

The first few centimeters near the ground belong to the viscous sub-layer. Is it assumed that molecular diffusion is the most important process in this layer. The layer was parameterized using equations proposed by Zilitinkevitch (1970) and Deardorff (1974), which were integrated into the soil module (section 2.3).

2.

The first meters above the ground form the surface layer. This layer is still strongly influenced by the soil. The momentum fluxes within this layer are assumed to be constant. The traditional similarity theory by Monin and Obukhov (1954), with experimental parameterizations by Businger et al. (1971) were used. To model the surface layer, Metphomod uses a special surface-layer grid point near the surface (Fig. 1).

3.

The rest of the model domain consists of the atmospheric boundary layer and the free atmosphere. They spread over all normal model layers. Two turbulence parameterization schemes are available in Metphomod :

(a)

a non-local first-order turbulence closure scheme according to Stull (1988), the so called transilient turbulence theory closure scheme which allows non-local mixing throughout the model domain in one or more time steps. Because this approach is based on a Markovian matrix operation it is numerically stable at any integration interval. Despite its easy application, the TTT-closure scheme works well in buoyant situations and situations with complex layering of the atmosphere. Unfortunately there does not seem to exist a very realistic parameterization for the evaluation of the transilient matrix coefficients. The parameterization used in Metphomod proposed by Stull (1988) overestimated vertical mixing in our calculations.

(b)

A k- $\epsilon$ turbulence closure scheme (e.g. Apsley and Castro; 1997; Duynkerke; 1987; Stull; 1988). This scheme includes prognostic differential equations for the turbulent kinetic energy E , and the turbulence dissipation rate $\epsilon$

$${\frac{\partial E}{\partial t}} \nabla \cdot \overline{\bf U} {\frac{\partial}{\partial z}} {\frac{K}{\sigma_E}} {\frac{\partial E}{\partial z}} \epsilon$$ = (12) (13)

$${\frac{\partial \epsilon}{\partial t}} \nabla \cdot \overline{\bf U} \epsilon {\frac{\partial}{\partial z}} {\frac{K}{\sigma_\epsilon}} {\frac{\partial \epsilon}{\partial z}} {\frac{\epsilon}{E}} \epsilon \epsilon \epsilon$$ = (14)

where P , the production of turbulent kinetic energy is the sum of energy production by wind shear S and buoyancy B :

P = S + B (15)

$$\left( \frac{\partial U_i}{\partial x_i} + \frac{\partial U_j}{\partial x_i} \right)\frac \partial \partial$$ S = (16)

$${\frac{K g}{\rho_0 \sigma_\theta}} {\frac{\partial U_i}{\partial z}}$$ B = (17)

The eddy diffusivity ( K ) is calculated as

 

$${\frac{c_\mu E^2}{\epsilon}}$$ (18)

These K parameters serve to solve the normal diffusion equation:

 

$$\overline{w'' A''} {\frac{\partial \overline A}{\partial z}}$$ (19)

where $\overline{w'' A''}$ is the turbulent vertical flux of entity A . The parameters for the equation set were taken from Apsley and Castro (1997):

c_$\mu$ = 0.09, C_$\epsilon$1 = 1.44, C_$\epsilon$2 = 1.92, $$\sigma_{E}^{}$$ = 1.0, $$\sigma_{\epsilon}^{}$$ = 1.11, $$\sigma_{\theta}^{}$$ = 0.9

The physical meaning and derivation of these parameters has been extensively discussed in literature (e.g. Apsley and Castro; 1997; Duynkerke; 1987). A fully implicit numerical solver was used to integrate equation (19). The k- $\epsilon$ turbulence is more complicated to apply and more time consuming than the TTT scheme. The results obtained from this parameterization, however, were found to give more realistic turbulent mixing, and therefore were used in the simulations presented here.

We see it as an advantage that these turbulence closure schemes do not rely on mixing height or any similar measure, because mixing height in complex terrain is often difficult to determine. For an evaluation of different turbulence schemes see Hurley (1997).

2.1.5 Boundary conditions

Any model which is working on domains smaller than global, have lateral boundaries which are not physical, and should, however, be treated as cleverly as possible. In most cases, values for all entities will be set by the model user. Metphomod has different numerical methods to deal with these boundary values, which can be used depending on the situation to simulate:

1.

The wind speed is set to a preset value. Temperature, humidity and concentrations of chemical species depend on the wind direction: If air enters the model, the border is fixed. At the downwind boundary of the model, a zero gradient outflow condition is used. To prevent from numerical distortion produced by the advection scheme, the advection anti-diffusion step is not applied at the downwind boundary (Clappier; 1998).

2.

Radiative boundaries are used. All Variables at the border are damped towards a user preset value:

$${\frac{\partial \overline A}{\partial t}} \underbrace{\frac{r}{\Delta t}(\overline A - A_0)}_{radiation\ term}^{}$$ (20)

where $\overline{A}$ is some entity, r the relaxation rate, $\Delta$t the integration time step and A₀ the preset value, respectively. This method is very stable. It, however, has the disadvantage of being expensive in computing time and memory consumption, because it needs several grid points near the border.

3.

Some additional options, such as cyclic or symmetric boundary conditions can be used to investigate theoretical situations.