the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Enforcing conservation of axial angular momentum in the atmospheric general circulation model CAM6
Thomas Toniazzo
Mats Bentsen
Cheryl Craig
Brian E. Eaton
Jim Edwards
Steve Goldhaber
Christiane Jablonowski
Peter H. Lauritzen
Numerical general circulation models of the atmosphere are generally required to conserve mass and energy for their application to climate studies. Here we draw attention to another conserved global integral, viz. the component of angular momentum (AM) along the Earth's axis of rotation, which tends to receive less consideration. We demonstrate the importance of global AM conservation in climate simulations with the example of the Community Atmosphere Model (CAM) with the finitevolume (FV) dynamical core, which produces a noticeable numerical sink of AM. We use a combination of mathematical analysis and numerical diagnostics to pinpoint the main source of AM nonconservation in CAM–FV. We then present a method to enforce global conservation of AM, and we discuss the results in a hierarchy of numerical simulations of the atmosphere of increasing complexity. In line with theoretical expectations, we show that even a crude, nonlocal enforcement of AM conservation in the simulations consistently results in the mitigation of certain persistent model biases.
 Article
(8576 KB)  Fulltext XML

Supplement
(1118 KB)  BibTeX
 EndNote
The atmosphere exchanges angular momentum (AM) with the material bodies at the surface, which are, to a good approximation, in a state of motion consisting of uniform rotation about the planetary axis connecting the poles. Per unit of mass, surface AM increases in quadratic proportion to its distance from the planetary axis of rotation, from zero at the poles to a maximum at the Equator. AM is a constant of motion of the dynamical (e.g. Newton's) equations so that as air travels meridionally, it carries a specific AM that increasingly differs from that of the Earth's surface. A variety of mechanisms redistribute atmospheric AM and eventually lead to an exchange of AM between the atmosphere and the surface, mainly as a result of lowlevel wind shear (“surface stress”) and of smallscale wave motions over steep surface topography (“form drag”).
The importance for the atmospheric circulation of the conservation of AM in the free troposphere and of AM exchange of air with the surface was recognised long ago. Already in 1735, George Hadley, Esq, F.R.S., noted that “without the Assistance of the diurnal Motion (i.e. rotation) of the Earth, Navigation (…) would be very tedious” (Hadley, 1735) due to the absence of the trade winds. This insight still lies at the core of modern conceptual models for the atmospheric circulation (Schneider, 1977; Held and Hou, 1980; Lindzen and Hou, 1988; Pauluis, 2004; Walker and Schneider, 2006). In the upper branch of the Hadley circulation (HC), the advection of planetary angular momentum determines a sharp acceleration of the zonal wind in the midlatitudes linked to a frontlike drop in air temperatures, marking the location of the subtropical jets (STJs). Partly by baroclinic instability, the midlatitude circulation redistributes atmospheric AM vertically and produces intense surface westerlies, whereby the air loses AM to the surface. The equatorward return flow in the surface branch of the HC in turn results in easterly “trade” winds, whereby surface stresses replenish atmospheric AM until air is lifted in cumulus convection within the intertropical convergence zone (ITCZ).
This circulation is the object of numerical simulations with general circulation models (GCMs) used in meteorological forecasting and in climate modelling. They describe the atmosphere as a thin, densitystratified, rotating gaseous spherical shell. These properties allow the introduction of a convenient set of approximations in the equations of motion, which result in a system known as the hydrostatic primitive equations (HPEs). The reader is referred to White et al. (2005) for a detailed analysis and discussion. Given suitable boundary conditions, the HPEs guarantee the global conservation of three fundamental physical quantities: mass, energy, and AM along the Earth's rotation axis. Analytic expressions of these laws can be found e.g. in Laprise and Girard (1990). The three conservation laws determine the fundamental character of the largescale circulation of the atmosphere, and virtually every climate application of GCMs is sensitive to their enforcement when the continuum equations are discretised in space and time. For example, the effects of changes in radiative forcing of ∼2 W m^{−2} (e.g. IPCC, 2013, chap. 8, p. 697) can only be simulated if the model's energy conservation is significantly better than 1 %. Estimates based on ECMWF reanalysis data suggest that the conservation of AM of a similar precision is desirable for an accurate representation of the annual cycle and of interannual variations of the atmospheric circulation in model simulations (e.g. Egger and Hoinka, 2005).
CAM, the Community Atmosphere Model developed and maintained at the National Center for Atmospheric Research (NCAR) in Boulder, Colorado, is one of the atmospheric general circulations models (AGCMs) in most widespread use today. It also constitutes the core atmospheric component of NorESM, the Norwegian Earth System Model. Although it offers a choice of dynamical cores, the finitevolume (FV) dynamical core (Lin, 2004) has been, and in many instances still is, the default option. The FV dynamical core is exactly mass and vorticity conserving, and it has been employed in all model integrations submitted by NCAR and by the Norwegian Climate Centre (NCC) for the 5th phase of the Coupled Model Intercomparison Project (CMIP) contributing to the Assessment Report (AR) of the Intergovernmental Panel for Climate Change (IPCC, 2013); it is also expected to be used for phase 6 of CMIP by both institutions. Due to its high numerical efficiency, FV also continues to be the code of choice for all uses for which the overall availability of supercomputing resources is a limiting factor. This includes long historical or palaeoclimate simulations, studies with coupled chemistry and/or carbon cycle, seasonaltodecadal coupled forecasts, academic research, and all model development efforts currently underway with NorESM.
In this paper, we employ CAM with the FV dynamical core at two standard CESM resolutions only, a coarser one of 1.9^{∘} × 2.5^{∘} in latitude and longitude, respectively (f19 for short), and a finer one of 0.9^{∘} × 1.25^{∘} (f09). In agreement with previous results (Lauritzen et al., 2014; Lebonnois et al., 2012), we find that all existing simulations with CAM FV, from CMIP5 to present development versions of CAM6, have a numerical sink of global AM with a magnitude of about 30 % of physical sources at f19 resolution and about 15 % at f09 resolution.
Figure 1 shows the spurious AM source in aquaplanet (AP; Neale and Hoskins, 2000; Blackburn et al., 2013) and Held–Suarez (HS; Held and Suarez, 1994) simulations with CAM FV and an otherwise identical simulation but with the global spectral dynamical core with T42 truncation. Although many other models also do not conserve AM, CAM FV is peculiar in producing a sink nearly everywhere, resulting in a particularly large global nonconservation.
First principles (e.g. Held and Hou, 1980; Einstein, 1926) suggest that the dissipation of AM, equivalent to a body force acting on the fluid as a sink of zonal momentum, forces a secondary circulation with the same sign as the Hadley circulation. As a result, the simulated Hadley circulation may become too vigorous. Reduced meridional advection of zonal momentum may lead to midlatitude westerlies that are too weak or displaced poleward. The zonal momentum lost to the nonphysical sink must be balanced by a matching additional eastward torque, for example in an expanded or excessively intense area of tropical easterly surface winds. Model simulations with CAM FV consistently tend to reflect such phenomenology: for example, Feldl and Bordoni (2016) and Lipat et al. (2017) show that among CMIP5 models, those based on the FV dynamical core (GFDLx, CCSM4, and NorESMx) simulate both relatively large overturning mass flux in the HC and a high latitude of its edge.
It is useful to illustrate these effects of AM nonconservation by means of idealised AGCM experiments that do not include complicating factors such as orographic form drag or parameterised bulk stresses associated with gravity waves. Figure 2 shows the surface torques resulting from four solutions for the mean circulation with CAM in AP mode. One of these is obtained directly from integrations of CAM using the FV dynamical core at f19 resolution (black line). An otherwise identical integration with the global spectraltransform dynamical core at T42 spectral truncation (green line) is chosen for comparison as a bona fide example of an AMconserving simulation (see Fig. 1).
The other two integrations, represented by the blue and red lines, are perturbed in an identical but opposite manner. First, the global total numerical torque due to the FV dynamical core was diagnosed at every time step of the reference FV simulation and averaged in time afterwards. This was converted into a solidbody axial rotation tendency that was applied continuously everywhere as a constant sink of AM in a new integration with the spectral dynamical core, resulting in the simulation represented by the red curve. Vice versa, the opposite additional solidbody rotation tendency was applied to a new FV integration, thus compensating for its internal numerical sink. This integration produced the physical torque represented by the blue curve. Comparing the different curves, it may be seen that equatorward of about 23^{∘} of latitude the simulated physical torque depends primarily on the global budget of atmospheric AM. In particular, notwithstanding the complications of interactive moist physics and the different spatial and temporal discretisations used in the two integrations, the stronger trade winds (in terms of surface stress) in the FV simulation compared with the T42 simulation can be explained entirely with the nonphysical, numerical torque of the FV dynamical core. The result is insensitive to how that torque is in fact applied. Even at subtropical and middle latitudes, half of the difference between the two simulations, in terms of surface stresses, can be explained in this way. Similar results are found for the zonalmean meridional circulation and for the surface pressure in the HC (Fig. S1 in the Supplement), confirming the strength and robustness of the Einstein (1926) “tea leaves” mechanism.
These results motivate us to address the issue of AM conservation in the CAM's FV dynamical core. One may speculate that systematic biases in surface stresses due to the numerical sink of AM must also impact coupled ocean–atmosphere climate simulations, with excessive Ekman and Sverdrup forcing of the subtropical gyres. The northward displacement of the midlatitude westerlies may also result in excessive mechanical and thermal forcing of the subpolar gyres with possible implications for the Atlantic meridional overturning circulation.
In this paper, we propose ways to address the numerical dissipation of AM in CAM–FV simulations. Section 2 describes our main hypotheses as to the root cause of the error and our approaches towards rectification. Section 3 presents the result of our corrections in a set of idealised simulations. The impact on realistic simulations of the atmospheric circulation is discussed in Sect. 4. Conclusions are finally offered in Sect. 5.
The FV dynamical core (Lin, 2004) solves the HPE by updating first the advective (C grid) and then the prognostic (D grid) winds in two steps. The first step represents pure advection, i.e. the increments associated with transport, including geometric and Coriolis terms. In this step, the scheme conserves absolute vorticity exactly for the Dgrid winds (Lin and Rood, 1997; hereafter LR97). The second step calculates the wind increments associated with hydrostatic pressure forces. These are computed in a special way (Lin, 1997) that differs from most Arakawa and Lamb (1981) type schemes. Violations of AM conservation may occur in either substep.
2.1 Pressuregradient force
We first analysed the Lin (1997) treatment of the pressuregradient terms for conservation. A general discussion is given by Simmons and Burridge (1981), who introduce a set of hybridlevel dimensionless variables, a_{k}, defined as ${a}_{k}:=({\mathit{\varphi}}_{k}{\mathit{\varphi}}_{k+\mathrm{1}/\mathrm{2}})/\mathrm{2}(\mathit{\alpha}p{)}_{k}$ (in Simmons and Burridge these variables are denoted by α_{k}; we change the notation here to avoid confusion), where ϕ is the geopotential, p the pressure, $\mathit{\alpha}:={\partial}_{\mathit{\eta}}\mathit{\varphi}/{\partial}_{\mathit{\eta}}p$ the specific volume, and η is the generalised or hybrid vertical coordinate. Here and in the following, the index k refers to the vertical level, or to halflevels as appropriate, and subscripts to the partial derivative symbol indicate differentiation with respect to the variable in subscript, ${\partial}_{X}\equiv \partial /\partial X$. The variables a_{k} need not be constants. Simmons and Burridge (1981) derive the discrete form that pressure and geopotential terms must take in general vertical coordinates in order to ensure conservation of axial angular momentum. Their Eq. (3.8) can be generalised to
where the symbol Δ is employed to represent a difference between vertical levels, $\mathrm{\Delta}\phantom{\rule{0.125em}{0ex}}{p}_{k}:={p}_{k+\mathrm{1}/\mathrm{2}}{p}_{k\mathrm{1}/\mathrm{2}}$ (and similarly for ϕ), and λ is the longitude.
Performing Lin's (1997) path integration around the finitevolume element on this expression yields the following form for the body force:
where δ_{λ} is the finitedifference operator in the zonal direction, and $\stackrel{\mathrm{\u203e}}{{\mathit{\varphi}}_{k\pm \mathrm{1}/\mathrm{2}}}$ is an average over λ. An expression identical in form to Lin's (1997) Eq. (11) is then recovered if the choices
are made, where i is the index corresponding to the longitude λ.
In other words, Lin's (1997) expression for the pressuregradient term is consistent with the Simmons and Burridge (1981) prescription for AM conservation, provided that the physical pressure variable p is used in the integration in place of the general pressure function indicated by the symbol π in Lin (1997). This can be directly verified algebraically by summing all expressions of the form of the numerator in the righthand side of Eq. (11) in Lin (1997) along all longitudes and levels. Provided ϕ is constant at one model boundary and p at the other, it always returns zero. This is the required result provided that the denominator on the righthand side of Eq. (11) in Lin (1997) represents the inertial mass associated with the velocity points. They do so if π is the hydrostatic pressure.
Accordingly, we performed tests in which the integration variable in the relevant section of CAM–FV's dynamical core was replaced with true interface pressure. The effect was generally seen to be very small on the dynamical core's momentum conservation properties.
We note, however, that in the CAM implementation there may be an additional problem associated with the use of the D grid. The application of Lin's (1997) method would strictly require a C grid, with zonal velocity points interleaving pressure (scalar) points along the same latitude. Thus, in CAM pressure is interpolated to the grid cell corners before use. While the formal expressions for the pressure forces do not change, thus ensuring the Simmons and Burridge (1981) total torque constraints, the inertial mass associated with each Dgrid zonal velocity point is in fact averaged over six scalar points surrounding it, with 121 weights along the zonal direction. This additional zonal smoothing effectively adds spurious terms to the zonal momentum equation of the form $u{\partial}_{x}^{\mathrm{2}}\mathrm{\Delta}p$. This is a potential source of nonconservation. However, it is not expected to be systematic.
2.2 Geometry, polar filtering, and FFSL extension
AM conservation may be affected by the treatment of geometric terms in latitude–longitude coordinates, especially near the poles where such terms become large. Furthermore, convergence of the meridians forces filtering of the solution, and additional approximations need to be made. In particular, LR97 implement a fluxform semiLagrangian extension of Colella and Woodward's (1984) piecewise parabolic method (PPM), which is used near the poles where Courant–Friedrichs–Lewy numbers (Courant et al., 1928) become large during the time integration. We performed several sensitivity tests on each of these aspects without being able to notice significant impacts on AM conservation.
Particularly compelling is the comparison with the performance of a prototype implementation in CAM of the FV scheme on a cubedsphere grid (FV3), which lacks any poles and does not require or use any of these special formulations (and is, in particular, run in pure Eulerian mode, i.e. without the fluxform semiLagrangian extension described in Lin and Rood, 1996). We ran an AP simulation on the C48 grid, viz. six pseudocubic faces with 48×48 grid cells each, for a total number of grid points identical to the standard 2^{∘} FV configuration but a 25 % higher resolution at the Equator. The AM sink (Fig. S2 in the Supplement) is nevertheless comparable, i.e. about 25 % smaller, consistent with the scaling with the resolution of simulations with standard FV. We conclude that FV and FV3 suffer from the same problem independent of geometry or the FFSL extension of LR97.
In order to minimise the impact of other minor (and partly intentional) numerical sources and sinks of AM, in all idealised numerical tests presented in this paper we applied the following modifications: (1) the order of the advection scheme is kept the same (fourth) for all model layers, instead of reducing it to first in the top layer and to second up to the eighth layer; (2) an additional conservation check is applied in the vertical remapping of zonal wind , and column momentum is conserved in the moistmass adjustment at the end of physics; (3) the surface stress residual resulting from closure of the diffusion operator (in physics) is applied in full rather than partially.
2.3 Discretisation of the kinetic energy term
The evidence from our theoretical and diagnostic analysis points at the advective, shallowwater part of the implementation of LR97 in CAM–FV as the root of the AM conservation error. Its “vectorinvariant” formulation (Arakawa and Lamb, 1981) allows for different forms of the divergence to be used in the momentum and in the mass and tracer equations, resulting in inconsistent values for the divergence of the flux of planetary AM (associated with mass divergence) and of the flux of relative AM (associated with momentum divergence). In the momentum equations, the divergence is contained in a kinetic energy (KE) gradient term, which due to the presence of a numerical symmetric instability (Hollingsworth et al., 1983) is expressed as the local gradient of a Lagrangianaverage KE. Its form violates the finitevolume approximations used for other quantities, e.g. vorticity. This feature is intrinsic to the LR97 numerical discretisation scheme and cannot be eliminated.
To address the resulting violation of AM conservation, we first note that even in AMconserving schemes, conservation can only be guaranteed in the zonal average (Simmons and Burridge, 1981). We therefore do not attempt a local correction to the scheme, which is liable to numerical instabilities (Hollingsworth et al., 1983), and instead formulate a zonalmean correction as follows. We enforce the AM conservation law,
by adding a zonalmean zonalwind tendency term to the vectorinvariant form:
Here, K is the KE plus the contribution from explicit divergence damping used in FV. In the continuum limit the expression on the righthand side simply reduces to the massweighted zonal average of the zonal gradient of $K({u}^{\mathrm{2}}+{v}^{\mathrm{2}})/\mathrm{2}$.
In discrete form, the last two terms must be approximated. In the C–D grid formulation of the LR97 scheme the second one is especially problematic. Various possibilities were explored, which resulted in various degrees of accuracy and stability. The best compromise is to discretise it as
where u is the prognostic Dgrid zonal wind and v is the advective (Cgrid) meridional wind. The details of the derivation are given in Appendix B. Using the mass conservation equation, this approximation allows us to discretise the two last terms together and write the zonalwind correction increment in a form consistent with LR97:
Here, ${\mathit{\zeta}}_{\mathit{\lambda}}:=\frac{\mathrm{1}}{a\mathrm{cos}\mathit{\phi}}{\partial}_{\mathit{\lambda}}v$, and the notation of LR97 is used for the discrete transport operators 𝒴 and ℱ for the meridional transport of ζ_{λ} and the zonal transport of mass, respectively. The first three terms in the integrand of Eq. (7) thus correspond to the first three terms (first and second lines) on the righthand side of Eq. (B10) in Appendix B. The last symbol on the righthand side of Eq. (7) represents higherorder terms (also detailed in Eq. B10). We will refer to this modification of the LR97 scheme as the “correction”.
2.4 Diagnostic tools and global conservation
Irrespective of whether the correction, as described above, is applied or not, for diagnostic purposes we calculate the apparent nonphysical torque associated with the FV dynamical core advective tendencies only, i.e. excluding the increments associated with pressure gradients. These tendencies are diagnosed separately for each layer at every advective substep and integrated horizontally to yield the apparent numerical global total torque during the substep. At the same time, the layer effective moment of inertia over the substep is also computed.
The opposite of the ratio of these quantities gives an angular acceleration that, applied to the zonal wind in each layer at every advective substep, enforces conservation of the AM of that layer under advection. The application of this solidbody rotation increment at each dynamical time step and for each layer independently is what we call the “level” fixer. The details of the computation are given in Appendix C.
Irrespective of whether they are actually applied, the fixer's velocity increments (Eq. C2), are vertically interpolated and accumulated over the entire dynamic time step and written out diagnostically. In addition to the fixer, partial wind and pressure tendencies arising from the dynamical core are separately diagnosed and written to the standard output streams, providing additional diagnostic tools for crosschecking.
A variant of the fixer was tested in CAM simulations. This variant is a “global” fixer, which still acts by applying an increment to the zonal wind at each time step. In this fixer, the apparent torque and the moment of inertia are integrated over all levels within the domain over which strict overall angular momentum conservation is desired. The zonalwind increments are then applied as a single solidbody rotational acceleration within this domain. Experimentation showed that such acceleration should not be applied in the stratosphere, where conservation errors are small and the impact of unphysical zonal accelerations large. The necessary limitation of the domain for the global fixer, however, introduces a certain degree of arbitrariness in its application. Although sometimes used for diagnostic purposes, we do not discuss this global fixer variant any further.
Lin's (2004) FV scheme conserves mass and absolute vorticity exactly. The AM modifications described above were explicitly designed not to alter the mass flux calculations and intervene only on the rotational component of the flow in the momentum equations. Other choices, involving alterations to the calculation for the divergent flow, would have been possible. However, we judged exact mass conservation more important for climate simulations than exact vorticity conservation. The AM modifications also change the kinetic energy of the flow and thus change the total energy budget of the model. However, the unmodified FV scheme does not conserve energy. CAM–FV therefore employs an energy “fixer” (analogous to out AM fixer), described e.g. in Williamson et al. (2015). The fixer diagnoses the energy nonconservation at each time step. This allowed us to monitor the impact of the AM mods on energy nonconservation in all our experiments. We found no systematic effect, either in sign or magnitude, of the AM modifications on the energy nonconservation of the model.
3.1 Dry baroclinic wave tests
Initial tests were carried out for adiabatic dynamics and flat bottom topography from baroclinically unstable initial conditions, as defined in Jablonowski and Williamson (2006; JW06). Figure 3 shows the result in terms of the conservation of global AM for CAM–FV integrations at f19 resolution (1.9×2.5^{∘} of latitude and longitude) and 30 hybrid levels.
It may be seen that both the correction and the fixer are effective in reducing the systematic numerical sink of AM in these integrations. In particular, the fixer appears to remove it almost completely; in other words, the integration with the fixer conserves global AM in the time mean. This result is central to this paper, and it proves its two main conclusions. The first is that the systematic nonconservation of global AM in the FV dynamical core indeed resides in the advective wind increments of the shallowwater part of the dynamical core. The second is that, by virtue of its effectiveness and its formulation that is entirely independent of the model configuration or parameterisations (topography, physical momentum sources, etc.), the fixer is a useful and accurate general diagnostic tool that allows us to quantify the numerical torque in any CAM–FV integration. By virtue of this quality, the diagnosed timeaveraged fixer tendencies were, for example, used for the perturbations in the experiments shown in Figs. 2 and S2.
The impact of the correction on conservation is generally smaller, and different dynamical regimes may be seen when the size and quality of that impact change. In the baroclinic instability tests in Fig. 3, the correction achieves good results in the linear and nonlinear stages of baroclinic growth (up to day 30; see JW06) but is not able to correct the slow drift that sets in after zonalisation of the global flow; then wind speed decreases everywhere as a result of numerical dissipation (there are no external sources or sinks of either momentum or energy in these adiabatic simulations). This is a partly desirable behaviour, as the action of the correction should not change the dissipation properties of the scheme.
Aside from the conservation properties they are designed for, both the correction and the fixer represent a perturbation of the numerical solutions of the FV dynamical core. By arbitrarily modifying the relative vorticity associated with the zonal wind, both destroy one of the fundamental numerical properties of the LR97 formulation, viz. the conservation of absolute vorticity under advection. (In the case of the fixer, the vorticity input has a rigid dependency on latitude, ∼sin φ.) Figure 4a shows their impact on the accuracy of the JW06 baroclinic wave test in terms of the root mean square (RMS) of the differences in surface pressure from a nominal reference solution with original FV dynamical core. The latter is obtained at f19 resolution (0.9^{∘} × 1.25^{∘}), which is sufficiently close to JW06's reference solution (see JW06, Section 5(e), points i and ii) for our purposes. It may be seen that on this measure the solutions with and without the AM corrections are virtually indistinguishable during the stages of both linear and nonlinear baroclinic growth. A similar result holds for the phase (not shown).
It may be noted that the largest impact on the RMS of surface pressure arises from the correction. Within the first 30 d this impact is formally always well below significance (as defined in JW06; see their Fig. 10), but it increases in time and eventually becomes appreciable as a full global meridional circulation is established. Similar results hold for the vorticity field, as seen in Fig. 4b.
Other aspects of the solution besides RMS differences also show limited sensitivity to the application of the correction and the fixer. Figure 5 shows the evolution of the minimum pressure in the developing baroclinic wave. By this measure, the solutions only start to diverge with the filling of the primary cyclone and the deepening of the secondary wave after day 17. The solution with the fixer deepens the secondary cyclone more quickly so that the minimum pressure is seen to jump from the first to the second wave minimum between days 18 and 19; this occurs 1 d later with the unmodified dynamical core. A third transition after day 25 has higher central pressure in the solutions with the fixer; by this time, however, rapid cyclogenesis is occurring in the jet stream of the Southern Hemisphere, attaining a similar minimum pressure, which is slightly deeper in the solutions with the fixer. In any case the pressure differences of the minima remain of the order of a few hectopascals (hPa), and there is no systematic difference in their position.
3.2 Other idealised tests
Even if the impacts of the modifications of the FV dynamical core are relatively small on local circulations over subseasonal timescales, as shown above, the rationale for introducing them is the hope of achieving a better simulation of the state of the atmosphere in integrations under specified forcings. As explained in the Introduction, one particular expectation is that the subtropical easterlies should weaken without affecting the circulation elsewhere too heavily. In particular, the role of the correction, which alone does not ensure AM conservation, must be clarified and its eventual use justified. Here we document the results of two sets of idealised simulations that still have a simplified, equipotential lower boundary but include nonvanishing physical torques and heating tendencies.
The first set of such simulations adhere to the benchmark test of Held and Suarez (1994; HS henceforth), whereby the forcing has the form of a relaxation towards a specified threedimensional atmospheric temperature field. Likewise, surface friction is represented by a damping of the winds within a set of levels near the bottom boundary. Apart from the small numerical diffusion, these stresses are communicated to the rest of the atmosphere by means of momentum advection in the mean circulation and of pressure fluctuation in resolved transient motions (including travelling waves). The second set of simulations follows the aquaplanet (AP) test first proposed by Neale and Hoskins (2000), wherein only a persistent field of bottomboundary temperatures is prescribed (the QOBS profile of Neale and Hoskins, 2000), and the full set of moist atmospheric physical parameterisations of CAM6 are used to force the circulation (except for those specific to orographic processes). The bottom boundary is a notional static ocean with unlimited heat and water capacity. Surface stresses are computed by the coupler and passed to the moist atmospheric boundarylayer parameterisation, which then distributes those stresses vertically. Momentum is also transported in moist convection, where active, and further adjustments are made when the moist mass of the atmospheric column changes due to precipitation and surface evaporation processes. To simplify the analysis, the gravitywave parameterisation of CAM6 was turned off in our AP tests. In both sets of tests, FV's advection scheme is used at PPM's standard fourth order at all levels; i.e. the numerical diffusion obtained in standard CAM–FV integrations by employing loworder calculations near the model top is avoided. For initial conditions, HS simulations are coldstarted with uniform surface pressure and geopotential, as well as vanishing wind fields except for a westerly perturbation identical to that used in the dry baroclinic wave tests (necessary in order to break zonal symmetry and to allow for a nonvanishing correction). The AP simulations all take the same instantaneous atmospheric state from a previous spunup run, even though this requires more adjustment for the corrected (fixed) simulations than for the control.
Figure 6a indicates that the global AM conservation properties of the simulations in these tests are broadly in line with the expectations from the previous discussion. Standard FV tests (black lines) show a steady loss of AM in the atmospheric circulation, with a magnitude of the order of 10 %–20 % of the physical flux of AM through the atmosphere. (We count eastward stress as positive, by which the atmosphere gains westerly momentum in the tropical surface easterlies and loses westerly momentum in the subtropical surface westerlies.) Use of the correction leads to an orderofmagnitude reduction of the numerical sink of AM in HS integrations, but it is of limited effectiveness in fullphysics AP integrations (blue lines). Integrations with the fixer, with or without the correction (orange and red lines, respectively), maintain atmospheric AM in the time mean. In HS simulations, there appears to be a very small residual drift of AM notwithstanding the fixer. This is due to a small inconsistency in the application of the stress terms, which are calculated and diagnosed in the “physics” part of the model time stepping but applied later as velocity tendencies in the physics–dynamics interface on updated layer masses. This is an intrinsic feature of the time stepping of CAM–FV that we have not modified. More notably, AP simulations differ from HS simulations in that they show obvious fluctuations of total AM around the time mean or around the longterm drift when there is one. Such fluctuations are similar in all AP integrations, with a magnitude of a few percent of the physical sources, and depend on nonconservation in CAM's physics parameterisations. Fortunately, they are not systematic and do not produce a noticeable longterm drift.
The effectiveness of the fixer in removing most of the AM drift confirms that the systematic sink of AM in CAM–FV integrations arises predominantly from the shallowwater advection calculations. The accuracy of the correction, by contrast, depends on the features of the circulation, with good accuracy for numerically wellresolved features, as in the HS tests, but a poorer one when gridscale forcing associated with the water cycle occurs. Figure 6b gives more details on the effect of the correction. Here, the timeaverage AM sink due to the dynamical core is diagnosed using the fixer increments for the zonal velocity at the Equator at each model level. This diagnostic is produced irrespective of whether such increments are applied during the integration. Apart from the smaller increments in HS integrations than in AP integrations, which partly depend on the slower circulation (“surface” stresses are 1 order of magnitude larger in the HS setup than in the AP setup), the advective AM sink has a distinctive shape in pressurelevel space, with a maximum in the upper troposphere and small values in the atmospheric boundary layer. This shape partly reflects the underlying global mean zonalwind field, but the maximum sink lies below the maximum wind (at around 250 hPa rather than around 150 hPa). The profile of the impact of the correction, i.e. the reduction in fixer increments when the correction is applied, has again a similar shape but with an even lower position of the maximum, which better corresponds to the maximum in the vertical profile of the levelintegral zonal momentum of the underlying flow. Combined with the offline diagnostic information for the apparent AM sink from Fig. 1, it can be deduced that the main loci of the timemean AM sink in these simulations are found near the subtropical jet streams, where large zonal asymmetries occur in both the mass fields and the wind fields.
The effect on the mean circulation of applying the correction and/or the fixer is shown in Figs. 7 and 8 for HS and AP simulations, respectively. The zonalmean zonal winds are shown, which is the quantity that both the correction and the fixer directly modify. Nevertheless, it should be remembered that the net effect is indirect, since the zonal winds remain in the time average close to geostrophic balance with the (equivalent) temperature field. In HS simulations, the local temperature differences between simulations are simply proportional to the difference in temperature advection by the meridional and vertical circulation, which is modified primarily through a “tea leaves” mechanism. As already seen in the Introduction, the leadingorder effect of the fixer is a weakening of this circulation and thus of the associated advective temperature tendencies. These tend to cool the lower troposphere in the subtropical easterlies, cool the upper troposphere near the Equator, and warm the troposphere poleward of the jet streams. The effect of the fixer on the zonalmean zonal wind shown in Fig. 7a is generally consistent with this expectation, with an equatorward retreat of the surface easterlies and weaker westerlies in the higher latitudes. There is, however, an additional large westerly difference near the equatorial tropopause; this is a direct consequence of the westerly forcing of the fixer, which is greatest at the Equator. This is clearly an undesirable effect of the fixer on the simulations. A more selective effect on the circulation is produced by the correction (Fig. 7b). As seen above, its main action is in where the greatest sink of AM is located, i.e. on the flanks of the subtropical jet stream. By correcting part of the AM nonconservation, it also acts to limit the action of the fixer (Fig. 7d). As a result, the combination of the correction and fixer together, as well as ensuring good global AM conservation, is less severe in terms of its upperlevel equatorial westerly effect (Fig. 7d). This suggests that the fixer is best employed in combination with the correction.
In AP simulations, a slowdown of the meridional circulation is still expected and found, but the interaction between dynamical forcing by the fixer or the correction and the physics tendencies is much more complex and difficult to predict. The fixer now produces large westerly differences near the Equator at all levels and a marked weakening of the subtropical jet stream (Fig. 8a). The equatorial winds above 300 hPa become westerly. The correction is less effective overall than in HS simulations, and its impacts are mostly confined to levels close to the model lid or to the high latitudes (Fig. 8b). Nonetheless, its use is still beneficial in terms of limiting the action of the fixer, at least in the troposphere (Fig. 8d). The result of the combined correction and fixer can be seen in Fig. 8c. In terms of tropospheric impacts, it appears acceptable; equatorial winds remain easterly below 200 hPa and weak above. The weakening of the equatorial and tropical easterlies compared with the control simulation implies greater similarity with simulations with AMconserving spectral models. Large changes, however, can be seen near the model lid, especially in the four model layers with pressures less than 25 hPa. This is a consequence of momentum accumulation within these layers. In CAM's default configuration, the order of FV's PPM advection scheme is reduced here, which results in large numerical dissipation. Effectively, these levels are used as sponge layers and are thus not part of the valid computational domain of the model. In fullmodel configurations it is therefore advised to keep the reduced order of advection and turn off both the correction and the fixer in these layers. The large meanstate changes seen near the top in Fig. 8d then vanish. Considering the troposphere only, the conclusion obtained from HS simulations can be seen to also hold for fullphysics AP model simulations in that the combined application of the fixer and the correction results in smaller overall meanstate changes in the solution compared to default FV without modifications, while ensuring good conservation of AM.
The relevance of the AM modifications to the FV dynamical core for CAM simulations in a realistic configuration is investigated here using F2000 cases, which are AMIPtype simulations (Gates, 1992) wherein sea surface temperatures (SSTs) and all compositional forcings are prescribed as a repeating annual cycle obtained from an observed climatology of the decade spanning the turn of the century. We test at two grid resolutions, one of 1.9^{∘} × 2.5^{∘} (f19) as in all integrations already discussed above and one of 0.9^{∘} × 1.25^{∘} (f09) to test the impacts of AM modifications in a case that is scientifically supported by NCAR at this time. The CESM model version used (here as above) is release 2.1.1^{1}.
Figure 9 illustrates the effects of the fixer and the correction on f19 simulations. The control simulation shows a characteristic easterly surface windstress bias throughout the tropics (Fig. 9a). In addition, there are excessive westerlies at southern high latitudes. The effect of the fixer is to reduce the tropical biases (Fig. 9b), with an evident westerly effect on the simulations nearly symmetrically about the Equator (Fig. 9d). By that same token, however, the highlatitude westerly errors are enhanced (Fig. 9b). The application of the correction in addition to the fixer not only brings further improvements in the tropics, but also corrects the westerly effect of the fixer in high latitudes (Fig. 9e). The result is a significant improvement in the simulation of the surface windstress field over the entire ocean domain.
In general, we obtain a similar conclusion as for the AP simulations. The impact of the correction on the global conservation of AM is modest, removing only about 15 % of the sink at f19 resolution. However, its action is stronger on upperlevel winds (see Fig. 6b), which leads to proportionally reduced fixer increments at those levels and thus to smaller impacts by the fixer on areas affected by baroclinic instability.
Figures 10 and S3 in the Supplement show the seasonally resolved impacts on the zonalmean zonal winds from applying the combination of the fixer and correction in F2000 simulations at both f19 and f09 resolutions (see also Fig. S3 in the Supplement for JJA). In all cases, the reduction of biases in both easterly and westerly wind regimes is noticeable, the latter especially at the subpolar latitudes of the winter hemisphere.
More in detail, it may be noted that the benefits of the AM modifications appear more clearly for the winds in the simulation at the lower resolution, for which the numerical sink of AM is indeed larger. These benefits, however, are not limited to the zonalmean zonal winds, and they are also appreciable at the f09 resolution. Most notable is the reduction in the strength of the Hadley circulations (see Fig. S4 in the Supplement), which is expected from the arguments set out in the Introduction. This has consequences for many aspects of the global circulation. Figure 11 shows a summary of the impacts on the quality of the simulations in relation to the observed climatology. The improvements at f09 seem particularly remarkable considering that the unmodified simulation is a scientifically supported case that has been fully tuned for a best match to observations. It may be noted that no additional tuning whatsoever is involved in the simulation with AM modifications shown here and that the AM modifications themselves have no free parameters as they follow directly from an effort to reduce the numerical sink stemming from the FV dynamical core. The better quality of this simulation thus follows entirely from better adherence of the solution to a fundamental property of the equations of motion.
AM conservation in CAM–FV has been substantially improved by means of a correction that reduces the zonalmean numerical sink of Lin and Rood's (1997) shallowwater scheme and a fixer that ensures the conservation of global angular momentum under advection. The effectiveness of these modifications in terms of AM conservation in the simulations presented here is summarised in Table 1. We show that aside from global AM conservation, they have other significant impacts on the simulations, consistent with the “tea leaves” mechanism (Einstein, 1926) that rapidly redistributes pressure forces in a rotating fluid in response to zonal accelerations. The most notable effect is a reduction of the excessive easterlies of the model, with a concomitant slowdown of the Hadley circulation. As a result of such changes, the simulations of the observed climatology show marked improvements.
The zonalmean correction of the shallowwater scheme is not necessary for enforcing global conservation, as this can be achieved by the fixer alone. Indeed, the correction is quite ineffective in realistic simulations of the atmosphere in terms of global conservation. Nevertheless, we find that its concomitant application with the fixer has positive impacts on the simulations. In particular, it reduces the effects of the fixer in the midlatitudes. This can be explained with the greater effectiveness of the correction in the baroclinically unstable regions around the subtropical jet streams, where the zonalmean numerical sink appears to be largest. Even so, because of its potentially large local effects, the utilisation of the correction under different setups should be tested on a casebycase basis according to its impacts on the results.
Improving the quality of the simulation of the global distribution of surface wind stress should be expected to bring particular benefits to coupled atmosphere–ocean simulations. An adequate discussion of such a coupled simulation would exceed the scope of the present paper, which is aimed primarily at presenting the method. In particular, due to their computational expense, at the present time it is not possible to produce well spunup coupled simulations that can provide an assessment of the impact of the AM modifications.
The modification to the FV dynamical core that we describe and utilise is relatively crude and causes local loss of accuracy due to violation of vorticity conservation under advection. Nevertheless, the associated detrimental impacts appear to be fairly limited, with insignificant differences under standard tests such as the Jablonowski and Williamson (2006) baroclinic wave test, which should be sensitive to local conservation. Even so, it is clear from the very same tests that simulations over weather timescales are not sensitive to AM conservation, so for such an application it is not advisable to trade enforcing such conservation for a loss of accuracy. On the longer timescales of climate simulations, by contrast, our results demonstrate the importance of the global conservation of atmospheric AM in order to obtain a realistic global circulation.
The diagnosis of the residual torque that violates AM conservation in CAM simulations follows from the hydrostatic primitive equations (see White et al., 2005). In our zonally and vertically integrated diagnostics such as in Fig. 1 the AM source is calculated as
where the first term on the righthand side represents the tendency of relative atmospheric AM, the second term represents the divergence of the flux of relative AM, the third the external torque (which in all simulations presented in Sects. 1, 2, and 3, when nonvanishing, is exclusively due to surface stresses or linear friction in the planetary boundary layer – PBL), and the last term is the tendency of planetary atmospheric AM due to the vertically integrated divergence of atmospheric mass. The formulas are as follows.
Here, a is the Earth's radius, φ the latitude, λ the longitude, g the gravitational acceleration in Earth's surface, Ω the angular speed of Earth's rotation, and u, v, p_{*}, and τ_{x} are the zonalwind component, the meridional wind component, the surface pressure, and the zonal component of the surface or frictional stress acting on the air in the model simulations. Note that to obtain C_{λ} the continuity equation was used. Note that for the timeaverage values of S_{M}, the time differentials become increments between the initial and the final state; terms T_{x} and C_{λ} are linear in the wind stress and the surface pressure, respectively. Terms L_{r} and D_{L} are bilinear and trilinear in the model prognostic quantities u, v, and p_{*}, so an online computation of the time averages of the integrands is required for these terms. CAM provides timemean diagnostics of the zonal wind u and of the product of the wind components uv conservatively interpolated onto standard pressure levels, and the integrals in Eq. (A1) are computed with their help.
The local conservation equation for the shallowwater equations is
where φ and λ are latitude and longitude, respectively, Δp is the layer thickness in terms of hydrostatic pressure, u and v are the zonal and meridional wind components, a is the Earth's radius, and Ω the Earth's angular velocity. Note that we are ignoring pressure and geopotential terms here, as we focus exclusively on the process of advection. Accordingly, Δp, i.e. the layer under consideration, may be arbitrary, except that it satisfies the shallowwater mass conservation equation; i.e. we follow Lin's (2004) “vertically Lagrangian” approach by following the vertical motion of the layer. Integrating Eq. (B1) over longitude, we obtain
where f is the Coriolis parameter. To address the FV scheme's violation of this conservation, we apply an additional, zonally uniform increment of the zonal wind, $\stackrel{\mathrm{\u203e}}{\mathit{\delta}u}$, such that over each shallowwater substep δt (we shall refer to this simply as the “time step” in this section) of the dynamical core,
Here, “old” prognostic quantities (i.e. valid at the beginning of the time step) and “new” prognostic quantities (i.e. valid at the end of the time step, before any correction) are indicated by the subscripts “o” and “n”, respectively; quantities without subscripts are intended to be timecentred, representing advective fluxes over the time step. To obtain the correction, we solve this equation for the required increment $\stackrel{\mathrm{\u203e}}{\mathit{\delta}u}$ and substitute for u_{n} the actual FV zonalwind increment over the time step:
where ξ is the absolute vorticity, and K is the kinetic energy term as discretised in LR97's scheme. The result is
The term in the second line on the righthand side representing the advection of planetary vorticity is written in a roundabout way for later convenience.
We note two aspects of this expression. First, there is a significant numerical cancellation between the second and the third lines on the righthand side. Second, all advective terms in the first two lines on the righthand side can be easily discretised according to the standard LR97 prescription and are thus automatically defined on Dgrid zonal velocity points, i.e. where required for $\stackrel{\mathrm{\u203e}}{\mathit{\delta}u}$. However, all mass factors are defined on scalar points, i.e. on the A grid. Furthermore, the integrand in the third line on the righthand side has no natural expression in LR97's discretisation, and both zonal and meridional winds in that expression need to be interpolated onto the A grid. Hence, additional interpolation is required for these terms. Notwithstanding these issues, we found that this correction, when implemented, gave accurate conservation of AM. However, it also proved to cause numerical instability such that the integration crashed within seven or eight time steps. Analysis suggested that the last term on the righthand side had to be recast in a different form.
We therefore chose to approximate the last term as follows:
The approximation here consists of using Cgrid (advective) fluxes in the partial differentials on the righthand side. Considering this as a calculation for the advective fluxes of zonal momentum, which is its physical meaning, this appears to be a valid interpretation for v. For the values of Δp and u outside the operators, we adopt the substitutions
where
and δ^{′′}u and δ^{′′}Δp are formally o(δt). The increments are still understood as advective only, i.e. they exclude pressure force terms. By further using the identities
we finally arrive at the expression for our approximate angularmomentumconserving zonalmean zonalwind correction:
where we have used the shorthand ${\mathit{\zeta}}_{\mathit{\lambda}o}:=\frac{\mathrm{1}}{a\mathrm{cos}\mathit{\phi}}{\partial}_{\mathit{\lambda}}{v}_{\mathrm{o}}$.
We note that setting the higherorder terms to zero implies that the correction has no effect on a zonally symmetric flow. If, in addition, the flow is in an exact steady state, then the correction always vanishes identically, regardless of these terms. It can further be shown that, if the term in K is the true gradient of the kinetic energy in the original scheme, for any values of δ^{′′}u and δ^{′′}Δp that are first order in δt or higher, the correction Eq. (B10) is formally third order in δt or higher. In other words, the correction will not affect solutions that are already locally angular momentum conserving.
In Eq. (B10), all mass terms must be averaged over φ; by contrast, all advective terms (in square brackets) represent fluxes as discretised according to the standard LR97 algorithm. The discretised expression of Eq. (B10) thus corresponds to Eq. (7). The only additional PPM calculation required to calculate this correction is the meridional advection of the partial relative vorticity, ζ_{λ}, with a minimal additional computational cost that is hardly detectable in CAM simulations.
As we explain in Sect. 2.4, the fixer is based on diagnosing the global change in atmospheric AM due to advective increments only, which should vanish identically according to the continuous equations. When applied, the fixer counteracts that change at every advective substep; irrespectively, its timemean increments can always be used to diagnose AM nonconservation in the simulations in a manner that is completely independent of the physics parameterisations or boundary conditions used and hence independent of the particular configuration of the simulations itself. All the calculations related to the fixer and the quantification of the numerical (advective) AM source are internal to the dynamical core only, indeed of its shallowwater part.
So, for each time step and at each level k, we require the advective shallowwater equation increments to satisfy
where the indices i and j refer to longitude and latitude, respectively, e_{j} represents the latitudes of the zonal velocity points of the D grid, and c_{j} represents the latitudes of the scalar points (A grid). The other symbols have the same meaning as in the previous section, and δ represents the purely advective increment obtained in the dynamical core, which may include the correction discussed above. The action of the fixer in this context is represented by an additional increment δϖ_{k} so that the total increment of the zonal wind becomes $\mathit{\delta}{u}_{i,j,k}+a\mathit{\delta}{\mathit{\varpi}}_{k}\mathrm{cos}{e}_{j}$. We obtain the following:
where the numerical torque is
and the moment of inertia is
In these expressions,
Equation (C2) gives the required angular acceleration of the entire atmospheric shell at model level k. The action of the level fixer is therefore to add an increment to the zonal wind:
In some regions of the model domain, it is not desirable to apply a fixer, since dissipation is explicitly built into in the dynamical core formulation. This is the case near the upper boundary of CAM's domain (the lower boundary in pressure space), where the fixer is accordingly switched off. In general, a weight w_{k}≤1 can be applied at each level so that Eq. (C2) becomes
where only a fraction w_{k} of the numerical torque at level k is compensated for by the fixer at that level.
The global fixer applies the same solidbody rotation increment to all levels within the domain where it is required. When all weights are unity, this is simply
when $\exists k:{w}_{k}<\mathrm{1}$, the vertical integrals must be weighted accordingly and the weights applied to the correction at each level so that
It can be seen that ${\sum}_{k}{I}_{k}\mathit{\delta}{\mathit{\varpi}}_{g,k}={\sum}_{k}{w}_{k}{T}_{k}$ so that the numerical torque associated with the domain of interest is also fully compensated for by this fixer. Experimentation has shown that tapering the global fixer so as to exclude its action from levels in the stratosphere was necessary in order to avoid distortions of the dynamics in layers where it is sensitive to small amounts of zonal acceleration and where, moreover, thanks to the predominance of solenoidal dynamics (before gravitywave drag, which is applied in the physics parameterisations), the dynamical core performs well in terms of AM conservation. For the latter reason, no tapering (i.e. any weights other than 1 in the valid domain and 0 in the filtered layers near the model lid) is in fact required for the level fixer.
For diagnostic purposes, fixer increments are always calculated as in Eq. (C2) and provided in output. Use of the increments in Eq. (C2) leads to the conservation of total AM in idealised spinup or spindown experiments with no physical sources or sinks of momentum (see Fig. 3), as well as an accurate balance of the surface torques in Held–Suarez or aquaplanet simulations wherein only surface stresses are present (and accurately diagnosed). Hence, we obtain two important conclusions. First, all numerical sources of AM indeed reside in the advective wind increments of the shallowwater part of the dynamical core; second, the fixer diagnostics return an accurate estimate of the apparent numerical AM source for any CAM–FV integration, irrespective of physics parameterisations or boundary fluxes (including orographic form drag).
The code used in the numerical simulations of this paper is available under: https://doi.org/10.5281/zenodo.3529202 (Goldhaber et al., 2019). CAM6 is published in the openaccess CESM ESCOMP git repository, which is freely available under: https://github.com/ESCOMP (Goldhaber and Craig, 2020). The AM options can be switched on by setting standard CAM namelist parameters to nondefault values (i.e. T instead of F; there are no free numerical parameters). Apart from these switches, all atmosphere model configurations presented in this paper are standard CESM cases that can be set up and run using the scripts provided in the repository. Users can obtain technical support if requested.
The supplement related to this article is available online at: https://doi.org/10.5194/gmd136852020supplement.
TT conceived the idea, proposed the work, made the calculations, implemented the code, ran the simulations, evaluated them, produced all figures, and wrote the paper. MB supported this activity through national infrastructure projects of the Norwegian Research Council. CC, BEE, and JE revised the code and included it in the official ESCOMP CESM repository. SG gave technical advice on CAM code and simulations. PHL, MB, and CJ were at hand for critical discussions of the scientific ideas and helped provide the initial impetus for this work. MB, JE, SG, and PHL also provided useful comments and suggestions on the draft paper.
The authors declare that they have no conflict of interest.
Warm thanks go to Christoph Heinze for his unbending dedication to model development in the NorESM consortium and in particular for allowing this work to go forward despite the lack of dedicated funding. We are grateful to Alok Gupta at NORCE and Cecile Hannay at NCAR for their assistance with NorESM and CESM development simulations.
This work was partially supported by the Norwegian Research Council (EVA (grant no. 229771) and INES (grant no. 270061)) as well as by the US National Science Foundation Cooperative Agreement (grant no. 1852977).
This paper was edited by Paul Ullrich and reviewed by two anonymous referees.
Arakawa, A. and Lamb, V. R.: A potential enstrophy and energy conserving scheme for the shallow water equations, Mon. Weather Rev., 109, 18–36, 1981.
Blackburn, M., Williamson, D. L., Nakajima, K., Ohfuchi, W., Takahashi, Y. O., Hayashi, Y.Y., Nakamura, H., Ishiwatari, M., McGregor, J. L., Borth, H., Wirth, V., Frank, H., Bechtold, P., Wedi, N. P., Tomita, H., Satoh, M., Zhao, M., Held, I. M., Suarez, M. J., Lee, M.I., Watanabe, M., Kimoto, M., Liu, Y., Wang, Z., Molod, A., Rajendran, K., Kitoh, A., and Stratton, R.: The AquaPlanet Experiment (APE): CONTROL SST simulation, J. Meteorol. Soc. Jpn., 91A, 17–56, https://doi.org/10.2151/jmsj.2013A02, 2013.
Colella, P. and Woodward, P. R.: The piecewise parabolic method (PPM) for gasdynamical simulations, J. Comput. Phys., 54, 174–201, 1984.
Courant, R., Friedrichs, K., and Lewy, H.: Über die partiellen Differenzengleichungen der mathematischen Physik, Math. Ann., 100, 3274, https://doi.org/10.1007/BF01448839, 1928.
Egger, J. and Hoinka, K.P.: The annual cycle of the axial angular momentum of the atmosphere, J. Climate, 18, 757–771, 2005.
Einstein, A.: Die Ursache der Mäanderbildung der Flußläufe und des sogenannten Baerschen Gesetzes, Die Naturwissenschaften, 11, 223–224, 1926.
Feldl, N. and Bordoni, S.: Characterizing the Hadley circulation response through regional climate feedbacks, J. Climate, 29, 613–622, https://doi.org/10.1175/JCLID150424.1, 2016.
Gates, W. L.: AMIP: The Atmospheric Model Intercomparison Project, B. Am. Meteorol. Soc., 73, 1962–1970, 1992.
Goldhaber, S. and Craig, C.: CAM: The Community Atmosphere Model, available at: https://github.com/ESCOMP/CAM, last access: 17 February 2020.
Goldhaber, S., Craig, C., and Toniazzo, T.: tto061/CAM: v6.1.029.GMD_AM_paper (Version cam6_1_029), Zenodo, https://doi.org/10.5281/zenodo.3529202, 2019.
Hadley, G.: Concerning the cause of the general tradewinds, Philos. T. R. Soc., 39, 58–62, 1735.
Held, I. M. and Hou, A. Y.: Nonlinear axially symmetric circulations in a nearly inviscid atmosphere, J. Atmos. Sci., 37, 515–533, 1980.
Held, I. M. and Suarez, M. J.: A proposal for the intercomparison of the dynamical cores of atmospheric general circulation models, B. Am. Meteorol. Soc., 75, 1825–1830, 1994.
Hollingsworth, A,, Kållberg, P., Renner, V., and Burridge, D. M.: An internal symmetric computational instability, Q. J. Roy. Meteor. Soc., 109, 417–428, 1983.
IPCC: Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change, edited by: Stocker, T. F., Qin, D., Plattner, G.K., Tignor, M., Allen, S. K., Boschung, J., Nauels, A., Xia, Y., Bex, V., and Midgley, P. M., Cambridge University Press, Cambridge, UK and New York, NY, USA, 1535 pp., 2013.
Jablonowski, C. and Williamson, D. L.: A Baroclinic Instability Test Case for Atmospheric Model Dynamical Cores, Q. J. Roy. Meteor. Soc., 132, 2943–2975, 2006.
Laprise, R. and Girard, C.: A spectral general circulation model using a piecewiseconstant finiteelement representation on a hybrid vertical coordinate system, J. Climate, 3, 32–52, 1990.
Lauritzen, P. H., Bacmeister, J. T., Dubos, T., Lebonnois, S., and Taylor, M. A.: HeldSuarez simulations with the community atmosphere model spectral element (CAMSE) dynamical core: a global axial angular momentum analysis using Eulerian and floating Lagrangean vertical coordinates, J. Adv. Model. Earth Sy., 6, 129–140, https://doi.org/10.1002/2013MS000268, 2014.
Lebonnois, S., Covey, C., Grossman, A., Parish, H., Schubert, G., Walterscheid, R., Lauritzen, P. H., and Jablonowski, C.: Angular momentum budget in general circulation models of superrotating atmospheres: A critical diagnostic, J. Geophys. Res., 117, E12004, https://doi.org/10.1029/2012JE004223, 2012.
Lin, S. J.: A finitevolume integration method for computing pressure gradient force in general vertical coordinates, Q. J. Roy. Meteor. Soc., 123, 1749–1762, 1997.
Lin, S.J.: A “vertically lagrangian” finitevolume dynamical core for global models, Mon. Weather Rev., 132, 2293–2307, 2004.
Lin, S. J. and Rood, R. B.: Multidimensional ﬂuxform semiLagrangian transport schemes, Mon. Weather Rev., 124, 2046–2070, 1996.
Lin, S. J. and Rood, R. B.: An explicit fluxform semiLagrangian shallowwater model on the sphere, Q. J. Roy. Meteor. Soc., 123, 2477–2498, 1997.
Lindzen, R. S. and Hou, A. Y.: Hadley circulations for zonally averaged heating centered off the equator, J. Atmos. Sci., 45, 2416–2427, 1988.
Lipat, B. R., Tselioudis, G., Grise, K. M., and Polvani, L. M.: CMIP5 models' shortwave cloud radiative response and climate sensitivity linked to the climatological Hadley cell extent, Geophys. Res. Lett., 44, 5739–5748, https://doi.org/10.1002/2017GL073151, 2017.
Neale, R. B. and Hoskins, B. J.: A standard test for AGCMs including their physical parameterizations. II: Results for the Met Office model, Atmos. Sci. Lett., 1, 108–114, https://doi.org/10.1006/asle.2000.0024, 2000.
Pauluis, O.: Boundary layer dynamics and crossequatorial Hadley circulation, J. Atmos. Sci., 61, 1161–1173, 2004.
Quilfen, Y., Chapron, B., Bentamy, A., Gourrion, J., Elfouhaily, T., and Vandemark, D.: Global ERS1/2 and NSCAT observations: Upwind/crosswind and upwind/downwind measurements, J. Geophys. Res., 104, 11459–11469, 1999.
Schneider, E. K.: Axially symmetric steadystate models of the basic state for instability and climate studies. Part II: Nonlinear calculations, J. Atmos. Sci., 34, 280–297, 1977.
Simmons, A. J. and Burridge, D. M.: An energy and angularmomentum conserving vertical finitedifference scheme and hybrid vertical coordinates, Mon. Weather Rev., 109, 758–766, 1981.
Taylor, K. E.: Summarizing multiple aspects of model performance in a single diagram, J. Geophys. Res., 106, 7183–7192, 2001.
Uppala, S. M., Kållberg, P. W., Simmons, A. J., Andrae, U., Bechtold, V. D. C., Fiorino, M., Gibson, J. K., Haseler, J., Hernandez, A., Kelly, G. A., Li, X., Onogi, K., Saarinen, S., Sokka, N., Allan, R. P., Andersson, E., Arpe, K., Balmaseda, M. A., Beljaars, A. C. M., Berg, L. V. D., Bidlot, J., Bormann, N., Caires, S., Chevallier, F., Dethof, A., Dragosavac, M., Fisher, M., Fuentes, M., Hagemann, S., Hólm, E., Hoskins, B. J., Isaksen, L., Janssen, P. A. E. M., Jenne, R., Mcnally, A. P., Mahfouf, J. F., Morcrette, J. J., Rayner, N. A., Saunders, R. W., Simon, P., Sterl, A., Trenberth, K. E., Untch, A., Vasiljevic, D., Viterbo, P., and Woollen, J.: The ERA40 ReAnalysis, Q. J. Roy. Meteor. Soc., 131, 2961–3012, https://doi.org/10.1256/qj.04.176, 2005.
Walker, C. C. and Schneider, T.: Eddy influences on Hadley circulations: simulations with an idealized GCM, J. Atmos. Sci., 63, 3333–3350, 2006.
White, A. A., Hoskins, B. J., Roulstone, I., and Staniforth, A.: Consistent approximate models of the global atmosphere: shallow, deep, hydrostatic, quasihydrostatic, and nonhydrostatic, Q. J. Roy. Meteor. Soc., 131, 2081–2107, 2005.
Williamson, D. L., Olson, J. G., Hannay, C., Toniazzo, T., Taylor, M., and Yudin, V.: Energy considerations in the Community Atmosphere Model (CAM), J. Adv. Model. Earth Sy., 7, 1178–1188, https://doi.org/10.1002/2015MS000448, 2015.
More precisely, we used a prerelease of CESM2.1.1 (no. 20, 22 March 2019). In terms of the simulations presented in this paper, the differences with the full 2.1.1 release only affect the F2000 cases at f19 resolution, for which slightly different emission datasets are used to force the simulations. The impacts of this are of negligible consequence for the results discussed in this section.
 Abstract
 Introduction
 Analysis of potential causes and approaches to correction
 Numerical simulations and results
 Simulations of the observed climatology
 Summary and conclusions
 Appendix A: Offline diagnostics of numerical torque in model simulations
 Appendix B: Formulation and approximations for the AM correction in CAM–FV
 Appendix C: Formulation and implementation of the AM fixer in CAM–FV
 Code and data availability
 Author contributions
 Competing interests
 Acknowledgements
 Financial support
 Review statement
 References
 Supplement
 Abstract
 Introduction
 Analysis of potential causes and approaches to correction
 Numerical simulations and results
 Simulations of the observed climatology
 Summary and conclusions
 Appendix A: Offline diagnostics of numerical torque in model simulations
 Appendix B: Formulation and approximations for the AM correction in CAM–FV
 Appendix C: Formulation and implementation of the AM fixer in CAM–FV
 Code and data availability
 Author contributions
 Competing interests
 Acknowledgements
 Financial support
 Review statement
 References
 Supplement