Abstract

Conservation of energy and the incorporation of parameterized heating in an atmospheric model are discussed. Energy conservation is used to unify the treatment of heating and kinetic energy dissipation within the Community Atmosphere Model, version 2 (CAM2). Dry static energy is predicted within the individual physical parameterizations and updated following each parameterization. Hydrostatic balance leads to an efficient method for determining the temperature and geopotential from the updated dry static energy. A consistent formulation for the heating due to kinetic energy dissipation associated with the vertical diffusion of momentum is also derived. Both continuous and discrete forms are presented. Tests of the new formulation verify that the impact on the simulated climate is very small.

1. Introduction

The solution of the equations of motion in an atmospheric model requires the parameterization of processes that occur on scales smaller than those explicitly represented by the model. Parameterizations are typically included for several processes (e.g., radiation, convection, vertical diffusion), which are treated as physically distinct, even though they are intimately coupled. Because of the difference between horizontal and vertical length scales, the parameterized processes typically operate entirely in the vertical direction, treating the atmosphere as a set of distinct columns, which can be treated independently. This paper discusses conservation of energy in the Community Atmosphere Model, version 2 (CAM2; more information available online at http://www.ccsm.ucar.edu/models/atm-cam/), particularly within the parameterization suite, and a unification of the application of parameterized thermodynamic processes in terms of dry static energy.

CAM2 is the successor to the National Center for Atmospheric Research (NCAR) Community Climate Model, Version 3 (CCM3; Kiehl et al. 1998). In common with many atmospheric models, the column parameterizations in CAM2 are clearly separated from the solution of the resolved adiabatic dynamics. In fact, CAM2 supports three “dynamical cores”: an Eulerian spectral core, as in CCM0 through CCM3 (CCM0–3); a semi-Lagrangian spectral core based on Williamson and Olson (1994); and a finite volume core based on Lin and Rood (1997). Horizontal diffusion is the only parameterized process that acts across columns in CAM2 and it is treated within the dynamical cores. Conservation of energy within the resolved dynamics is the responsibility of the dynamical core and the latter two (which contain implicit diffusion) include ad hoc fixers to enforce conservation.

Energy conservation was not a significant consideration in developing the original CCM0, which contained an energy imbalance of ∼10 W m–2 (Williamson 1988), primarily due to inconsistencies in the vertical numerical approximations. Energy conservation was taken more seriously in CCM1–3, which conserved energy to ∼0.4 W m–2 (e.g., Boville and Gent 1998). This level of conservation was partly due to a cancellation of errors since individual processes conserved energy only to ∼1 W m–2. Energy conservation is now a serious issue because of the application of CAM2 in multicentury coupled ocean–atmosphere simulations in which even imbalances of 0.4 W m–2 can cause spurious long-term trends.

For historical reasons, the parameterizations of thermodynamic processes in CCM0–3 have involved evaluating tendencies for a mixture of temperature, potential temperature and dry static energy. Both computational constraints and the evolution of the conceptual framework for the individual parameterizations have contributed to this diversity. In addition, the forcing terms have been applied using a mixture of time and process splitting (see below). Recently, Williamson (2002) has reformulated the CCM3 parameterization suite entirely in time-split form. In this formalism, the parameterizations operate sequentially, on profiles modified by the previous parameterization.

In CCM0–3, heating was considered equivalent to temperature tendency within the parameterization suite, so that energy was not conserved for the intermediate states which follow time split parameterizations. The total heating due to all parameterizations was correctly reported to the dynamics and applied in a nearly conservative fashion. This manuscript describes a method for applying heating to the dry static energy so that the energy is locally conserved at all stages within the parameterization suite, regardless of the splitting method. The methods are equivalent when height is the vertical coordinate but the distinction becomes important for pressure-based vertical coordinates, as used in most global models. The meaning of energy conservation is stated precisely, including the effects of variable gas mixtures, although CAM2 is still somewhat inconsistent in its treatment of variations in water vapor.

The vertical diffusion in CAM2 has been reformulated using dry static energy and both the kinetic energy flux and the heating due to kinetic energy dissipation are derived. These terms depend correctly on the stress tensor as discussed by Fiedler (2000) and Becker (2001). In CCM1–3, the kinetic energy fluxes were ignored.

We begin by reviewing the total energy equation in section 2 and its application to the physical parameterizations, together with the hydrostatic approximation in section 3. The splitting procedures used for parameterizations are summarized in section 4, the procedure for deriving temperature and geopotential from dry static energy is discussed in section 5, and vertical diffusion and dissipation of kinetic energy are discussed in section 6. Energy conservation considerations for the dynamical cores are commented on briefly in section 7, the impact of the changes on the CAM2 simulation is discussed in section 8, and conclusions appear in section 9.

2. Total energy conservation

Energy conservation is a important property of the equations of motion, which should be maintained when approximations are made in order to solve them. Energy conservation can be guaranteed by forming an evolution equation for total energy and verifying that it is satisfied by the approximate equations. The total energy equation is obtained by adding the kinetic energy equation to the first law of thermodynamics. The derivations of these equations may be found in Gill (1982), among other sources, and will only be summarized here.

A fundamental statement of the first law for a fluid is dI/dt + dW/dt = Q, where I is the internal energy per unit mass, W is the work done by the fluid through expansion, Q is the heating rate, and d/dt represents a substantive derivative. For a variable mixture of ideal gases, each with a fixed specific heat at constant volume cυ, = RT, and IcυT, where R and cυ are the apparent (mass weighted) constants for the mixture, α is specific volume, p is pressure, and T is temperature. Noting that the work term is dW/dtp(/dt) = (d/dt)α(dp/dt), ωdp/dt and cpcυ + R, the first law is

 
formula

or, equivalently,

 
formula

Here Q may have contributions from the divergence of radiative fluxes, molecular conduction, phase change and precipitation processes, chemical reactions, and viscous dissipation. In the lower atmosphere, molecular conduction is usually replaced by turbulent transport.

For a variable mixture of gases, cυ, R, and cp are determined by their mass-weighted average over the gases in the mixture. In CAM2, we make the common approximation that dry air is well mixed and that water vapor is the only additional substance present in sufficient abundance to affect the total pressure. Let superscripts ω and d denote water vapor and dry air, respectively. Then the apparent constants are

 
formula

where q is specific humidity. Water vapor is sufficiently abundant in the Tropics so that vertical variations in R may change the sign of the static stability, while horizontal variations may produce pressure gradients even with uniform temperatures. Variations in cp are often considered to be less important, since they affect the response to heating, but do not directly affect the density and pressure gradients. CAM2 uses cp = cdp, so that (2) becomes dT/dtRTω/pcdp = Q/cdp. This is the thermodynamic equation solved by both the spectral (Eulerian and semi-Lagrangian) cores, although cp is used in the second term due to an oversight. The finite volume core converts (2) to a virtual potential temperature equation. We retain a variable cp in the derivations below, but properly accounting for variations in cp in a model will probably require predicting an enthalpy (cpT) related variable. This will become increasingly important as the thermosphere is incorporated into CAM2, since diffusive separation of constituents results in large variations of cp and R above ∼120 km.

The kinetic energy equation is derived from the momentum equation,

 
formula

where τ is the stress tensor due to molecular or turbulent viscosity and other notation is standard. Defining kinetic energy KV · V/2, and taking the scalar product of V with (5) we obtain

 
formula

where the last two terms are the divergence of the diffusive flux (FK = V · τ) and the dissipation (D = –αV · τ) of kinetic energy, as in Fiedler (2000).

We obtain the total energy equation by adding (6) and (1), using mass continuity (∇;th·;thV = α–1/dt), and noting that wgdΦ/dt, where Φ ≡ gz is the geopotential:

 
formula

The heating Q0 = QD in (7) excludes kinetic energy dissipation, which merely converts one form of energy to another. Energy conservation requires that Q in (2) include all kinetic energy dissipation, whether it arises explicitly from viscous processes in (5), or implicitly from imperfections in the numerical approximations of the dynamical core. The local value of the implicit dissipation is usually difficult to determine, but the globally averaged value may be determined by integrating (7).

Converting (7) to fixed volume elements and assuming that pV and FK vanish on the top and bottom boundaries, the global integral of (7) is

 
formula

where dA is an area element and zt is the height of the model top. The globally integrated total energy may change only due to the globally integrated nondissipative heating, which can be expressed as the net boundary fluxes of sensible, latent, and radiative energy. Implicit dissipation due to the numerical approximations will result in an imbalance in (8), with the left-hand side usually being less than the right-hand side. In order to maintain energy balance, the globally averaged implicit dissipative heating, determined by this imbalance, can be added to (2) as either a globally uniform heating, or other ad hoc function. The implicit heating for the CAM2 dynamical cores is discussed in section 7.

The form of the total energy equation used to construct conservative parameterizations is a variation of (7). Applying the product rule for a divergence and using mass continuity, α∇·(pV) = d(αp)/dtαp/∂t, which is substituted into (7) to obtain

 
formula

where scpT + Φ is the dry static energy. The dry static energy includes the internal energy (cυT), the potential energy (Φ), and the work done (RT = ) to raise the volume of a fluid element from 0 to α at constant p.

3. Hydrostatic approximation

Up to this point, we have not made the hydrostatic approximation:

 
formula

However, CAM2 is a hydrostatic model and we require that the state be in hydrostatic balance in addition to conserving energy. Furthermore, each atmospheric column is treated independently within the physical parameterization suite and the parameterizations are time split. A sequence of intermediate states is constructed and we require both energy conservation and hydrostatic balance for each state. In fact, within the CAM2 parameterization suite, we require local conservation for both energy and constituents within each layer of every column.

Satisfying these requirements begins by forming a dry static energy equation from (2). Expanding the substantive derivative in the second term of (2) and applying (10),

 
formula

where Vh is the horizontal velocity. We emphasize that (11), which is adopted in CAM2, is an exact statement of the first law of thermodynamics in a hydrostatic system. A similar dry static energy equation, in which the heating Q0 does not include D, may be obtained from (9) by neglecting K, as discussed in Gill (1982). This approximation is not made in CAM2.

From the hydrostatic form of (5), we obtain the hydrostatic kinetic energy equation,

 
formula

and the total energy equation is still (9). Note that in the hydrostatic approximation, the kinetic energy is redefined to include only the horizontal wind, K = Vh · Vh/2. The pressure tendency and all advective and horizontal gradient terms are treated within the dynamical core, so that within the parameterizations (11) and (12) are just

 
formula

where the subscript c refers to the column parameterizations and FK is the vertical flux of kinetic energy due to parameterized stresses. We distinguish between Qc = Q0 + Dc and Q = Q0 + D, because kinetic energy dissipation within the dynamical core, whether implicit or explicit, is treated separately. Within the parameterizations, the total energy is defined as s + K and its tendency is given by (13) + (14). Total energy is locally conserved by using (13) to add all parameterized heating, including Dc. Hydrostatic balance is maintained by using (10) to obtain T and Φ from s.

In the CAM2 column parameterizations, only vertical diffusion and orographic gravity wave drag dissipate kinetic energy by acting on momentum. The dissipation was computed in CCM2–3 by neglecting FK in (14) so D = –∂K/∂t = –V · ∂V/∂t for each process. This is actually the correct heating for orographic (zero phase speed) gravity waves, which have no vertical energy flux. The correct vertical integral is still obtained for molecular and turbulent diffusion processes, despite neglecting FK. The correct forms of FK and D for the vertical diffusion parameterization in CAM2 are derived below.

Latent energy is conserved by conserving constituents individually and computing the energy associated with phase changes. As in previous versions of the CCM, the latent heat of fusion is neglected in CAM2, but not in the land and sea ice models that determine the surface fluxes. This inconsistency results in a globally and annually averaged conservation error of ∼0.2 W m–2, as mentioned in Boville and Gent (1998). Of course, the error is systematic and is much greater over ice surfaces. The latent heat of fusion has been implemented in CAM2 and will be discussed in a separate manuscript.

A statement of vertically integrated energy conservation within the physical parameterizations acting on each column can be obtained by adding the vertical integrals of (13) and (14) and writing the heating Qc as the net latent heating due to moist processes plus the convergence of an upward total sensible heat flux Fs due to radiation and vertical transport:

 
formula

where Lυ is the latent heat of vaporization, zs and zt are the surface and model top heights, Fs is the heat flux, and C is the vertically integrated net condensation rate. In CCM0–3, C was just the precipitation rate. CAM2 includes a prognostic cloud water that acts as a reservoir of liquid water, so that C is not balanced instantaneously by precipitation. We have assumed that the stress vanishes at zt and that V(zs) = 0, so FK(zs) = FK(zt) = 0. Even if the motion of the sea and ice surfaces is accounted for when computing the surface stress, one must still assume that FK(zs) = 0 unless the kinetic energy flux is also accounted for in the ocean and sea ice. For most applications, Fs(zt) is just the net radiative flux, while Fs(zs) also includes the sensible heat flux. In the thermosphere, emission of longwave radiation is very inefficient and absorption of solar radiation is largely balanced by downward molecular diffusion of heat. Therefore, in upper-atmosphere applications Fs(zt) may include a downward diffusive flux of heat due to solar heating above zt.

The hydrostatic approximation is used to pose the parameterizations in pressure coordinates in CAM2. Since the pressure tendency is treated within the dynamical core, the pressure is assumed to be invariant over the parameterization step. This procedure is not strictly correct in the presence of parameterized water fluxes, which should change the local pressure. Omitting this effect implies a compensating flux of dry air. The net surface flux of water (precipitation minus evaporation) should change the surface pressure, but instead implies a surface flux of dry air. For the finite volume core, the dry pressure in each layer is corrected after the parameterizations, while conserving the masses of all constituents. An energy and momentum conserving form has now been tested, but is not in CAM2. The Eulerian and semi-Lagrangian cores require a fixed relationship between the layer pressures and the surface pressure, so a simple correction to the dry mass does not work and only a global correction to the dry air mass is applied.

4. Process and time splitting

The set of parameterizations in a model may be logically conceived of as operating simultaneously (process splitting) or sequentially (time splitting) within a time step. In process splitting, each parameterization computes a heating based on the same input state,

 
formula

where Qi is the heating due to a particular parameterization. Although the same states enter each parameterization, the parameterizations are not independent. For example, radiative heating depends on cloud properties, which depend on convection and large-scale condensation. Therefore, the solution may depend on the order of the parameterizations. In time splitting, the order dependence of the parameterizations is made explicit. The partial time derivative in (13) is evaluated after each parameterization and the state is updated sequentially:

 
formula

After the application of all I parameterizations, the net heating is calculated as

 
Q = (sIs)/δt.
(18)

In CCM3, the first three parameterizations (penetrative convection, shallow convection, large-scale condensation), which involve condensation and latent heat release, were time split while the last three parameterizations (radiation, vertical diffusion, and gravity wave drag) were process split. Time splitting the moist processes allows precise enforcement of the constraint that relative humidity ≤100%. In CAM2, all parameterizations are time split. The surface turbulent fluxes are computed after radiation and are applied as boundary conditions in vertical diffusion.

Application of (17) is straightforward if z is the vertical coordinate, since Φ is invariant. However, CAM2 actually solves the parameterizations in pressure coordinates and the sequence defined by (17) requires that Ti and Φi be derived from si after each parameterization, using (10). The time-split parameterizations in CCM3 (and earlier versions) were not actually posed as in (17). Instead, heating was equated with temperature tendency, so that

 
Ti = Ti–1 + δtQi/cp,
(19)

and hydrostatic balance was imposed by integrating (10) to determine Φi. While this procedure is apparently more straightforward, it does not conserve energy within the parameterization suite.

Energy was conserved outside of the parameterizations by reporting the net heating from the time-split parameterization as

 
formula

with Qi = 0 for time-split parameterizations. Williamson (2002) time split all parameterizations in CCM3 using this method.

An efficient method to solve for Ti and Φi from si is given below, allowing (17) to be used at the same cost as (19). The difference is shown in Fig. 1 for a uniform heating of Q/cp = 1 K day–1 and for the net convective heating from CAM2 over the western equatorial Pacific Ocean. Applying a uniform heating in the vertical results in a temperature tendency that decreases with height. Heating increases both the temperature and the thickness of lower layers, therefore increasing Φ in layers above. In the absence of heating at upper levels, conservation of s would require ∂T/∂t < 0 as Φ increases. This can be clearly seen for the convective heating, which is >0 below 200 mb (except in the surface layer) and ∼0 above. The resulting temperature tendency is smaller than Q/cp at all levels and is negative above 200 mb.

Fig. 1.

Heating rate Q/cp (dash, K day–1) and temperature tendency ∂T/∂t, (solid), (left) for an idealized heating Q/cp = 1 and (right) for the actual net tropical convective heating from CAM2. Both solutions for ∂T/∂t use the average tropical profile of T and q from CAM2

Fig. 1.

Heating rate Q/cp (dash, K day–1) and temperature tendency ∂T/∂t, (solid), (left) for an idealized heating Q/cp = 1 and (right) for the actual net tropical convective heating from CAM2. Both solutions for ∂T/∂t use the average tropical profile of T and q from CAM2

We note that the time-splitting method does not allow all properties of the original equations to be preserved for intermediate states. We have chosen energy conservation and hydrostatic balance as the properties to preserve, while potential temperature conservation has not been preserved within the parameterization suite.

5. Discrete equations for s, T, and Φ

We define a compact notation for use in the discrete equations below. For an arbitrary variable ψ, let a subscript denote a discrete time level, with current step ψn and next step ψn+1. The model has L layers in the vertical, with indices running from top to bottom. Let ψk denote a layer midpoint quantity and let ψk denote the value on the upper interface of layer k while ψk+ denotes the value on the lower interface. The relevant quantities are then

 
formula

Using the above notation, the dry static energy at step n and level k is

 
formula

which can be calculated from Tn by integrating (10):

 
formula

where Φs is the geopotential at the earth’s surface and ps is the surface pressure. A fairly arbitrary discretization of (22) can be represented using a triangular hydrostatic matrix 𝗛kl:

 
formula

Note that (23) is often written in terms of the virtual temperature Tυ = TR/Rd. Using (23) in (21),

 
formula

The interface geopotential in (25) is defined as

 
formula

and cp and R are evaluated from (3) and (4), using qn.

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 (see, e.g., Williamson and Olson 1994). The definitions of 𝗛 for the three dynamical methods used in CAM2 are given in appendix A.

If sn is modified by diabatic heating in a time-split process, then the new sn+1 = sn + Qnδt can be converted into Tn+1 and Φn+1 using (25):

 
formula

with ckp and Rk evaluated from (3) and (4) using qkn+1 (recall that ckpcdp in CAM2). Once 𝗛 is defined, (26) and (27) can be solved for Tn+1 and Φn+1. Calculating T and Φ from s involves the same amount of computation as calculating Φ and s from T.

The only significant constraint imposed on the numerical algorithm to obtain (27) is that 𝗛 be triangular, which implies that Φ(p) depends only on T(p′≥ p). This constraint is a property of the continuous system represented by (22) and is obeyed by most numerical algorithms in use today. In fact, for “well-behaved” algorithms (including those in CAM2), 𝗛 is an upper-triangular matrix with the special property that all above-diagonal entries in each column are identical. That is, the thickness of a layer l in (23) is independent of the level k to which the sum is taken. Therefore, the solution of (27) can proceed from the bottom up and is much cheaper than the matrix notation might lead one to believe.

6. Vertical diffusion

A general vertical diffusion parameterization for a hydrostatic model can be written in terms of the divergence of diffusive fluxes:

 
formula

where D is the heating rate due to the dissipation of resolved kinetic energy in the diffusion process. In the notation of section 2, (Fu, Fυ) are the stress tensor components τ31, τ32 corresponding to horizontal stresses on interfaces normal to the vertical direction. These are the only stresses relevant for hydrostatic column processes. We note that the vertical diffusion of heat is written here in terms of s, which is the variable diffused in CAM2. The vertical diffusion in CCM3 operated on θ, which is very similar to operating on s, although it is easier to write formally conservative operators for s.

Following Holtslag and Boville (1993), the diffusive fluxes are defined as

 
formula

The viscosity Km and diffusivities Kq,s are the sums of turbulent components Ktm,q,s, which dominate below the mesopause; and molecular components Kmm,q,s, which dominate above ∼120 km. The nonlocal transport terms γq,s apply in the boundary layer and represent the effects of horizontally subgrid-scale eddies with vertical scale comparable to the boundary layer depth.

We form the equation for total energy, equivalent to (9), from (28)–(29):

 
formula

The diffusive kinetic energy flux in (34) is

 
formula

and the kinetic energy dissipation is

 
formula

These definitions correspond precisely to the more general expressions of FKE and D in terms of the stress tensor introduced after (5). To show that D is positive definite, we use (30) to expand for Fu and Fυ:

 
formula

In CAM2, (28)–(31) are converted to pressure coordinates using (10) and solved using a Euler backward time step. Like the continuous equations, the discrete equations are required to conserve momentum, total energy, and constituents. The discrete forms of (28)–(29) are

 
formula

For interior interfaces, 1 ≤ kL – 1,

 
formula

Surface fluxes FL+u,υ,q,s are provided explicitly at time n by separate surface models for land, ocean, and sea ice while the top boundary fluxes are usually F1–u,υ,q,s = 0. The turbulent diffusion coefficients Ktm,q,s and nonlocal transport terms γq,s are calculated for time n by the turbulence model, which is identical to CCM3. The molecular diffusion coefficients are only included if the model top is above ∼90 km, in which case nonzero top boundary fluxes may be included for heat and some constituents. The formulation of thermospheric processes in a vertically extended version of CAM2 will be described elsewhere.

Similarly to the continuous form (36), Dk 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:

 
formula

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

 
formula

According to (43), the internal dissipation of kinetic energy in each layer Dk is the average of of the dissipation on the bounding interfaces dk±u,υ , given by (44) and (45). Expanding (44) using (40) and recalling that un+ = (un+1 + un)/2,

 
formula

for 1 ≤ kL – 1 and similarly for dk+υ . The discrete kinetic energy dissipation is not positive definite, because the last term in (46) is the product of the vertical difference of momentum at two time levels. The surface layer is heated by the frictional dissipation associated with generating the surface stress, since the surface stress is opposed to the bottom-level wind. However, if  | FL+u |  or  | FL+υ |  is large enough to change the sign of uLn+ or υLn+ in (45), the dissipative heating may be negative.

A series of time-split operators is actually defined by (38)–(41) and (43–45). First, γq,s are used to update the input profiles. Although qn > 0, the updated profile for a constituent may occasionally contain negative values, in which case γq is discarded for that profile. As mentioned in Holtslag and Boville (1993), this problem usually arises under rapidly changing conditions for which the boundary layer formulation in the turbulence model is not strictly appropriate. Second, FL+u,υ,s,q and F1–u,υ,s,q are used to update the bottom and top model layers, respectively. Third, (38) is inverted for (u, υ)n+1, using the method given in appendix B. Fourth, the diffusive heating, Dk is added to sn before solving (39) for sn+1. Finally, (38) is solved for qn+1.

7. Energy conservation in the dynamical cores

The conservation properties of the numerical approximations in the dynamical cores are beyond the scope of this manuscript and we only comment on issues that are closely related to conservation within the physical parameterizations. Following Williamson (2002), the parameterization suite may be either time or process split from the dynamics, regardless of the method used within the parameterization suite. The Eulerian and semi-Lagrangian dynamical cores are both process split from the parameterizations, with heating and specific forces being applied within the dynamics. The finite volume (fv) core is entirely adiabatic and is time split from the parameterizations, taking updated T, u, υ as input. The current CAM2 implementation of the fv core takes Tn+1 = Tn + Q/cp as the input state, leading to a small energy imbalance. This issue will be addressed in a future model revision.

The dissipation of kinetic energy into heat D must be calculated explicitly and included in the heating Q in the first law of thermodynamics (2), in order to conserve energy. In Newtonian fluids, D is a positive definite quantity given by the product of the stress tensor and the velocity gradient, as discussed by Fiedler (2000) and Becker (2001). The spectral Eulerian core in CAM2 includes a biharmonic horizontal diffusion operator that cannot be represented by a symmetric stress tensor and therefore the kinetic energy dissipation cannot be correctly defined, as noted by Becker (2001). Instead, FK is ignored and D = –∂K/∂t = –V · (∂V/∂t)d, where (∂V/∂t)d is the specific force from the diffusion process. Note that D < 0 (cooling) if ∂( | V | /∂t)d > 0.

Both the finite volume and semi-Lagrangian dynamical cores in CAM2 have some diffusion implicit in their numerical approximations, for which D is unknown. In principle, D could be determined as a residual in (6), by predicting K in addition to V. In practice, (8) is used to define the global integral of D, which is then applied in (2) as a uniform heating. For all three dynamical cores, the globally averaged value of D is ∼2 W m–2 (D. Williamson 2002, personal communication).

Conservation of water vapor within the transport component of the dynamical core is also required in order to conserve the latent energy included in Q. The finite-volume core conserves constituents explicitly while the semi-Lagrangian transport algorithm, used by both the spectral Eulerian and semi-Lagrangian cores, contains an ad hoc mass fixer to ensure global conservation, as discussed in Rasch and Williamson (1990). Sources in the constituent equations are always time split in CAM2.

8. CAM2 results

The formalism discussed above is implemented in CAM2. The impact of the changes compared to CCM3 have been evaluated with a series of 5-yr simulations using the standard CAM2 configuration of spectral Eulerian dynamics with a triangular truncation at total wavenumber 42 (T42) and 26 layers in the vertical, extending from the surface to a top interface at ∼2 mb. (The tests actually used a prerelease version of CAM2 that had slightly different adjustable constants for some parameterizations.) All of the simulations used climatological sea surface temperatures. The comparison to “CCM3” has been made by returning specific processes to the CCM3 form while retaining the CAM2 form elsewhere. The CCM3 form requires updating T instead of s after parameterizations, computing D from ∂K/∂t instead of (37), mixing θ instead of s, and all of the above changes together. In addition to the time-splitting method and vertical diffusion parameterization discussed above, the primary differences between CAM2 and CCM3 are the inclusion of: the prognostic cloud water variable of Rasch and Kristjánsson (1998), the generalized cloud overlap treatment in the radiative transfer of Collins (2001), the updated treatment of the infrared radiative effect of water vapor of Collins et al. (2002), and a simple parameterization for evaporation of convective rainfall.

Changing between updating T and updating s had no significant impact on the simulation. The change does alter the conceptual framework for incorporating parameterizations and energy conservation within the parameterization suite can be verified. In fact, a few minor conservation errors in CAM2 were identified during this work that will be fixed in the next version. Most notably, the radiative heating rates (flux divergences) are computed hourly, saved, and applied to multiple time steps. The fluxes should be saved instead to allow for pressure changes. On an annual average basis, the standard configuration of CAM2 conserves energy to <0.3 W m–2, of which ∼0.2 W m–2 is due to the neglect of the latent heat of fusion within the atmosphere. The remainder is due to lack of conservation within the spectral dynamical core, including the use of cp from (3) instead of cdp.

For a model that does not extend at least into the mesosphere, such as the standard version of CAM2, almost all of the vertical diffusion takes place in the boundary layer. Figure 2 shows the annually averaged heating from vertical diffusion over the bottom 100 mb of the atmosphere. Note that 1 m W kg–1 corresponds to a temperature tendency of 0.0864 K day–1. The diffusion heats throughout the boundary layer on average, except in polar regions. Turbulent mixing of s will cool the upper part of the boundary layer, since s increases upward in a stably stratified atmosphere. However, the input of sensible heat at the surface dominates the turbulent cooling over most of the globe. The mixing in polar regions is primarily mechanically driven and the surface sensible heat flux is small, resulting in cooling in the upper part of the boundary layer.

Fig. 2.

Annual and zonal average of the heating from vertical diffusion (mW kg–1) for (left, contour interval 10) the CAM2 diffusion of s and (right, contour interval 1) the difference between CAM2 and CCM3 diffusion. Negative contours are dashed

Fig. 2.

Annual and zonal average of the heating from vertical diffusion (mW kg–1) for (left, contour interval 10) the CAM2 diffusion of s and (right, contour interval 1) the difference between CAM2 and CCM3 diffusion. Negative contours are dashed

The difference in vertical diffusion heating between the CAM2 and CCM3 forms is very small, as shown in Fig. 2, where the contour interval for the difference is 0.1 of that for the heating itself. The difference due to diffusing s rather than θ is negligible. Most of the change comes from using the correct definition (36) for kinetic energy dissipation. The dissipative heating for the two definitions is shown in Fig. 3. Kinetic energy dissipation is distributed through the boundary layer in the CCM3 form, while it is mostly confined in the surface layer in the CAM2 form. The turbulent mixing of momentum in the boundary layer results in a relatively smooth decrease of kinetic energy in the vertical. However, most of the kinetic energy change results from the kinetic energy flux in (34), rather than dissipation. The CCM3 form ignores this distinction. The vertical integral of the heating is identical for the two definitions, given the same input profile, and is indistinguishable in climate simulations.

Fig. 3.

Annual and zonal average of the heating from kinetic energy dissipation (mW kg–1) for (left) the CAM2 diffusion of s and (right) the CCM3 diffusion. The contour interval is 1 mW kg–1 and negative contours are dashed

Fig. 3.

Annual and zonal average of the heating from kinetic energy dissipation (mW kg–1) for (left) the CAM2 diffusion of s and (right) the CCM3 diffusion. The contour interval is 1 mW kg–1 and negative contours are dashed

The global and annual average of the explicit dissipative heating due to surface stress and vertical diffusion is 2.0 W m–2 in CAM2. The surface layer accounts for 1.65 W m–2 of the global mean, while the rest of the atmosphere accounts for only 0.35 W m–2. Although the contributions of dL+u,υ and dLu,υ were not separated in the CAM2 output, the dissipative heating in the surface layer results almost entirely from the surface stress terms. In fact, the heating is mostly due to surface stress in the oceanic storm tracks with somewhat more dissipation in the Southern Hemisphere than in the Northern Hemisphere (Fig. 3). There is a distinct seasonal cycle, maximizing in winter in each hemisphere, when the dissipation in the storm track exceeds 8 W m–2 (not shown).

9. Conclusions

The parameterized heating in CAM2 has been reformulated around the unifying concept of conservation of energy. Vertical diffusion and gravity wave drag, which alter kinetic energy, conserve total energy by including heating due to kinetic energy dissipation. This heating was included in an ad hoc form in previous versions (CCM1–3), rather than depending correctly on the stress tensor and velocity gradient.

Within the parameterization suite, heating is always applied to the dry static energy s in CAM2, rather than being applied to T, as in CCM3 and earlier versions. Therefore, CAM2 conserves energy at each stage within a time step, except for any nonconservation within the dynamical core. Earlier versions approximately conserved energy in aggregate, but not within the parameterization suite.

For time-split parameterizations, applying heating to s requires inverting the hydrostatic equation to determine T and Φ from the modified s. The inversion method developed for CAM2 is actually as efficient as evaluating the hydrostatic equation to determine Φ and s from T. We note that time splitting is not necessarily the best way of treating the physical parameterizations. As part of the unification of the formulation discussed here, the CAM2 code has been structured so that the parameterization suite can easily be converted to any combination of time and process splitting.

Incorporating the “new” formalism in CAM2 did not change the simulation significantly. The largest systematic change is that most of the heating from kinetic energy dissipation is applied in the surface layer, rather than being distributed through the boundary layer. The importance of the formalism lies in the precise definition of conservation of energy, which can be monitored at each stage within a time step and enforced within the numerical approximations.

The net energy conservation in CAM2 is very similar to that in CCM3 but does not arise from cancellation of errors. The remaining conservation errors in the parameterization suite arise from the omission of the latent heat of fusion and from the application of a fixed radiative heating rate over multiple time steps. In the next generation of CAM, the parameterizations will conserve precisely and the net conservation will depend only on the properties of the dynamical core.

Acknowledgments

The authors thank Prof. David Randall and two anonymous reviewers for their comments, which led to significant improvements on the original version of the manuscript. CFB also thanks NCAR for supporting his sabbatical, during which much of this work was completed.

REFERENCES

REFERENCES
Becker
,
E.
,
2001
:
Symmetric stress tensor formulation of horizontal momentum diffusion in global models of atmospheric circulation.
J. Atmos. Sci.
,
58
,
269
282
.
Boville
,
B. A.
, and
P. R.
Gent
,
1998
:
The NCAR Climate System Model, version one.
J. Climate
,
11
,
1115
1130
.
Collins
,
W. D.
,
2001
:
Parameterization of generalized cloud overlap for radiative calculations in general circulation models.
J. Atmos. Sci.
,
58
,
3224
3242
.
Collins
,
W. D.
,
J. K.
Hackney
, and
D. P.
Edwards
,
2002
:
An updated parameterization for infrared emission and absorption by water vapor in the National Center for Atmospheric Research Community Atmosphere Model.
J. Geophys. Res.
,
107
.
4664, doi:10.1029/2001J0001365
.
Fiedler
,
B. H.
,
2000
:
Dissipative heating in climate models.
Quart. J. Roy. Meteor. Soc.
,
126
,
925
939
.
Gill
,
A. E.
,
1982
:
Atmosphere–Ocean Dynamics.
Academic Press, 662 pp
.
Holtslag
,
A. A. M.
, and
B. A.
Boville
,
1993
:
Local versus nonlocal boundary-layer diffusion in a global climate model.
J. Climate
,
6
,
1825
1842
.
Kiehl
,
J. T.
,
J. J.
Hack
,
G. B.
Bonan
,
B. A.
Boville
,
D. L.
Williamson
, and
P. J.
Rasch
,
1998
:
The National Center for Atmospheric Research Community Climate Model: CCM3.
J. Climate
,
11
,
1131
1149
.
Lin
,
S-J.
, and
R. B.
Rood
,
1997
:
An explicit flux-form semi-Lagrangian shallow-water on the sphere.
Quart. J. Roy. Meteor. Soc.
,
123
,
2477
2498
.
Rasch
,
P. J.
, and
D. L.
Williamson
,
1990
:
Computational aspects of moisture transport in global models of the atmosphere.
Quart. J. Roy. Meteor. Soc.
,
116
,
1071
1090
.
Rasch
,
P. J.
, and
J. E.
Kristjánsson
,
1998
:
A comparison of the CCM3 model climate using diagnosed and predicted condensate parameterizations.
J. Climate
,
11
,
1587
1614
.
Williamson
,
D. L.
,
1988
:
The effect of vertical finite difference approximations on simulations with the NCAR community climate model.
J. Climate
,
1
,
40
58
.
Williamson
,
D. L.
,
2002
:
Time-split versus process-split coupling of parameterizations and dynamical core.
Mon. Wea. Rev.
,
130
,
2024
2041
.
Williamson
,
D. L.
, and
J. G.
Olson
,
1994
:
Climate simulations with a semi-Lagrangian version of the NCAR Community Climate Model.
Mon. Wea. Rev.
,
122
,
1594
1610
.

APPENDIX A

Hydrostatic Matrix

CAM2 supports three separate numerical methods for solving the dynamical equations: spectral Eulerian, semi-Lagrangian, and finite volume dynamics. In both the Eulerian and semi-Lagrangian dynamics the hydrostatic matrix 𝗛 in (23)–(27) is defined in terms of an “energy conversion” matrix 𝗖 as

 
formula

The matrix 𝗖 is used in the evaluation of ω = dp/dt in (2). Here 𝗛 and 𝗖 are required to obey (A2) in order to conserve energy in the conversion of potential into kinetic energy.

For the finite volume dynamics, the hydrostatic matrix 𝗛 is defined as

 
formula

The rather unusual looking diagonal (l = k) term comes from deriving the layer average value of Φ

 
formula

where we have integrated by parts and applied the hydrostatic equation between (A4) and (A5). Expanding δk(pΦ) and assuming that T is constant within layers,

 
formula

Substituting Φk = Φk+ + RkTkδk lnp from (A3),

 
formula

and the second term on the rhs of (A7) gives the k = l element of (A3).

APPENDIX B

Solution of Implicit Vertical Diffusion Equations

Equations (38)–(41) constitute a set of four tridiagonal systems of the form

 
formula

where ψn indicates u, υ, q, or s after updating from time n values with the nonlocal and boundary fluxes. The superdiagonal (Ak), diagonal (Bk), and subdiagonal (Ck) elements of (B1) are

 
formula

The solution of (B1) has the form

 
formula

Substituting (B6) into (B1) gives

 
formula

Comparing (B5) and (B7), we find

 
formula

The terms Ek and Fk can be determined upward from k = L, using the boundary conditions:

 
EL+1 = FL+1 = AL = 0.
(B10)

Finally, (B7) can be solved downward for ψkn+1, using the boundary condition:

 
C1 = 0 ⇒ E1 = 0.
(B11)

CCM1–3 used the same solution method, but with the order of the solution reversed, which merely requires writing (B6) 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 (B8)–(B9) are determined from the bottom upward, it is only necessary to recalculate Ak, Ck, Ek, and 1/(BkAkEk+1) for each constituent within the region where molecular diffusion is important.

Footnotes

Corresponding author address: Dr. Byron Boville, NCAR, P. O. Box 3000, Boulder, CO 80307-3000. Email: boville@ucar.edu

*

The National Center for Atmospheric Research is sponsored by the National Science Foundation.