- 4.11.1 Free atmosphere turbulent diffusivities
- 4.11.2 ``Non-local'' atmospheric boundary layer scheme
- 4.11.3 Discretization of the vertical diffusion equations
- 4.11.4 Solution of the vertical diffusion equations
- 4.11.5 Discrete equations for , , and

4.11 Vertical Diffusion and Boundary Layer Processes

In the near future, the gravity wave parameterization will also be called from within the vertical diffusion. This will allow the turbulent and, especially, the molecular diffusivity to be passed to the gravity wave parameterization to damp vertically propagating waves. The gravity wave parameterization may return additional diffusivities and tendencies to be applied before the actual diffusion is applied.

As in CCM2 and CCM3, the turbulence parameterization in CAM 3.0 includes computation of diffusivities for the free atmosphere, based on the gradient Richardson number, and an explicit, non-local Atmospheric Boundary Layer (ABL) parameterization. The ABL parameterization includes a determination of the boundary layer depth. In practice, the free atmosphere diffusivities are calculated first at all levels. The ABL scheme then determines the ABL depth and diffusivities and replaces the free atmosphere values for all levels within the ABL, returning both the updated diffusivities and the non-local transport terms. The implementation of the ABL parameterization in CCM2 is discussed in [74], while the formalism only is discussed here. Following the ABL scheme, molecular diffusivities are computed if the model top extends above 90 km (0.1 Pa).

As described in Boville and Bretherton [25], a general vertical diffusion
parameterization can be written in terms of the divergence of
diffusive fluxes:

where is the dry static energy, is the geopotential height above the local surface (does not include the surface elevation) and is the heating rate due to the dissipation of resolved kinetic energy in the diffusion process. The diffusive fluxes are defined as:

The viscosity and diffusivities are the sums of: turbulent components , which dominate below the mesopause; and molecular components , which dominate above 120 km. The turbulent diffusivities are the sum of two components, free atmosphere and boundary layer diffusivities, defined below. In the future, these terms also may include effective diffusivities from the gravity wave parameterization. The non-local transport terms are given by the ABL parameterization. Note that , as defined in (4.451) and implemented in CAM 3.0, does not include the term which causes diffusive separation of constituents of differing molecular weights. The molecular diffusion in CAM 3.0 is currently incomplete and should be used with caution. The molecular viscosity and diffusivities are all currently defined as . A more complete form, allowing separation of constituents, will be implemented later.

The kinetic energy dissipation term in (4.449) is
determined by forming the equation for total energy from
(4.448-4.449):

The diffusive kinetic energy flux in (4.454) is

(4.455) |

and the kinetic energy dissipation is

To show that is positive definite, we use (4.450) to expand for and :

We show that energy is conserved in the column by integrating (4.454) in the vertical, from the surface (), to the top of the model ( ):

(4.458) |

Therefore, the vertically integrated energy will only change because of the boundary fluxes of energy, of which only the surface heat flux, , is usually nonzero. It is typically assumed that the surface wind vanishes, even over oceans and sea ice, giving . Then, the surface stress does not change the total energy in the column, but does result in kinetic energy dissipation and heating near the surface (see below) For coupled models, nonzero surface velocities can be accommodated by including on both sides of the surface interface.

Here is the local shear, defined by

and the mixing length is generally given by

where is the Von Karman constant, and is the so-called asymptotic length scale, taken to be 30 m above the ABL. Since the lowest model level is always greater than 30 m in depth, is simply set to 30 m in CAM 3.0. Furthermore, denotes a functional dependence of on the gradient Richardson number:

where is the virtual potential temperature,

For simplicity, in the free atmosphere, we specify the same stability functions for all . For unstable conditions we choose

This means that no distinction is made between turbulent vertical diffusion of heat, scalars and momentum outside the boundary layer. However, separate coefficient arrays are maintained and other parameterizations (such as gravity wave drag) may provide distinct diffusivities. We also note the the turbulent diffusivity is the same for all constituents, even within the ABL. However, the molecular diffusivities differ for each constituent since they depend on it's molecular weight.

The free atmosphere turbulent diffusivities, described above, are an
example of the local diffusion approach. In such an approach, the
turbulent flux of a quantity is proportional to the local gradient of
that quantity (*e.g.* (4.450)-(4.451)). In addition,
the eddy diffusivity depends on local gradients of mean wind and mean
virtual temperature (see (4.459)). These are reasonable
assumptions when the length scale of the largest turbulent eddies is
smaller than the size of the domain over which the turbulence
extends. In the Atmospheric Boundary Layer (ABL) this is typically
true for neutral and stable conditions only. For unstable and
convective conditions, however, the largest transporting eddies may
have a size similar to the boundary layer height itself, and the flux
can be counter to the local gradient [76,48].
In such conditions a local diffusion approach is no longer
appropriate, and the eddy diffusivity is better represented with
turbulent properties characteristic of the ABL. We will refer to such
an approach as non-local diffusion.

To account for ``non-local'' transport by convective turbulence in the ABL, the local diffusion term for constituent is modified as in (4.451):

where is the non-local eddy diffusivity for the quantity of interest. The term is a ``non-local'' transport term and reflects non-local transport due to dry convection. Eq. (4.466) applies to static energy, water vapor, and passive scalars. No countergradient term is applied to the wind components, so (4.450) does not contain these terms. For stable and neutral conditions the non-local term is not relevant for any of the quantities. The eddy diffusivity formalism is, however, modified for all conditions.

In the non-local diffusion scheme the eddy diffusivity is given by

where is a turbulent velocity scale and is the boundary layer height. Equation (4.467) applies for heat, water vapor and passive scalars. The eddy diffusivity of momentum , is also defined as (4.467) but with replaced by another velocity scale . With proper formulation of (or ) and , it can be shown that equation (4.467) behaves well from very stable to very unstable conditions in horizontally homogeneous and quasi-stationary conditions. For unstable conditions and are proportional to the so-called convective velocity scale , while for neutral and stable conditions and are proportional to the friction velocity .

The major advantage of the present approach over the local eddy diffusivity approach is that large eddy transport in the ABL is accounted for and entrainment effects are treated implicitly. Above the ABL, so (4.466) reduces to a local form with given by (4.459). Near the top of the ABL we use the maximum of the values by (4.459) and (4.467), although (4.467) almost always gives the larger value in practice.

The non-local transport term in (4.466), , represents non-local influences on the mixing by turbulence [48]. As such, this term is small in stable conditions, and is therefore neglected under these conditions. For unstable conditions, however, most heat and moisture transport is achieved by turbulent eddies with sizes on the order of the depth of the ABL. In such cases, a formulation for consistent with the eddy formulation of (4.466) is given by

where is a constant and is the surface flux (in kinematic units) of the transported scalar. The form of (4.468) is similar to the one proposed in Holtslag and Moeng [76]. The non-local correction vanishes under neutral conditions, for which .

The formulations of the eddy-diffusivity and the non-local terms are dependent on the boundary layer height . The CCM2 configuration of this non-local scheme made use of a traditional approach to estimating the boundary layer depth by assuming a constant value for the bulk Richardson number across the boundary layer depth so that was iteratively determined using

where is a critical bulk Richardson number for the ABL; and are the horizontal velocity components at ; is the buoyancy parameter and is the virtual temperature at . The quantity is a measure of the surface air temperature, which under unstable conditions was given by

where is a constant, is the virtual heat flux at the surface, is a virtual temperature in the atmospheric surface layer (nominally 10 m), represents a temperature excess (a measure of the strength of convective thermals in the lower part of the ABL) and unstable conditions are determined by . The quantity was calculated from the temperature and moisture of the first model level and of the surface by applying the procedure in Geleyn [60]. The value of the critical bulk Richardson number in (4.469), which generally depends on the vertical resolution of the model, was chosen as = 0.5 for the CCM2.

Vogelezang and Holtslag [178] have recently studied the suitability of this
formulation in the context of field observations, large-eddy
simulations [128], and an
turbulence closure
model [53]. They propose a revised formulation which
combines shear production in the outer region of the boundary layer
with surface friction, where the Richardson number estimate is based
on the differences in wind and virtual temperature between the top of
the ABL and a lower height that is well outside the surface layer (*i.e.*
20 m - 80 m). In addition to providing more realistic estimates of
boundary layer depth, the revised formulation provides a smoother
transition between stable and neutral boundary layers. Consequently,
CAM 3.0 employs the Vogelezang and Holtslag [178] formulation for estimating the
atmospheric boundary layer height, which can be written as

The quantities , , and represent the horizontal wind components and virtual potential temperature just above the surface layer (nominally 0.1). In practice, the lowest model level values for these quantities are used to iteratively determine for all stability conditions, where the critical Richardson number, , is assumed to be 0.3. The disposable parameter has been experimentally determined to be equal to 100 (see Vogelezang and Holtslag [178]). The computation starts by calculating the bulk Richardson number between the level of and subsequent higher levels of the model. Once exceeds the critical value, the value of is derived by linear interpolation between the level with and the level below.

Using the calculated value for and the surface fluxes, we calculate the velocity scales, the eddy diffusivities with (4.467), and the countergradient terms with (4.468), for each of the transported constituents. Subsequently, the new profiles for , , , and are calculated using an implicit diffusion formulation.

The turbulent velocity scale of (4.467) depends primarily on the relative height ( is boundary layer height), and the stability within the ABL. Here stability is defined with respect to the surface virtual heat flux . Secondly, the velocity scales are also generally dependent on the specific quantity of interest. We will assume that the velocity scales for mixing of passive scalars and specific humidity are equal to the one for heat, denoted by . For the wind components, the velocity scale is different and denoted by . The specification of and is given in detail by Troen and Mahrt [174]. Holtslag et al. [75] have rewritten the velocity scale, in terms of the more widely accepted profile functions of Dyer [54], and have given a new formulation for very stable conditions. Below we follow the latter approach.

For stable ( ) and neutral surface conditions ( ), the velocity scale for scalar transport is

where is the friction velocity defined by

Furthermore, is the dimensionless vertical temperature gradient given by Dyer [54],

for . Here is the Obukhov length, defined by

For ,

which matches (4.474) for . Equation (4.476) is a simple means to prevent from becoming too large (and too small) in very stable conditions. In stable conditions, the exchange coefficients for heat and momentum are often found to be similar. Therefore we may use =.

For unstable conditions , we have that and differ in the surface layer and in the outer layer of the ABL . For the surface layer, is given by (4.472) with

Similarly, is written as

where is the dimensionless wind gradient given by

In the surface layer, the scalar flux is normally given by

Comparison with (4.466) and (4.467) shows that, in the surface layer, we should have in (4.468) for consistency.

For the outer layer, and are given by

is the convective velocity scale. Furthermore, is the turbulent Prandtl number and is a constant. The latter is obtained by evaluating the dimensionless vertical wind gradient by (4.479) at the top of the surface layer, as discussed by Troen and Mahrt [174]. This results in . For very unstable conditions ( or , it can be shown with (4.481) that is proportional to 0.85 , while for the neutral case . The turbulent Prandtl number of (4.481) is evaluated from

for . Equation (4.484) arises from matching (4.466), (4.467), (4.468), and (4.480) at the top of the surface layer. As in Troen and Mahrt we assume that is independent of height in the unstable outer layer. Its value decreases from for the neutral case ( and ), to for in very unstable conditions.

In very unstable conditions, the countergradient term of (4.468) approaches

where , because for very unstable conditions we obtain . Since typically Troen and Mahrt [174], we have . Similarly, the temperature excess of (4.470) reads in this limit as . This leads to (= 0.85 ) = 8.5 in (4.470).

Finally, using the velocity scales described above, the flux equation (4.466) is continuous in relative height () and in the boundary layer stability parameter ( or ).

and discretized in a time-split form using an Euler backward time step. Before describing the numerical solution of the diffusion equations, we define a compact notation for the discrete equations. For an arbitrary variable , let a subscript denote a discrete time level, with current step and next step . The model has layers in the vertical, with indexes running from top to bottom. Let denote a layer midpoint quantity and let denote the value on the upper interface of layer while denotes the value on the lower interface. The relevant quantities, used below, are then:

(4.487) | |||

Like the continuous equations, the discrete equations are required to
conserve momentum, total energy and constituents. The discrete forms
of (4.448-4.449) are:

For interior interfaces, ,

Surface fluxes are provided explicitly at time by separate surface models for land, ocean, and sea ice while the top boundary fluxes are usually . The turbulent diffusion coefficients and non-local transport terms are calculated for time by the turbulence model described above, which is identical to CCM3. The molecular diffusion coefficients, described earlier, are only included if the model top is above km, in which case nonzero top boundary fluxes may be included for heat and some constituents.

The free atmosphere turbulent diffusivities , given by (4.459-4.465), are discretized as

The stability function is:

The neutral is calculated by

with m. The Richardson number in the free atmosphere is calculated from

Similarly to the continuous form (4.456), is
determined by separating the kinetic energy change over a time step
into the kinetic energy flux divergence and the kinetic energy
dissipation. The discrete system is required to conserve energy
exactly:

where we have assumed zero boundary fluxes for kinetic energy. This leads to

According to (4.498), the internal dissipation of kinetic energy in each layer is the average of the dissipation on the bounding interfaces , given by (4.499) and (4.500). Expanding (4.499) using (4.490) and recalling that ,

for and similarly for . The discrete kinetic energy dissipation is not positive definite, because the last term in (4.501) is the product of the vertical difference of momentum at two time levels. Although will almost always be , values may occur occasionally. The kinetic energy dissipation at the surface is

(4.502) |

Since the surface stress is opposed to the bottom level wind, the surface layer is heated by the frictional dissipation. However, is not guaranteed to be positive, since it involves the bottom level wind at two time levels.

Note that it has been assumed that the pressure does not change within the vertical diffusion, even though there are boundary fluxes of constituents, including water. This assumption has been made in all versions of the CCM and is still made in CAM 3.0. This assumption will be removed in a future version of CAM 3.0, since the implied horizontal fluxes of dry air, to compensate for the boundary flux of water, cause implied fluxes of other constituents.

A series of time-split operators is actually defined by (4.488-4.491) and (4.498-4.500). Once the diffusivities () and the non-local transport terms ( ) have been determined, the solution of (4.488-4.491), proceeds in several steps.

- update the and profiles using
- update the bottom level values of , , and using the surface fluxes;
- invert (4.488) and (4.490) for ;
- compute and use to update the profile;
- invert (4.488,4.489) and (4.491) for and ;

Note that since all parameterizations in CAM 3.0 return tendencies rather than modified profiles, the actual quantities returned by the vertical diffusion are .

The non-local transport terms, , given by (4.468), cannot be treated implicitly because they depend on the surface flux, the boundary layer depth and the velocity scale, but not explicitly on the profile of the transported quantity. Therefore, application of is not guaranteed to give a positive value for and negative values may not be removed by the subsequent implicit diffusion step. This problem is not strictly numerical; it arises under highly non-stationary conditions for which the ABL formulation is not strictly applicable. In practice, we evaluate

(4.503) |

and check the profile for negative values (actually for , where may be ). If any negative values are found, we set for that constituent profile (but not for other constituents at the same point).

Equations (4.488-4.491) constitute a set of four tridiagonal systems of the form

where indicates , , , or after updating from time values with the nonlocal and boundary fluxes. The super-diagonal (), diagonal () and sub-diagonal () elements of (4.504) are:

(4.505) | |||

(4.506) | |||

(4.507) |

The solution of (4.504) has the form

or,

Substituting (4.509) into (4.504),

Comparing (4.508) and (4.510), we find

The terms and can be determined upward from , using the boundary conditions

(4.513) |

Finally, (4.510) can be solved downward for , using the boundary condition

(4.514) |

CCM1-3 used the same solution method, but with the order of the solution reversed, which merely requires writing (4.509) for instead of . The order used here is particularly convenient because the turbulent diffusivities for heat and all constituents are the same but their molecular diffusivities are not. Since the terms in (4.511-4.512) are determined from the bottom upward, it is only necessary to recalculate , , and for each constituent within the region where molecular diffusion is important. Note that including the diffusive separation term for constituents (which will be in the next version of CAM 3.0) adds additional terms to the definitions of , , and , but does not otherwise change the solution method.

The dry static energy at step and level is

which can be calculated from by integrating the hydrostatic equation using the perfect gas law.

where is the geopotential, is the geopotential at the Earth's surface and is the surface pressure. A fairly arbitrary discretization of (4.516) can be represented using a triangular hydrostatic matrix ,

Note that (4.517) is often written in terms of the virtual temperature . The apparent gas constant includes the effect of water vapor and is defined as

where is the apparent gas constant for dry air and is the gas constant for water vapor.

Using (4.517) in (4.515), we have

The interface geopotential in (4.520) is defined as

and is evaluated from (4.518), using . Although the correct boundary condition on (4.521) is , within the parameterization suite it is usually sufficient to take .

The definition of the hydrostatic matrix depends on the numerical method used in the dynamics and is subject to constraints from energy and mass conservation. The definitions of for the three dynamical methods used in CAM 3.0 are given in the dynamics descriptions.

After is modified by diabatic heating in a time split process,
the new
can be converted into and
using (4.520):

with evaluated from using . Once is defined, (4.521) and (4.523) can be solved for and from the bottom up. Since the latter must normally recalculated if is modified, calculating and from involves the same amount of computation as calculating and from .