Potential-based thermodynamics with consistent conservative cascade transport for implicit large eddy simulation: PTerodaC³TILES version 1.0

A new computational fluid dynamics code for large eddy simulation (LES) of the atmospheric boundary layer and convection is presented and made available. A key novelty is that moist thermodynamics is formulated in terms of thermodynamic potentials, ensuring thermodynamic consistency. Despite the apparent complexity of the thermodynamic potential approach, the model's performance demonstrates that it is feasible and effective at reasonable computational cost for three-dimensional simulations. Semi-implicit semi-Lagrangian numerical methods are used; such methods are unusual for simulating boundary layer and convective flows and are more typical of global atmospheric models. Moreover, the model includes no explicit scheme to represent subgrid-scale fluxes of scalars and momentum but relies instead on the mixing and dissipation resulting from the numerical methods used; in other words, it employs implicit LES (ILES). Sample results from several standard LES test cases show that the model's ability to capture the main aspects of the flows is comparable to other LES models. At the same time, the results highlight limitations of the ILES approach near the bottom boundary and suggest that ILES might need to be augmented in some way, for example, by distributing the convergence of surface fluxes over several model layers. Also, results for a marine stratocumulus case show a significant sensitivity to different options for the numerical methods and parameters used. Further development and application of the code would benefit from a deeper understanding of both the bottom boundary behaviour and the sensitivities to numerics.

1 Introduction

Large eddy simulation (LES) has been an invaluable computational tool in atmospheric science since the early 1970s, both for advancing our understanding of complex atmospheric processes such as boundary layer turbulence and convection and for informing the development of parameterizations of those processes for use in weather and climate models. This article presents the formulation of a new LES code developed by the author with a threefold motivation:

1. to demonstrate the feasibility of using thermodynamic potentials to achieve a consistent representation of moist thermodynamics in a three-dimensional fluid dynamics code,

2. to evaluate the use in LES of the sort of numerical methods more usually used in synoptic- and global-scale models, and

3. to provide a modelling tool that is easy to set up and run and in which it is easy to set up new test cases and diagnostics.

The thermodynamics of moist air is complicated, and atmospheric models often make approximations. Some common approximations introduce inconsistencies between different aspects of the system or with the laws of thermodynamics – see Thuburn (2017b) for examples. Thermodynamic consistency can be ensured by deriving all thermodynamic quantities from one of the four standard thermodynamic potentials, internal energy, enthalpy, Helmholtz free energy, or Gibbs function, expressed as a function of its natural variables. Provided any approximations are made directly to the thermodynamic potential before deriving other quantities, consistency is maintained. This approach has been proposed for use in ocean modelling for consistent treatment of the complex equation of state of seawater (IOC et al., 2010). More recently, the approach has been advocated for deriving consistent equation sets for the thermodynamics of moist air in atmospheric models (Vallis, 2017; Eldred et al., 2022; Staniforth, 2022). As an alternative, a model can be formulated and implemented directly in terms of the thermodynamic potential and its derivatives (Thuburn, 2017b; Bowen and Thuburn, 2022a, b); in this way different approximations to the thermodynamics can be implemented by modifying the minimal number of routines.

For the case of thermodynamic equilibrium (no supersaturation, no condensate in subsaturated air, all components and phases at the same temperature) and in the absence of ice, Thuburn (2017b) presented a semi-implicit semi-Lagrangian model solving the compressible Euler equations in which the moist thermodynamics was formulated in terms of the Gibbs function for moist air. Whilst encouraging, this initial implementation suffered several limitations. First, because the formulation works in terms of the total Gibbs function for moist air, the water content must be partitioned into vapour and condensate whenever the Gibbs function is evaluated, complicating the calculation. More importantly, since the natural variables for the Gibbs function pressure p and temperature T are intensive variables, knowledge of p and T (and total specific humidity q) alone is insufficient to completely determine the equilibrium state, particularly the partition of water into its three phases, at the triple point. Thus, the implementation is restricted to two phases: vapour and liquid. Finally, the formulation of Thuburn (2017b) could not represent important nonequilibrium effects such as the delayed freezing of supercooled cloud droplets or the evaporation of rain in subsaturated air.

These limitations were overcome by Bowen and Thuburn (2022a, b). They formulated the thermodynamics in terms of the internal energy, whose natural variables are the extensive variables specific volume α and specific entropy η, avoiding the difficulty at the triple point. Moreover, by working with the individual internal energy potentials for dry air, water vapour, liquid water, and ice, rather than a combined thermodynamic potential for the air parcel, they were able to separate the calculation of the potentials from the calculation of the air parcel equilibrium state, simplifying the formulation. By expressing the evolution of a subset of variables in terms of thermodynamic forces and a resistivity matrix, they were able to account for departures from thermodynamic equilibrium, while seamlessly approaching the equilibrium case in the limit of zero resistivity.

The codes developed by Thuburn (2017b) and Bowen and Thuburn (2022a, b) were two-dimensional vertical slice models, and they were applied to simple buoyant bubble test problems. The apparent complexity of the thermodynamic potential approach, with accompanying concerns about its computational cost, might discourage model developers from pursuing the approach in three-dimensional models. A primary goal of the work described here is to demonstrate that the approach can be applied successfully, without excessive expense, in a three-dimensional model suitable for studying complex boundary layer and convective flows.

Traditionally, LES models for atmospheric applications are often based on relatively simple numerical methods supplemented by more or less sophisticated subgrid models (e.g. Siebesma et al., 2003). Typically, advection schemes are Eulerian and time stepping is explicit, so that stability requires the advective and diffusive Courant numbers to be less than some threshold value of order 1. Global weather and climate models, on the other hand, often use sophisticated (and relatively expensive) advection schemes that are stable for large advective Courant numbers (e.g. Temperton et al., 2001; Lin, 2004; Wood et al., 2014; Melvin et al., 2024), though they require the deformational Courant number to be bounded. Bartello and Thomas (1996) have argued that such large-time-step advection schemes are no longer cost-effective in flow regimes where the energy spectrum is shallower, the Lagrangian and Eulerian timescales become more comparable, and the deformational Courant number is much closer to the advective Courant number. Nevertheless, traditional LES codes are often run with time steps more than an order of magnitude smaller than could be used by a semi-implicit semi-Lagrangian scheme at the same resolution (e.g. Stevens et al., 2005, and compare section 5 below). This observation, combined with recent progress in improving the efficiency of conservative semi-implicit semi-Lagrangian solvers (Thuburn, 2024), encouraged the author to revisit the question by implementing a semi-implicit semi-Lagrangian LES model.

On a closely related point, global models with sophisticated advection schemes sometimes do not include a subgrid model to handle the turbulent downscale cascades of potential enstrophy and energy but rely instead on the dissipative nature of the advection scheme to play that role (e.g. Walters et al., 2017; ECMWF, 2023). In other words, they use a form of implicit large eddy simulation or ILES (e.g. Margolin et al., 2006; Grinstein et al., 2007). Although there have been some pioneering attempts to use or evaluate ILES for boundary layer and convective-scale atmospheric flows (Margolin et al., 1999; Brown et al., 2000; Smolarkiewicz and Prusa, 2002), and there is growing interest (e.g. Pressel et al., 2017; Souza et al., 2023), the approach is still far from mainstream, and there remain many open questions about its strengths and weaknesses and how they depend on details of the numerics. The need to answer these questions is becoming increasingly pressing as global prediction models begin to be used at the kilometre scale (Satoh et al., 2008; Stevens et al., 2019; Hohenegger et al., 2023; Tomassini et al., 2023). PTerodaC³TILES includes no subgrid model and uses the ILES approach. An important secondary goal for its development is to provide a tool to facilitate the study of the ILES properties of semi-implicit, semi-Lagrangian schemes for simulating convective-scale flows.

Finally, the author perceived a need for an LES tool that could be set up to run production science with minimal effort on a desktop machine or even a laptop. PTerodaC³TILES v1.0 comprises a single standalone fortran code and a namelist file. No external packages are needed other than a fortran compiler (with OpenMP shared memory parallel capability if desired). Just by selecting namelist options, the code can generate initial data and implement forcing terms for a number of standard test cases; new cases can easily be implemented by using the routines for existing cases as templates. To avoid large volumes of output, diagnostics are calculated online at run time. Many standard diagnostics are available via namelist switches; others can easily be implemented and output using existing routines as templates. Further details are given in the user manual, available from Zenodo; see the “Code availability” section.

Some aspects of the model formulation are novel or noteworthy; these are introduced briefly here and discussed in more detail in Sect. 2.

Unlike most atmospheric models, changes of phase of water and latent heat release are not treated as separate physics source terms but are fully integrated within the dynamical core's semi-implicit semi-Lagrangian time stepping (Thuburn, 2017b; Thuburn et al., 2022; Bowen and Thuburn, 2022a, b, Sect. 2.1 and 2.3 below). To maintain consistency of the thermodynamics, surface fluxes of water imply surface fluxes of mass. The same is true for the somewhat artificial sources of water specified in the domain interior for some test cases. Because of the model's Charney–Phillips vertical staggering, with total density stored at p levels and specific humidities stored at w levels, special care is needed to maintain that consistency (Sect. 2.7). Mass sources accompanying water sources are neglected in most atmospheric models.

In contrast to Thuburn (2017b) and Bowen and Thuburn (2022a, b), in PTerodaC³TILES conservative options are available for the advection of moisture and entropy variables. Although the numerical methods do not exactly conserve energy, energy conservation is significantly improved as a side effect of a conservative treatment of entropy and water (Thuburn, 2022).

Even with the use of cheap advection updates during the main solver iterations (Thuburn, 2024), advection remains one of the most expensive components of the model (along with the elliptic solver). Some modifications are made to the SLICE conservative semi-Lagrangian advection scheme (Zerroukat et al., 2009) to improve its efficiency, particularly to minimize the use of conditional code by avoiding searching. Geometrical calculations of coordinate line intersections to determine intermediate departure points are replaced by additional trajectory calculations, and information generated in remapping volume and mass is reused in remapping other fields (Sect. 2.6).

Linearization of the thermodynamics leads to an 11×14 linear subsystem at each model grid point that must be diagonalized in order to build the Helmholtz problem and to enable back substitution for the semi-implicit time stepping. To reduce what would otherwise potentially be a significant computational expense, the equations and unknowns are reordered to exploit the moderate sparsity of the thermodynamic subsystem, which is essentially the same at all grid points. The cost of the diagonalization is thereby reduced to about 30 % of the cost of a full Gaussian elimination for the case of equilibrium thermodynamics. The cost is further reduced by carrying out the diagonalization on all matrices in a grid column at the same time, with the inner loop over vertical levels, to improve vectorization (Sect. 2.10).

Some of the test cases implemented specify the use of Monin–Obukhov theory to compute surface momentum fluxes. Appendix A4 presents a slight reformulation of Monin–Obukhov theory that, with the aid of some curve fitting, enables the friction velocity U_* to be obtained without the need for an iterative calculation and guarantees the existence of a unique solution for U_* even when the validity of Monin–Obukhov theory breaks down.

Section 4 summarizes some of the ways in which the correctness of the formulation and implementation have been verified. Section 5 presents some sample results from standard LES test cases to demonstrate the performance of the model. The conclusions and areas where further work is needed are discussed in Sect. 6.

2 Model formulation

The formulation of PTerodaC³TILES is inspired by semi-implicit semi-Lagrangian dynamical cores such as ENDGame (Wood et al., 2014) and GungHo (Melvin et al., 2024) that iterate towards a (possibly off-centred) Crank–Nicolson time discretization along trajectories. Such a formulation has been found to be stable and robust at large time steps for synoptic-scale flow. A significant departure, though, is in the treatment of moist thermodynamics. First, to guarantee consistency, the thermodynamics is expressed in terms of thermodynamic potentials, in this case the internal energy. Second, processes such as phase changes and latent heat release are not treated as separate physics source terms but directly couple to the dynamics through the dependence of pressure and buoyancy on the thermodynamic state. The wide range of timescales associated with the thermodynamics, ranging from instantaneous for processes in equilibrium to many minutes for some nonequilibrium processes, is naturally handled by the semi-implicit time discretization.

A feature of the consistent treatment of thermodynamics is that any surface or interior source of moisture implies a corresponding source of mass (Sect. 2.7). Such a mass source is neglected in most atmospheric models.

PTerodaC³TILES version 1.0 includes only equilibrium thermodynamics, but almost all of the machinery is in place to deal with the nonequilibrium case. Also, this version includes no representation of microphysical processes; any condensate is simply carried along with the flow. The inclusion of a simple microphysics scheme is a priority for future work.

2.1 Continuous governing equations

The continuous equations to be solved are the same as those solved by Bowen and Thuburn (2022a), extended to three spatial dimensions, and with the inclusion of Coriolis terms and external source terms for mass, momentum, and water. See that article for a derivation and in-depth discussion. For the purpose of exposition it is convenient to split the governing equations into dynamics and thermodynamics.

2.1.1 Dynamics

The subset of governing equations involving material derivatives is the following:

$$\begin{array}{}\text{(1)}& {\displaystyle }& {\displaystyle \frac{\mathrm{1}}{\mathcal{V}}}{\displaystyle \frac{\mathrm{D}}{\mathrm{D}t}}\left(\mathcal{V}\mathit{\rho }\right)=S_{\mathit{\rho }},\text{(2)}& {\displaystyle }& {\displaystyle \frac{\mathrm{D}\mathit{u}}{\mathrm{D}t}}+\mathrm{2}\mathbf{\Omega }\times \mathit{u}+\mathit{\alpha }\mathrm{\nabla }p+\mathrm{\nabla }\mathrm{\Phi }={\mathit{S}}_{\mathit{u}},\text{(3)}& {\displaystyle }& {\displaystyle \frac{\mathrm{1}}{\mathcal{V}}}{\displaystyle \frac{\mathrm{D}}{\mathrm{D}t}}\left(\mathcal{V}\mathit{\rho }q\right)=\mathit{\rho }{S}_q,\text{(4)}& {\displaystyle }& {\displaystyle \frac{\mathrm{1}}{\mathcal{V}}}{\displaystyle \frac{\mathrm{D}}{\mathrm{D}t}}\left(\mathcal{V}\mathit{\rho }\mathit{\eta }\right)=\mathit{\rho }{\mathit{J}}^{\mathrm{T}}\mathit{P}+\mathit{\rho }{S}_{\mathit{\eta }},\text{(5)}& {\displaystyle }& {\displaystyle \frac{\mathrm{D}\mathit{X}}{\mathrm{D}t}}=\mathit{J}+{\mathit{S}}_{\mathit{X}}.\end{array}$$

Here, ρ is the total fluid density, $\mathit{u}=(u,v,w)$ is the fluid (barycentric) velocity, q is the total specific humidity, η is the total specific entropy, $\mathit{\alpha }=\mathrm{1}/\mathit{\rho }$ is the total specific volume, p is the pressure, and Ω is the rotation vector of the frame of reference. Φ=gz is the geopotential with g the gravitational acceleration. S and S terms indicate external sources.

The vector $\mathit{X}=(q^{\mathrm{l}},q^{\mathrm{f}},q^{\mathrm{v}}{\mathit{\eta }}^{\mathrm{v}},q^{\mathrm{l}}{\mathit{\eta }}^{\mathrm{l}},q^{\mathrm{f}}{\mathit{\eta }}^{\mathrm{f}}{)}^{\mathrm{T}}$ (superscript T meaning transpose) encodes the additional thermodynamic information that needs to be predicted in the nonequilibrium case. Superscripts d, v, l, and f indicate dry air, water vapour, liquid water, and frozen water, respectively. J is a vector of thermodynamic fluxes, with P the corresponding thermodynamic forces, defined below. The expression J^TP in Eq. (4) gives the entropy source per unit mass due to nonequilibrium processes¹.

$\mathrm{D}/\mathrm{D}t$ is the material derivative. In Eqs. (1), (3), and (4), 𝒱 is the material volume element. These equations are written in a form that lends itself to numerical solution using a conservative semi-Lagrangian scheme. The material derivative in Eq. (5) is written in a form that anticipates discretization using an interpolating semi-Lagrangian scheme, on the assumption that it will be sufficient to advect the total specific humidity and entropy conservatively, with a cheaper non-conserving scheme for the components of X. Nevertheless, an option for conservative advection is available. In the equilibrium case, however, J is not needed (see Eq. 13 below), so advection of X becomes superfluous and is automatically switched off.

The user can also choose to include additional tracers stored either at p levels or at w levels (see Sect.  2.2). Since these do not feed back on the dynamics or thermodynamics, we largely omit them from further discussion.

2.1.2 Thermodynamics

The diagnostic equations describing the thermodynamics are as follows:

$$\begin{array}{}\text{(6)}& {\displaystyle }& {\displaystyle }{q}^{\mathrm{v}}{\mathit{\alpha }}^{\mathrm{v}}-q^{\mathrm{d}}{\mathit{\alpha }}^{\mathrm{d}}=\mathrm{0},\text{(7)}& {\displaystyle }& {\displaystyle }\left(\mathrm{1}-{\mathit{\delta }}^{\mathrm{l}}\right){\mathit{\lambda }}^{\mathrm{l}}-{\mathit{\delta }}^{\mathrm{l}}{q}^{\mathrm{l}}=\mathrm{0},\text{(8)}& {\displaystyle }& {\displaystyle }\left(\mathrm{1}-{\mathit{\delta }}^{\mathrm{f}}\right){\mathit{\lambda }}^{\mathrm{f}}-{\mathit{\delta }}^{\mathrm{f}}{q}^{\mathrm{f}}=\mathrm{0},\text{(9)}& {\displaystyle }& {\displaystyle }{q}^{\mathrm{v}}+q^{\mathrm{l}}+q^{\mathrm{f}}-q=\mathrm{0},\text{(10)}& {\displaystyle }& {\displaystyle }{q}^{\mathrm{v}}{\mathit{\alpha }}^{\mathrm{v}}+q^{\mathrm{l}}{\mathit{\alpha }}^{\mathrm{l}}+q^{\mathrm{f}}{\mathit{\alpha }}^{\mathrm{f}}-\mathit{\alpha }=\mathrm{0},\text{(11)}& {\displaystyle }& {\displaystyle }{q}^{\mathrm{d}}{\mathit{\eta }}^{\mathrm{d}}+q^{\mathrm{v}}{\mathit{\eta }}^{\mathrm{v}}+q^{\mathrm{l}}{\mathit{\eta }}^{\mathrm{l}}+q^{\mathrm{f}}{\mathit{\eta }}^{\mathrm{f}}-\mathit{\eta }=\mathrm{0},\text{(12)}& {\displaystyle }& {\displaystyle }p+e_{\mathit{\alpha }}^{\mathrm{d}}+e_{\mathit{\alpha }}^{\mathrm{v}}=\mathrm{0},\text{(13)}& {\displaystyle }& {\displaystyle }\mathbf{R}\mathit{J}-(\mathit{P}-\mathit{C})\equiv -{\mathit{R}}_{\mathrm{PE}}=\mathrm{0}.\end{array}$$

Where needed, the mass fraction of dry air is given by $q^{\mathrm{d}}=\mathrm{1}-q$.

Equation (6) states that the water vapour and dry air occupy the same volume within an air parcel. Equations (9)–(11) state that the total specific humidity equals q, the total specific volume equals α, and the total specific entropy equals η.

Equation (12) expresses Dalton's law of partial pressures and is key for coupling the thermodynamics to the dynamics. Here, e^d(α^d,η^d), e^v(α^v,η^v), e^l(η^l), and e^f(η^f) are the internal energy potentials for dry air, water vapour, liquid water, and frozen water, respectively, expressed as functions of their natural variables (Appendix A2). As in Bowen and Thuburn (2022a, b), condensate is assumed to be incompressible, with α^l and α^f specified constants, so we suppress the dependence of e^l on α^l and of e^f on α^f, and the pressure within any condensate (for example, as needed to compute the Gibbs function) is equal to the total pressure of the surrounding gas. Subscripts α and η on the species internal energies indicate partial derivatives with respect to the respective natural variables.

Equation (13) is a set of phenomenological equations relating thermodynamic fluxes J to the thermodynamic forces P that push the system towards equilibrium (de Groot and Mazur, 1984). R_PE is a residual that should be driven towards zero by the iterative solver. In the present system P is given by

$$\begin{array}{}\text{(14)}& \begin{array}{rl}\mathit{P}=& {\mathrm{\nabla }}_{\mathit{X}}\mathit{\eta }={\displaystyle \frac{\mathrm{1}}{T^{\mathrm{d}}}}\left(g^{\mathrm{v}}-g^{\mathrm{l}},g^{\mathrm{v}}-g^{\mathrm{f}},T^{\mathrm{d}}-T^{\mathrm{v}},T^{\mathrm{d}}-T^{\mathrm{l}},\right.\\ & {\left.T^{\mathrm{d}}-T^{\mathrm{f}}\right)}^{\mathrm{T}},\end{array}\end{array}$$

where the gradient ∇_Xη is taken at constant mass and energy, Tⁱ is the temperature of species i, and gⁱ is the Gibbs function for species i, $i\in \mathit{\{}\mathrm{d},\mathrm{v},\mathrm{l},\mathrm{f}\mathit{\}}$ (Appendix A2). R is a symmetric and positive semi-definite resistivity matrix controlling the rate at which an air parcel approaches thermodynamic equilibrium. $\mathit{C}={\mathit{\eta }}_{\mathrm{ref}}({\mathit{\delta }}^{\mathrm{l}}{\mathit{\lambda }}^{\mathrm{l}},{\mathit{\delta }}^{\mathrm{f}}{\mathit{\lambda }}^{\mathrm{f}},\mathrm{0},\mathrm{0},\mathrm{0}{)}^{\mathrm{T}}$ is a vector of switched Lagrange multipliers, with ${\mathit{\delta }}^{\mathrm{l}},{\mathit{\delta }}^{\mathrm{f}}\in \mathit{\{}\mathrm{0},\mathrm{1}\mathit{\}}$, used in enforcing the constraints that q^l and q^f must be non-negative. For example, if the thermodynamic forces imply evaporation of liquid but q^l is already zero, then the liquid constraint is switched on (δ^l=1) and the solution for λ^l balances the relevant component of P. Finally, Eqs. (7) and (8) are a convenient way to express the complementarity conditions (e.g. Nocedal and Wright, 2006) that either q^l or λ^l must be zero and either q^f or λ^f must be zero. η_ref is an arbitrary reference value, here equal to 1000 J kg⁻¹ K⁻¹, introduced to ensure that Eqs. (7) and (8) are dimensionally correct.

Bowen and Thuburn (2022b) show how R can be related to the thermal conductivity of air and the molecular diffusivity of water vapour in air for cloud droplets of a given radius. In PTerodaC³TILES v1.0 the resistivity R is set to zero, imposing local thermodynamic equilibrium. Nevertheless, almost all of the machinery is in place to handle the nonequilibrium case.

Note that PTerodaC³TILES has no explicit representation of subgrid variability in temperature or humidity, hence in condensate. Each grid cell is either entirely saturated or entirely unsaturated. In other words, it is an all-or-nothing representation of saturation and condensate.

2.2 Domain, grid, and discretization

The domain of PTerodaC³TILES v1.0 is rectangular and doubly periodic in the horizontal, with flat rigid boundaries at the bottom and top. Cartesian coordinates (x, y, z) are aligned with the domain. When a full three-dimensional Coriolis force is used, then the x and y directions are assumed to be east and north, respectively, but otherwise the model is agnostic about which direction is north.

A C-grid staggering (Arakawa and Lamb, 1977) is used in the horizontal, and a Charney–Phillips grid staggering (Charney and Phillips, 1953) is used in the vertical (Fig. 1). The choice of a Charney–Phillips vertical grid is unusual for LES models but is more common in global models. It avoids computational modes – oscillatory vertical profiles of thermodynamic variables that spuriously satisfy discrete hydrostatic balance – and gives more accurate coupling between vertical velocity and buoyancy on small vertical scales but makes it difficult to formulate an exactly energy-conserving scheme.

Figure 1 Schematic showing the placement of model variables on the vertically staggered Charney–Phillips grid. Superscript i can stand for any of d, v, l, and f.

In the following, integer indices are used for p points, with an offset of $\mathrm{1}/\mathrm{2}$ in the relevant direction for velocity points. Vertical indices 1 and N_z correspond to the lowermost and uppermost p levels; vertical indices $\mathrm{1}/\mathrm{2}$ and $N_z+\mathrm{1}/\mathrm{2}$ correspond to the w levels at the bottom and top model boundaries. Horizontal grid spacings Δx and Δy are uniform but need not be equal to each other. The vertical grid spacing Δz may be uniform or stretched. See Appendix A1 for the grid specification in the stretched case. All of the results shown in Sects. 4 and 5 use uniform Δz.

Simple finite-difference or finite-volume approximations are used for gradient and divergence operators. The components of the pressure gradient (and similarly the geopotential gradient) are needed at velocity points:

$$\begin{array}{}\text{(15)}& \begin{array}{rl}& {\left.{\displaystyle \frac{\partial p}{\partial x}}\right|}_{i+\mathrm{1}/\mathrm{2}jk}={\displaystyle \frac{p_{i+\mathrm{1}jk}-p_{ijk}}{\mathrm{\Delta }x}},\\ & {\left.{\displaystyle \frac{\partial p}{\partial y}}\right|}_{ij+\mathrm{1}/\mathrm{2}k}={\displaystyle \frac{p_{ij+\mathrm{1}k}-p_{ijk}}{\mathrm{\Delta }y}},\\ & {\left.{\displaystyle \frac{\partial p}{\partial z}}\right|}_{ijk+\mathrm{1}/\mathrm{2}}={\displaystyle \frac{p_{ijk+\mathrm{1}}-p_{ijk}}{\mathrm{\Delta }{z}_{k+\mathrm{1}/\mathrm{2}}}},\end{array}\end{array}$$

where $\mathrm{\Delta }{z}_{k+\mathrm{1}/\mathrm{2}}=z_{k+\mathrm{1}}-z_k$ and with the obvious modifications to allow for periodic boundary conditions. The divergence of the velocity is needed at p points and is given by

$$\begin{array}{}\text{(16)}& \begin{array}{rl}{\left.\mathrm{\nabla }\cdot \mathit{u}\right|}_{ijk}& ={\displaystyle \frac{u_{i+\mathrm{1}/\mathrm{2}jk}-u_{i-\mathrm{1}/\mathrm{2}jk}}{\mathrm{\Delta }x}}+{\displaystyle \frac{v_{ij+\mathrm{1}/\mathrm{2}k}-v_{ij-\mathrm{1}/\mathrm{2}k}}{\mathrm{\Delta }y}}\\ & +{\displaystyle \frac{w_{ijk+\mathrm{1}/\mathrm{2}}-w_{ijk-\mathrm{1}/\mathrm{2}}}{\mathrm{\Delta }{z}_k}},\end{array}\end{array}$$

where $\mathrm{\Delta }{z}_k=z_{k+\mathrm{1}/\mathrm{2}}-z_{k-\mathrm{1}/\mathrm{2}}$, with an analogous expression for the divergence of mass flux increments.

Averaging between velocity points and p points is needed at various places in the discretization. Horizontal averaging uses a simple two-point average with weights $\mathrm{1}/\mathrm{2}$, indicated by ${\stackrel{\mathrm{‾}}{(.)}}^x$ or ${\stackrel{\mathrm{‾}}{(.)}}^y$. For example, the values of α used in the horizontal components of Eq. (2) are given by $\mathrm{1}/{\stackrel{\mathrm{‾}}{\mathit{\rho }}}^x$ or $\mathrm{1}/{\stackrel{\mathrm{‾}}{\mathit{\rho }}}^y$.

Because the vertical grid may be stretched, four different vertical averaging operators are possible:

• ${\stackrel{\mathrm{‾}}{(.)}}^w$, linear interpolation from p levels to w levels;

• ${\stackrel{\mathrm{‾}}{(.)}}^r$, piecewise constant conservative remapping from p levels to w levels;

• ${\stackrel{\mathrm{‾}}{(.)}}^p$, linear interpolation from w levels to p levels; and

• ${\stackrel{\mathrm{‾}}{(.)}}^s$, piecewise constant conservative remapping from w levels to p levels.

Complete expressions for these four operators along with some useful conservation and discrete product rule properties are given in Appendix B of Thuburn et al. (2022). For example, the averaging of velocity components to compute the Coriolis terms (Eq. 2) and departure points (Sect. 2.6) uses a combination of ${\stackrel{\mathrm{‾}}{(.)}}^x$, ${\stackrel{\mathrm{‾}}{(.)}}^y$, ${\stackrel{\mathrm{‾}}{(.)}}^p$, and ${\stackrel{\mathrm{‾}}{(.)}}^w$. When a value for p is needed at the bottom and top boundaries, a variant of the ${\stackrel{\mathrm{‾}}{(.)}}^w$ operator is employed that uses linear extrapolation rather than the default constant extrapolation.

In order to ensure consistent and conservative transport of water and entropy, a key aspect of the model formulation is that a mass budget for a dual w-level density ${\stackrel{\mathrm{‾}}{\mathit{\rho }}}^r$ that is consistent with the p-level ρ budget should be satisfied (Konor and Arakawa, 2000; Thuburn, 2022; Bendall et al., 2023). Achieving this consistency requires some care in the advection of w-level scalars (Sect. 2.6), both for SLICE and for the cheap advection updates, and also in the discrete formulation of surface sources and interior sources (Sect. 2.7).

In the following sections, for clarity, details of the spatial discretization are suppressed except to indicate where different vertical averaging operators are used.

2.3 Overview of information flow

Before delving further into details, it is useful to take an overview of the flow of information between the dynamics and thermodynamics parts of the model formulation (Fig. 2) to help clarify the organization of the linearization and the derivation of the Helmholtz problem in Sect. 2.8–2.10.

Figure 2 Schematic showing the flow of information between the dynamics, thermodynamics at p levels, and thermodynamics at w levels for the case of equilibrium thermodynamics. (For the nonequilibrium case the w-level thermodynamics retains some memory of qⁱ and ηⁱ.) Superscript i stands for any of d, v, l, or f. An overbar indicates a vertical averaging operation from p levels to w levels or vice versa.

For the continuous equations, the role of the thermodynamics is to return the pressure p, given the density ρ, total water q, and total entropy η. To do so, the thermodynamics must determine how q and η are partitioned among the different components and phases q^v, q^l, q^f, η^d, η^v, η^l, and η^f.

The situation is more complicated when discretized with a Charney–Phillips vertical grid staggering, since p and ρ are stored at p levels while q and η are stored at w levels. Moreover, in order to obtain optimal coupling between vertical velocity and buoyancy, the w-level specific volume that appears in the vertical pressure gradient term in Eq. (2) must be calculated not simply from a vertically averaged density $\mathrm{1}/{\stackrel{\mathrm{‾}}{\mathit{\rho }}}^r$, but from the w-level η and q and a vertically averaged pressure ${\stackrel{\mathrm{‾}}{p}}^w$ (Thuburn, 2017a). It is convenient to express this requirement through an additional equation,

$$\begin{array}{}\text{(17)}& {\stackrel{\mathrm{‾}}{p}}^w-p^{\left(w\right)}\equiv -R_{\stackrel{\mathrm{‾}}{p}}=\mathrm{0},\end{array}$$

where p^(w) is the pressure appearing in Eq. (12) on w levels, and $R_{\stackrel{\mathrm{‾}}{p}}$ is a residual that should be driven to zero by the iterative solver. Thus, on w levels a full set of thermodynamic equations (Eqs. 6–13) is solved, while on p levels the following subset of thermodynamic equations is solved:

$$\begin{array}{}\text{(18)}& {\displaystyle }& {\displaystyle }{\stackrel{\mathrm{‾}}{q^{\mathrm{v}}}}^p{\mathit{\alpha }}^{\mathrm{v}\left(p\right)}+{\stackrel{\mathrm{‾}}{q^{\mathrm{l}}}}^p{\mathit{\alpha }}^{\mathrm{l}}+{\stackrel{\mathrm{‾}}{q^{\mathrm{f}}}}^p{\mathit{\alpha }}^{\mathrm{f}}-{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}=\mathrm{0},\text{(19)}& {\displaystyle }& {\displaystyle }{\stackrel{\mathrm{‾}}{q^{\mathrm{v}}}}^p{\mathit{\alpha }}^{\mathrm{v}\left(p\right)}-{\stackrel{\mathrm{‾}}{q^{\mathrm{d}}}}^p{\mathit{\alpha }}^{\mathrm{d}\left(p\right)}=\mathrm{0},\text{(20)}& {\displaystyle }& {\displaystyle }p+e_{\mathit{\alpha }}^{\mathrm{d}}+e_{\mathit{\alpha }}^{\mathrm{v}}=\mathrm{0}.\end{array}$$

Equation (18) determines the p-level specific volume of water vapour α^v(p), and then Eq. (19) determines the p-level specific volume of dry air α^d(p). Finally, the pressure is computed from the dry air and water vapour internal energy potentials in terms of their p-level natural variables α^d(p), ${\stackrel{\mathrm{‾}}{{\mathit{\eta }}^{\mathrm{d}}}}^p$, α^v(p), and ${\stackrel{\mathrm{‾}}{{\mathit{\eta }}^{\mathrm{v}}}}^p$. Note that the specific volumes α^d(p) and α^v(p) are calculated directly from ρ on p levels rather than using vertical averages of the w-level values ${\stackrel{\mathrm{‾}}{{\mathit{\alpha }}^{\mathrm{d}}}}^p$, ${\stackrel{\mathrm{‾}}{{\mathit{\alpha }}^{\mathrm{v}}}}^p$, ensuring that p responds correctly and locally to changes in ρ.

2.4 Semi-implicit semi-Lagrangian scheme

Splitting the momentum equation into its horizontal and vertical components with $\mathit{v}=(u,v,\mathrm{0})$ being the horizontal velocity, a semi-implicit semi-Lagrangian discretization of Eqs. (1)–(5) is

$$\begin{array}{}\text{(21)}& {\left[\mathit{\rho }-a\mathrm{\Delta }t{S}_{\mathit{\rho }}\right]}^{n+\mathrm{1}}-{\left[\mathit{\rho }+b\mathrm{\Delta }t{S}_{\mathit{\rho }}\right]}_{\mathrm{T}}^n\equiv -R_{\mathit{\rho }}=\mathrm{0},\end{array}$$

$$\begin{array}{}\text{(22)}& {\displaystyle }\begin{array}{rl}\left[\mathit{v}\right.& {\left.+a\mathrm{\Delta }t\left((\mathrm{2}\mathbf{\Omega }\times \mathit{u}{)}_{\mathbf{H}}+\mathit{\alpha }{\mathrm{\nabla }}_{\mathbf{H}}p-{\mathit{S}}_{\mathit{v}}\right)\right]}^{n+\mathrm{1}}\\ & -{\left[\mathit{v}-b\mathrm{\Delta }t\left((\mathrm{2}\mathbf{\Omega }\times \mathit{u}{)}_{\mathbf{H}}+\mathit{\alpha }{\mathrm{\nabla }}_{\mathbf{H}}p-{\mathit{S}}_{\mathit{v}}\right)\right]}_{\mathrm{D}}^n\\ & \equiv -{\mathit{R}}_{\mathit{v}}=\mathbf{0},\end{array}\text{(23)}& {\displaystyle }\begin{array}{rl}\left[w\right.& {\left.+a\mathrm{\Delta }t\left((\mathrm{2}\mathbf{\Omega }\times \mathit{u}{)}_{\mathbf{V}}+\mathit{\alpha }{\displaystyle \frac{\partial p}{\partial z}}+{\displaystyle \frac{\partial \mathrm{\Phi }}{\partial z}}-S_w\right)\right]}^{n+\mathrm{1}}\\ & -{\left[w-b\mathrm{\Delta }t\left((\mathrm{2}\mathbf{\Omega }\times \mathit{u}{)}_{\mathbf{V}}+\mathit{\alpha }{\displaystyle \frac{\partial p}{\partial z}}+{\displaystyle \frac{\partial \mathrm{\Phi }}{\partial z}}-S_w\right)\right]}_{\mathrm{D}}^n\\ & \equiv -R_w=\mathrm{0},\end{array}\text{(24)}& {\displaystyle }\begin{array}{rl}{\left[{\stackrel{\mathrm{‾}}{\mathit{\rho }}}^r\left(\mathit{\eta }-a\mathrm{\Delta }t{S}_{\mathit{\eta }}\right)\right]}^{n+\mathrm{1}}& -{\left[{\stackrel{\mathrm{‾}}{\mathit{\rho }}}^r\left(\mathit{\eta }+b\mathrm{\Delta }t{S}_{\mathit{\eta }}\right)\right]}_{\mathrm{T}}^n\\ & -\mathrm{\Delta }t{\stackrel{\mathrm{‾}}{\mathit{\rho }}}^r{\mathit{J}}^{\mathrm{T}}{\mathit{P}}^{n+\mathrm{1}}\equiv -R_{\mathit{\rho }\mathit{\eta }}=\mathrm{0},\end{array}\text{(25)}& {\displaystyle }\begin{array}{rl}{\left[{\stackrel{\mathrm{‾}}{\mathit{\rho }}}^r\left(q-a\mathrm{\Delta }t{S}_q\right)\right]}^{n+\mathrm{1}}& -{\left[{\stackrel{\mathrm{‾}}{\mathit{\rho }}}^r\left(q+b\mathrm{\Delta }t{S}_q\right)\right]}_{\mathrm{T}}^n\\ & \equiv -R_{\mathit{\rho }q}=\mathrm{0},\end{array}\text{(26)}& {\displaystyle }{\left[\mathit{X}-a\mathrm{\Delta }t{\mathit{S}}_{\mathit{X}}\right]}^{n+\mathrm{1}}-{\left[\mathit{X}+b\mathrm{\Delta }t{\mathit{S}}_{\mathit{X}}\right]}_{\mathrm{D}}^n-\mathrm{\Delta }t\mathit{J}=\mathbf{0}.\end{array}$$

Here, superscripts n and n+1 are time step indices. Subscripts H and V indicate horizontal and vertical components, respectively. Subscript D indicates a quantity interpolated or remapped to a semi-Lagrangian trajectory departure point or cell, while subscript T indicates a conservatively transported quantity defined by ${\mathcal{V}}^{n+\mathrm{1}}{\mathit{\psi }}_{\mathrm{T}}=[\mathcal{V}\mathit{\psi }{]}_{\mathrm{D}}^n$.

Δt is the time step, and a and $b=\mathrm{1}-a$ are off-centring parameters for the dynamics. To damp acoustic waves that might be generated by initial perturbations imposed to trigger turbulence or by the switching on of forcing terms at the initial time, a is smoothly adjusted from 1 to a user-specified value over the first 900 s of a model run. A value of a=0.51 (after the initial adjustment) is used for the results shown in Sect. 5.

The terms R_ρ, R_v, R_w, R_ρη, and R_ρq represent residuals in their respective equations. In the target discretization these residuals should be zero. The iterative solver described in the following sections attempts to drive those residuals to zero. No residual appears in Eq. (26) because that equation is used to diagnose the fluxes J so it is always satisfied exactly.

The subsystem of thermodynamics equations (Eqs. 6–13) is solved at step n+1, and it is Pⁿ⁺¹ that appears in Eq. (24), effectively giving a backward Euler treatment of the thermodynamics. The thermodynamic processes of interest typically involve relaxation towards equilibrium rather than oscillation about equilibrium. A backward Euler treatment should be sufficiently accurate for nonequilibrium processes whose timescale is much longer than Δt, such as evaporation of falling rain. For processes whose timescale is shorter than about $\mathrm{\Delta }t/\mathrm{2}$, a backward Euler step is preferable to a Crank–Nicolson step since Crank–Nicolson can overshoot the equilibrium solution. If the form of the resistivity matrix implies that any part of the system is in equilibrium (e.g. T^v=T^d), then a backward Euler step is essential, since the equilibrium must be imposed at step n+1 and not as a time average.

2.5 Time-stepping algorithm

The semi-implicit semi-Lagrangian scheme described in Sect. 2.4 is both nonlinear in the unknown step n+1 values and nonlocal because unknown values at neighbouring grid points are coupled. The equations are solved using an iterative quasi-Newton algorithm (Algorithm 1). By elimination of unknowns, the linear system for the Newton update is reduced to a more or less standard Helmholtz problem, which is solved using a multigrid method (Sect. 2.11). To avoid the expense of computing conservative semi-Lagrangian transport multiple times per step, a single full advection calculation is made once at the start of the time step, and relatively cheap updates to the transport are made at each solver iteration (Zerroukat and Allen, 2020). These cheap transport updates use simple upwind or centred schemes and are made in such a way that the transport of scalars remains conservative and bounded (assuming conservative and bounded options have been chosen by the user) and consistent with the transport of mass at every solver iteration (Thuburn, 2024).

Algorithm 1 Computations performed to take one model time step.

Compute time step n terms

Initialize step n+1 state variables to step n state

Full advection calculation: compute [.]_D and [.]_T terms

for ℓ=1 to N_ℓ do

Compute step n+1 terms based on iteration ℓ−1 values

Compute residuals

Build and solve the Helmholtz problem

Back-substitute, updating all transport terms and state variables

end for

The number of solver iterations takes a default value N_ℓ=3. This was found to be sufficient for all test cases simulated except DYCOMS, for which N_ℓ=4 was needed. This exceptional case is discussed in Sects. 4 and 5.

When cheap transport updates are employed, the resulting time integration scheme at solver convergence is not quite as written in Sect. 2.4, with semi-Lagrangian advection by the trajectory-average velocity. Rather, the net transport results from semi-Lagrangian advection using the first-guess trajectory-average velocity followed by a sequence of small corrections. Nevertheless, the end result is very close to the target discretization and appears to work well in practice.

An attractive aspect of the use of the consistent and conservative cheap transport updates is that, after the first solver iteration, the residuals in Eqs. (21), (24), and (25) and the equations for any advected tracers are very small and result only from changes in source terms between one iteration and the next (Thuburn, 2024). It is likely that this helps solver convergence (Sect. 4.4).

2.6 Advection

All advection terms are computed using semi-Lagrangian schemes. Different advection options are available depending on the advected variable and its location on the grid. All momentum equation components are advected with an interpolating (non-conservative) semi-Lagrangian scheme. Terms in the total fluid density equation and the mass fractions of any p-level tracers are advected using the mass-conserving SLICE scheme (Zerroukat et al., 2009). Advected terms in the water and entropy equations (Eqs. 24–26), as well as any w-level tracers, may be advected with either an interpolating semi-Lagrangian scheme or with SLICE.

Both the interpolating semi-Lagrangian scheme and SLICE use the cascade idea (Purser and Leslie, 1991) to replace a three-dimensional interpolation or remapping by a sequence of one-dimensional interpolations or remappings (Fig. 3).

Figure 3 Departure points and cascade remapping in two dimensions. The red line D_f to A_f is an example full trajectory to a v-point arrival point A_f from the corresponding departure point D_f. The x and y coordinates of A_f are known, and the x and y coordinates of D_f are to be calculated. The blue line D_i to A_i is an example intermediate trajectory to A_i from D_i. The x coordinate of A_i and the y coordinate of D_i are known, and the y coordinate of A_i and the x coordinate of D_i are to be calculated. Once the intermediate and full departure points are found, fields are conservatively remapped in the x direction, from the model grid cells (straight black lines) to intermediate cells (dashed blue lines), and then in the s direction to the departure cells (red curves). The remapped field is then transported from its departure cells (e.g. pink cell) to the corresponding arrival cells (e.g. grey cell).

The interpolating semi-Lagrangian scheme is based on cubic Lagrange interpolation. For the vertical interpolation, modifications are needed near the upper and lower boundaries. Between the uppermost pair of data points and between the lowermost pair of data points a modified cubic interpolation is used that uses only the two nearest data points and is very close to a linear interpolation. For terms in the u and v equations, when extrapolation above the uppermost data point is needed constant extrapolation is used. When extrapolation below the lowest data point is needed, the scheme uses either constant extrapolation or linear interpolation between the data value at level 1 and zero at the surface, according to whether the user selects a freeslip or noslip$\bigsqcup \bigsqcup$ boundary condition for the advection.

For the SLICE scheme, for each advected field the user may choose between piecewise constant, piecewise parabolic (Colella and Woodward, 1984), and parabolic spline (Zerroukat et al., 2007) remapping schemes. Experimentation suggests that, for most purposes, piecewise constant remapping is adequate for cell volume, divergence, and density. The results shown in Sect. 5 use this option, with parabolic spline method remapping for entropy and water. For each advected field the user has the option to use a simple limiter ensuring boundedness of the advected field. The limiter is redundant in the case of SLICE advection with piecewise constant remapping. In Sect. 5 the limiter is used for advection of entropy and water but not for advection of velocity components.

Some modifications to previous implementations of cascade advection schemes are made to reduce expensive searching and conditional code and to maximize reuse of information. First, in addition to the trajectory departure points, three-dimensional cascade interpolation or remapping requires two sets of intermediate departure points. Rather than construct these intermediate departure points by computing intersections between the arrival coordinate system defined by the model grid and the departure or Lagrangian coordinate system, as in previous work (e.g. Purser and Leslie, 1991; Nair et al., 2002; Zerroukat et al., 2002), here they are computed by separate trajectory calculations. The idea is illustrated in Fig. 3 for the two-dimensional case. A first guess followed by a single fixed point iteration is used for all trajectory calculations.

Second, we want to remap the density field using a volume coordinate to ensure consistency with the trajectory-average divergence (see below), and we want to remap water, entropy, and tracers in a mass coordinate to ensure that their advection is conservative, bounded if desired, and consistent with the density advection. However, neither the volume coordinate nor the mass coordinate is a simple function of cell index, so it might appear, at first glance, that expensive searching is needed to determine the necessary origin grid indices in these coordinates. This difficulty is circumvented with the aid of the following insight. To compute a one-dimensional remapping of a field f from an origin grid to a destination grid, the information needed comprises the origin grid cell average values of f, the origin grid coordinate intervals Δs_k, the origin grid indices k_i corresponding to destination grid cell edges i, and the cell fractions ξ_i in the s coordinate (Fig. 4). Fields are then remapped in the following sequence.

1. Remap cell volume in a geometrical coordinate x, y, or z; the origin grid indices k_i and cell fractions ξ_i are easily determined.

2. Remap density in a volume coordinate; the origin grid indices k_i are unchanged while the cell fractions ξ_i are obtained as a by-product of the volume remapping. For example, if Fig. 4 represents the remapping of volume in a geometrical coordinate, then the volume-coordinate cell fraction needed for the density remapping is given by the ratio of the dark shaded area to the full shaded area.

3. Remap p-level tracer mass fractions in a mass coordinate; the origin grid indices k_i are again unchanged while the cell fractions ξ_i are obtained as a by-product of the density remapping.

4. Remap w-level scalars in a mass coordinate. The origin grid indices k_i and cell fractions ξ_i must correspond to a mass coordinate based on ${\stackrel{\mathrm{‾}}{\mathit{\rho }}}^r$; they can be computed from the mass coordinate k_i and ξ_i at the p levels immediately below and above the w level in question.

Figure 4 Information used for one-dimensional conservative remapping of a field between an origin grid (cell edges indicated in blue) and a destination grid (cell edges indicated in red). The remapping coordinate s may be a geometrical coordinate (x, y, or z), volume, or mass. Δs_k is an origin grid cell size in the s coordinate. k_i is the origin grid cell index in which the destination grid edge i sits. ξ_i is the fraction of the origin grid cell k_i, measured in the s coordinate, to the left of destination edge i.

Thuburn et al. (2010) found that the accuracy of a semi-implicit semi-Lagrangian scheme with conservative semi-Lagrangian advection of density relies on the semi-Lagrangian departure volumes (departure areas in the shallow water context) being consistent with the trajectory-average divergence. Subsequent work showed that this divergence consistency criterion is crucial for stability too, unless a strong off-centring is employed. However, it appears very difficult to enforce this criterion directly within the trajectory calculations. Thuburn et al. (2010) solved this problem by conservatively advecting the divergence field, hence constructing the required departure cell areas, and using an area coordinate for the final SLICE remapping sweep in the advection of mass. Here a slightly different approach is taken to enable the remapping sequence discussed above. First, the divergence field is conservatively advected, and the information is used to compute the required departure cell volumes. These required departure cell volumes are compared with the actual departure cell volumes returned by SLICE, and the difference is used to compute a (small) divergent velocity increment sufficient to correct the discrepancy. The semi-Lagrangian advection of all other advected fields is then carried out, following which the machinery for making cheap advection updates is used to update the advection of all fields using this velocity increment and so satisfy divergence consistency.

One of the benefits of the Charney–Phillips vertical grid is that the colocation of w and η permits a tight coupling between vertical motion and buoyancy. However, when the entropy transport is computed in flux form, that tight coupling is lost unless the horizontal fluxes are corrected to account for the vertical gradients of η, $F_x={\stackrel{\mathrm{‾}}{\mathit{\rho }}}^{x}u$, and $F_y={\stackrel{\mathrm{‾}}{\mathit{\rho }}}^{y}v$ (Thuburn, 2022). Thuburn (2022) showed how the horizontal remapping of w-level fields in SLICE could be modified to implement the required correction. Here, since the correction is generally small, a much simpler approach is taken. Before carrying out the main SLICE advection of entropy, a small transport correction is made to the entropy using the horizontal fluxes,

$$\begin{array}{}\text{(27)}& \begin{array}{rl}& F_x^{\mathit{\eta }}=a^w{b}^w\mathrm{\Delta }{z}^{\mathrm{2}}{\displaystyle \frac{\partial F_x}{\partial z}}{\stackrel{\mathrm{‾}}{{\displaystyle \frac{\partial {\stackrel{\mathrm{‾}}{\mathit{\eta }}}^x}{\partial z}}}}^w,\\ & F_y^{\mathit{\eta }}=a^w{b}^w\mathrm{\Delta }{z}^{\mathrm{2}}{\displaystyle \frac{\partial F_y}{\partial z}}{\stackrel{\mathrm{‾}}{{\displaystyle \frac{\partial {\stackrel{\mathrm{‾}}{\mathit{\eta }}}^y}{\partial z}}}}^w,\end{array}\end{array}$$

where a^w and b^w are the coefficients associated with the ${\stackrel{\mathrm{‾}}{(.)}}^w$ operator. All quantities are evaluated at step n, and the fluxes are applied over the time step Δt. Subscripts $k+\mathrm{1}/\mathrm{2}$ indicating the vertical level have been omitted for clarity. Equation (27) is a slight generalization of Eq. (7) of Thuburn (2022) to allow for a vertically stretched grid. An analogous correction is applied to the transport of all conservatively transported w-level scalars. However, because this simple formulation does not guarantee boundedness of advected scalars, the user may choose whether or not to include the correction. The correction is included for the results shown in Sect. 5.

In the back-substitution stage of the time step, cheap updates are made to all transported terms [.]_D and [.]_T to account for the velocity increments u^′, applied over the backward part of the time step aΔt. The transported terms in the momentum equation are updated using an advective-form first-order upwind scheme. The transported term in the density equation is updated using a flux-form scheme with mass flux increments $\mathit{F}{}^{\prime }=\widehat{\mathit{\rho }}\mathit{u}{}^{\prime }$, where the cell edge values $\widehat{\mathit{\rho }}$ are given by a second-order centred scheme. The transported terms in any p-level tracer equations are updated using a flux-form scheme with tracer flux increments $\mathit{F}{}^{\prime }{}^{\mathit{\chi }}=\mathit{F}{}^{\prime }\widehat{\mathit{\chi }}$, where the cell edge values $\widehat{\mathit{\chi }}$ are given by a first-order upwind scheme. Finally, the transport terms for any w-level scalars are updated using a flux-form scheme with scalar flux increments $\mathit{F}{}^{\prime }{}^{\mathit{\chi }}={\stackrel{\mathrm{‾}}{\mathit{F}{}^{\prime }}}^d\widehat{\mathit{\chi }}$, with $\widehat{\mathit{\chi }}$ again given by a first-order upwind scheme. The notation ${\stackrel{\mathrm{‾}}{(.)}}^d$ here indicates that horizontal mass flux increments have been mapped to w levels using the ${\stackrel{\mathrm{‾}}{(.)}}^r$ operator and vertical mass flux increments have been mapped to p levels using the ${\stackrel{\mathrm{‾}}{(.)}}^p$ operator for consistency with the dual mass budget (Thuburn et al., 2022). This way of constructing the scalar flux increments ensures that the scalar transport remains consistent with density transport and conservative and bounded if the full semi-Lagrangian transport was conservative and bounded.

2.7 Boundary and interior forcing terms

Each simulated case requires the implementation of appropriate boundary and interior forcing terms. For clarity of the code, these boundary and forcing terms are implemented on a case-by-case basis. For the surface fluxes, the cases that are already implemented include examples of the setups most commonly used in LES: heat and moisture fluxes as specified functions of time, heat and moisture fluxes determined from a specified surface temperature and relative humidity using a simple bulk model, momentum fluxes determined by a simple bulk model, and momentum fluxes determined by Monin–Obukhov similarity theory.² A case using Monin–Obukhov theory to determine heat and moisture fluxes has not yet been implemented.

This section highlights some aspects of how the boundary and interior forcing terms are implemented for consistency with the rest of the model formulation.

2.7.1 Consistent surface fluxes of mass and water

Under a careful and consistent treatment of the moist thermodynamics, a surface flux of water implies a surface flux of mass. Most atmospheric models neglect that flux of mass. Because of the Charney–Phillips vertical staggering and the fact that the predicted density is the total density rather than the dry density, care is needed in computing the density and water increments to ensure that the changes in both total mass and total water within the model agree with the time-integrated surface fluxes.

First consider the forward part of the time step, i.e. the […]ⁿ terms in Eqs. (21) and (25). Let ℱ^ρ be the surface flux of mass and ${\stackrel{\mathrm{‾}}{{\mathit{\rho }}^n}}_{\mathrm{1}/\mathrm{2}}^r{\mathcal{F}}^q$ be the surface flux of water (both in kg m⁻² s⁻¹). Recall that the w-level density at level $k+\mathrm{1}/\mathrm{2}$ must equal the conservatively remapped p-level density ${\stackrel{\mathrm{‾}}{\mathit{\rho }}}_{k+\mathrm{1}/\mathrm{2}}^r$. Thus, any change in total density at level 1 affects the density of water at level $\mathrm{3}/\mathrm{2}$, ${\stackrel{\mathrm{‾}}{{\mathit{\rho }}^n}}_{\mathrm{3}/\mathrm{2}}^r{q}_{\mathrm{3}/\mathrm{2}}$, as well as the density of water at level $\mathrm{1}/\mathrm{2}$, ${\stackrel{\mathrm{‾}}{{\mathit{\rho }}^n}}_{\mathrm{1}/\mathrm{2}}^r{q}_{\mathrm{1}/\mathrm{2}}$. Nevertheless, imposing the constraints that

• only the density at level 1 may change, and the total change in column mass per unit area must equal ℱ^ρbΔt, and

• only the specific humidity at level $\mathrm{1}/\mathrm{2}$ may change, and the total change in column water mass per unit area must equal ${\stackrel{\mathrm{‾}}{{\mathit{\rho }}^n}}_{\mathrm{1}/\mathrm{2}}^r{\mathcal{F}}^{q}b\mathrm{\Delta }t$,

leads to unique solutions for the increment in density at level 1,

$$\begin{array}{}\text{(28)}& \mathit{\delta }{\mathit{\rho }}_{\mathrm{1}}={\displaystyle \frac{{\mathcal{F}}^{\mathit{\rho }}b\mathrm{\Delta }t}{\mathrm{\Delta }{z}_{\mathrm{1}}}},\end{array}$$

and the increment in specific humidity at level $\mathrm{1}/\mathrm{2}$,

$$\begin{array}{}\text{(29)}& \mathit{\delta }{q}_{\mathrm{1}/\mathrm{2}}={\displaystyle \frac{\left({\mathit{\rho }}_{\mathrm{1}}^n{\mathcal{F}}^q-{\stackrel{\mathrm{‾}}{q^n}}_{\mathrm{1}}^s{\mathcal{F}}^{\mathit{\rho }}\right)b\mathrm{\Delta }t}{{\mathit{\rho }}_{\mathrm{1}}^{+}\mathrm{\Delta }{z}_{\mathrm{1}/\mathrm{2}}}},\end{array}$$

where ${\mathit{\rho }}_{\mathrm{1}}^{+}={\mathit{\rho }}_{\mathrm{1}}^n+\mathit{\delta }{\mathit{\rho }}_{\mathrm{1}}$, and we have used the fact that ${\stackrel{\mathrm{‾}}{\mathit{\rho }}}_{\mathrm{1}/\mathrm{2}}^r\equiv {\mathit{\rho }}_{\mathrm{1}}$.

For the backward part of the time step, i.e. the […]ⁿ⁺¹ terms in Eqs. (21) and (25), an analogous argument leads to

$$\begin{array}{}\text{(30)}& \mathit{\delta }{\mathit{\rho }}_{\mathrm{1}}={\displaystyle \frac{{\mathcal{F}}^{\mathit{\rho }}a\mathrm{\Delta }t}{\mathrm{\Delta }{z}_{\mathrm{1}}}}\end{array}$$

and

$$\begin{array}{}\text{(31)}& \mathit{\delta }{q}_{\mathrm{1}/\mathrm{2}}={\displaystyle \frac{\left({\mathit{\rho }}_{\mathrm{1}}^{n+\mathrm{1}}{\mathcal{F}}^q-{\stackrel{\mathrm{‾}}{q^{n+\mathrm{1}}}}_{\mathrm{1}}^s{\mathcal{F}}^{\mathit{\rho }}\right)a\mathrm{\Delta }t}{{\mathit{\rho }}_{\mathrm{1}}^{-}\mathrm{\Delta }{z}_{\mathrm{1}/\mathrm{2}}}},\end{array}$$

where ${\mathit{\rho }}_{\mathrm{1}}^{-}={\mathit{\rho }}_{\mathrm{1}}^{n+\mathrm{1}}-\mathit{\delta }{\mathit{\rho }}_{\mathrm{1}}$. Terms at time step n+1 are approximated by the latest available estimate.

The treatment of surface entropy fluxes is analogous, with q replaced by η and ℱ^q replaced by ℱ^η.

2.7.2 Surface drag

The various test cases implemented specify how the vertical flux of horizontal momentum is to be parameterized at the surface. Since there are no parameterized subgrid fluxes in the interior of the domain, the convergence of the parameterized momentum flux occurs entirely in the lowest model layer; i.e. it is applied to the u and v components at model level 1.

Some of the test cases implemented specify that the surface momentum flux is to be computed using Monin–Obukhov similarity theory. However, there are two difficulties. First, Monin–Obukhov similarity theory expresses the mean flow speed U(z) at some height z as a nonlinear function of the friction velocity U_*. What we require in a numerical model is to express the friction velocity, hence the surface momentum flux, in terms of the known flow speed at the lowest model level; however, in the usual approach, this requires the iterative solution of a nonlinear problem in each grid column. Second, and more seriously, in the stable case (negative surface buoyancy flux) in light winds Monin–Obukhov similarity theory is valid only over a limited height range that might be less than the height of the lowest model level z₁. In this case, with commonly used stability functions, the theory can produce either no solution or multiple solutions for U_*, given U(z₁). To avoid these difficulties, a slight reformulation of Monin–Obukhov theory is used to compute surface momentum fluxes (Appendix A4).

2.7.3 Consistent interior forcing terms

For test cases that require a source of moisture in the interior of the model domain, consistency between the p-level mass budget and the w-level dual mass budget requires that

$$\begin{array}{}\text{(32)}& {\left(\mathit{\rho }{S}_q\right)}_{k+\mathrm{1}/\mathrm{2}}={\left({\stackrel{\mathrm{‾}}{S_{\mathit{\rho }}}}^r\right)}_{k+\mathrm{1}/\mathrm{2}}.\end{array}$$

The simplest way to satisfy this condition is to specify the required moisture/mass source first on p levels and then remap conservatively to w levels. For simplicity, interior entropy sources are specified in the same way.

2.8 Linearized dynamics

To take a model time step, a quasi-Newton method is used to solve Eqs. (21)–(26) along with the p-level and w-level thermodynamic equations. After some number of quasi-Newton iterations, using the latest available estimates to evaluate step n+1 terms, the residuals R_ρ, R_v, R_w, R_ρη, R_ρq, R_PE, and $R_{\stackrel{\mathrm{‾}}{p}}$ will generally be non-zero. We seek increments or updates to the step n+1 model variables that will reduce those residuals. This is done through an approximate linearization of the governing equations (this section and Sect. 2.9), leading to a large linear system for the increments. Despite the apparent complexity of this system, systematic elimination of unknowns, partly manually and partly numerically, leads to a nearly standard Helmholtz problem for a single unknown per model grid point: the pressure increment p^′ (Sect. 2.10). The Helmholtz problem can be solved efficiently using well-established methods; here a multigrid method is used. Once p^′ is found, the other increments can be found through back substitution (Sect. 2.11).

The approximate linearization of Eqs. (21)–(26) is

$$\begin{array}{}\text{(33)}& {\displaystyle }{\mathit{\rho }}^{\prime }+a\mathrm{\Delta }t\mathrm{\nabla }\cdot {\mathit{F}}^{\prime }=R_{\mathit{\rho }},\text{(34)}& {\displaystyle }{\mathit{v}}^{\prime }+a\mathrm{\Delta }t\left({\mathit{\alpha }}_{*}{\mathrm{\nabla }}_{\mathbf{H}}{p}^{\prime }-{\mathit{S}}_{\mathit{v}}^{\prime }\right)={\mathit{R}}_{\mathit{v}},\text{(35)}& {\displaystyle }{w}^{\prime }+a\mathrm{\Delta }t\left({\mathit{\alpha }}_{*}{\displaystyle \frac{\partial p^{\prime }}{\partial z}}+{\mathit{\alpha }}^{\prime }{\displaystyle \frac{\partial p_{*}}{\partial z}}\right)=R_w,\end{array}$$

$$\begin{array}{}\text{(36)}& {\displaystyle }\begin{array}{rl}{\stackrel{\mathrm{‾}}{{\mathit{\rho }}^{(\mathrm{\ell }+\mathrm{1})}}}^r{\mathit{\eta }}^{\prime }+{\stackrel{\mathrm{‾}}{{\mathit{\rho }}^{\prime }}}^r{\left(\mathit{\eta }-a\mathrm{\Delta }t{S}_{\mathit{\eta }}\right)}^{\left(\mathrm{\ell }\right)}& -{\left[{\stackrel{\mathrm{‾}}{\mathit{\rho }}}^r\left(\mathit{\eta }+b\mathrm{\Delta }t{S}_{\mathit{\eta }}\right)\right]}_{\mathrm{T}}^{\prime }\\ & =R_{\mathit{\rho }\mathit{\eta }},\end{array}\text{(37)}& {\displaystyle }\begin{array}{rl}{\stackrel{\mathrm{‾}}{{\mathit{\rho }}^{(\mathrm{\ell }+\mathrm{1})}}}^r{q}^{\prime }+{\stackrel{\mathrm{‾}}{{\mathit{\rho }}^{\prime }}}^r{\left(q-a\mathrm{\Delta }t{S}_q\right)}^{\left(\mathrm{\ell }\right)}& -{\left[{\stackrel{\mathrm{‾}}{\mathit{\rho }}}^r\left(q+b\mathrm{\Delta }t{S}_q\right)\right]}_{\mathrm{T}}^{\prime }\\ & =R_{\mathit{\rho }q},\end{array}\text{(38)}& {\displaystyle }{\mathit{X}}^{\prime }-\mathrm{\Delta }t{\mathit{J}}^{\prime }=\mathrm{0}.\end{array}$$

An asterisk on any variable indicates a reference value for the linearization. Here the most recent estimate for the step n+1 value is used as the reference value, avoiding the need to store additional three-dimensional fields. Certain terms in Eqs. (36) and (37) are explicitly evaluated at iteration number ℓ or ℓ+1. This specific way of writing the linearization ensures that, when Eqs. (36) and (37) are used in the back substitution, η and q are incremented exactly according to the consistent and conservative transport updates. Since the source terms S_ρ are generally non-stiff, their linearizations are mostly omitted. The exception is a linearization of the surface drag at model level 1 ${\mathit{S}}_{\mathit{v}}^{\prime }=-{\mathit{v}}^{\prime }/{\mathit{\tau }}_{\mathrm{drag}}$, where τ_drag is a surface drag timescale computed alongside the surface momentum flux. This term is omitted from the Helmholtz problem (in principle it could be included) but is included in the back substitution.

In order to arrive at a standard Helmholtz problem for the pressure increments we must make a w^′N² term appear in the linearization, where N² is an appropriately defined buoyancy frequency squared. This requires us to work with a linearization of advective-form transport equations for η and q. Consider the η equation (Eq. 36). In the back substitution the transport increment is computed as

$$\begin{array}{}\text{(39)}& {\left[{\stackrel{\mathrm{‾}}{\mathit{\rho }}}^r\left(\mathit{\eta }+b\mathrm{\Delta }t{S}_{\mathit{\eta }}\right)\right]}_{\mathrm{T}}^{\prime }=-a\mathrm{\Delta }t\mathrm{\nabla }\cdot \left({\stackrel{\mathrm{‾}}{{\mathit{F}}^{\prime }}}^d\widehat{\mathit{\eta }}\right)\end{array}$$

for some (upwind) cell edge values $\widehat{\mathit{\eta }}$. Thus, Eq. (36) becomes

$$\begin{array}{}\text{(40)}& {\stackrel{\mathrm{‾}}{{\mathit{\rho }}^{(\mathrm{\ell }+\mathrm{1})}}}^r{\mathit{\eta }}^{\prime }+{\stackrel{\mathrm{‾}}{{\mathit{\rho }}^{\prime }}}^r{\left(\mathit{\eta }-a\mathrm{\Delta }t{S}_{\mathit{\eta }}\right)}^{\left(\mathrm{\ell }\right)}+a\mathrm{\Delta }t\mathrm{\nabla }\cdot \left({\stackrel{\mathrm{‾}}{{\mathit{F}}^{\prime }}}^d\widehat{\mathit{\eta }}\right)=R_{\mathit{\rho }\mathit{\eta }}.\end{array}$$

Subtracting $\widehat{\mathit{\eta }}$ times a vertical average of Eq. (33) and approximating ${\stackrel{\mathrm{‾}}{{\mathit{\rho }}^{(\mathrm{\ell }+\mathrm{1})}}}^r$ by a generic reference density ρ_* gives

$$\begin{array}{}\text{(41)}& {\mathit{\eta }}^{\prime }+a\mathrm{\Delta }t{w}^{\prime }{\displaystyle \frac{\partial \mathit{\eta }}{\partial z}}={\displaystyle \frac{R_{\mathit{\rho }\mathit{\eta }}-\widehat{\mathit{\eta }}\stackrel{\mathrm{‾}}{R_{\mathit{\rho }}}}{{\mathit{\rho }}_{*}}}\equiv R_{\mathit{\eta }}.\end{array}$$

Proceeding in a similar way for the q equation (Eq. 37) gives

$$\begin{array}{}\text{(42)}& q^{\prime }+a\mathrm{\Delta }t{w}^{\prime }{\displaystyle \frac{\partial q}{\partial z}}={\displaystyle \frac{R_{\mathit{\rho }q}-\widehat{q}\stackrel{\mathrm{‾}}{R_{\mathit{\rho }}}}{{\mathit{\rho }}_{*}}}\equiv R_q.\end{array}$$

2.9 Linearized thermodynamics

When the w-level thermodynamic state is updated, certain fields are diagnosed, ensuring that Eqs. (6)–(11) are satisfied exactly: q^d is set equal to 1−q, Eq. (9) gives q^v, Eq. (10) gives α^v, and Eq. (6) gives α^d. Also, the updating of terms in Eqs. (7) and (8) during back substitution is done in such a way that those equations remain exactly satisfied. Thus, no residual term appears in the linearized versions of these equations.

The linearized forms of Eqs. (6)–(11) are then as follows:

$$\begin{array}{}\text{(43)}& {\displaystyle }{q}^{\mathrm{v}}{\mathit{\alpha }}^{\mathrm{v}}{}^{\prime }+{\mathit{\alpha }}^{\mathrm{v}}{q}^{\mathrm{v}}{}^{\prime }-q^{\mathrm{d}}{\mathit{\alpha }}^{\mathrm{d}}{}^{\prime }+{\mathit{\alpha }}^{\mathrm{d}}{q}^{\prime }=\mathrm{0},\text{(44)}& {\displaystyle }\left(\mathrm{1}-{\mathit{\delta }}^{\mathrm{l}}\right){\mathit{\lambda }}^{\mathrm{l}}{}^{\prime }-{\mathit{\delta }}^{\mathrm{l}}{q}^{\mathrm{l}}{}^{\prime }=\mathrm{0},\text{(45)}& {\displaystyle }\left(\mathrm{1}-{\mathit{\delta }}^{\mathrm{f}}\right){\mathit{\lambda }}^{\mathrm{f}}{}^{\prime }-{\mathit{\delta }}^{\mathrm{f}}{q}^{\mathrm{f}}{}^{\prime }=\mathrm{0},\text{(46)}& {\displaystyle }{q}^{\mathrm{v}}{}^{\prime }+q^{\mathrm{l}}{}^{\prime }+q^{\mathrm{f}}{}^{\prime }-q^{\prime }=\mathrm{0},\text{(47)}& {\displaystyle }{q}^{\mathrm{v}}{\mathit{\alpha }}^{\mathrm{v}}{}^{\prime }+{\mathit{\alpha }}^{\mathrm{v}}{q}^{\mathrm{v}}{}^{\prime }+{\mathit{\alpha }}^{\mathrm{l}}{q}^{\mathrm{l}}{}^{\prime }+{\mathit{\alpha }}^{\mathrm{f}}{q}^{\mathrm{f}}{}^{\prime }-{\mathit{\alpha }}^{\prime }=\mathrm{0},\text{(48)}& {\displaystyle }\begin{array}{rl}{q}^{\mathrm{d}}{\mathit{\eta }}^{\mathrm{d}}{}^{\prime }& -{\mathit{\eta }}^{\mathrm{d}}{q}^{\prime }+q^{\mathrm{v}}{\mathit{\eta }}^{\mathrm{v}}{}^{\prime }+{\mathit{\eta }}^{\mathrm{v}}{q}^{\mathrm{v}}{}^{\prime }+q^{\mathrm{l}}{\mathit{\eta }}^{\mathrm{l}}{}^{\prime }+{\mathit{\eta }}^{\mathrm{l}}{q}^{\mathrm{l}}{}^{\prime }+q^{\mathrm{f}}{\mathit{\eta }}^{\mathrm{f}}{}^{\prime }\\ & +{\mathit{\eta }}^{\mathrm{f}}{q}^{\mathrm{f}}{}^{\prime }-{\mathit{\eta }}^{\prime }=\mathrm{0}.\end{array}\end{array}$$

Note that $q^{\mathrm{d}}{}^{\prime }=-q^{\prime }$ has been used to eliminate increments of the dry mass fraction. The linearization of Eq. (17) is

$$\begin{array}{}\text{(49)}& {\stackrel{\mathrm{‾}}{p^{\prime }}}^w-p^{\left(w\right)}{}^{\prime }=R_{\stackrel{\mathrm{‾}}{p}}.\end{array}$$

With the aid of Eq. (38), the linearized phenomenological equations become

$$\begin{array}{}\text{(50)}& {\displaystyle \frac{\mathrm{1}}{\mathrm{\Delta }t}}\mathbf{R}{\mathit{X}}^{\prime }-{\mathit{P}}^{\prime }+{\mathit{C}}^{\prime }={\mathit{R}}_{\mathrm{PE}};\end{array}$$

for brevity the details of P^′ and C^′ are suppressed.

Two modifications are then made to the linear system. First, a change of variable is made, introducing ${\stackrel{\mathrm{̃}}{\mathit{\lambda }}}^{\mathrm{l}}{}^{\prime }={\mathit{\lambda }}^{\mathrm{l}}{}^{\prime }+q^{\mathrm{l}}{}^{\prime }$ and ${\stackrel{\mathrm{̃}}{\mathit{\lambda }}}^{\mathrm{f}}{}^{\prime }={\mathit{\lambda }}^{\mathrm{f}}{}^{\prime }+q^{\mathrm{f}}{}^{\prime }$. Second, η_ref times Eq. (44) is added to the q^l′ phenomenological equation and η_ref times Eq. (45) is added to the q^f′ phenomenological equation. These two modifications guarantee that the coefficient of q^l′ in Eq. (44) and the coefficient of q^f′ in Eq. (45) are nonzero and that the coefficients of ${\stackrel{\mathrm{̃}}{\mathit{\lambda }}}^{\mathrm{l}}{}^{\prime }$ and ${\stackrel{\mathrm{̃}}{\mathit{\lambda }}}^{\mathrm{f}}{}^{\prime }$ in the relevant phenomenological equations are nonzero. In this way, the sparsity pattern of the system matrix is known irrespective of the state of the switches δ^l and δ^f and is then effectively the same at all grid points.

The resulting system of linear equations may be compactly written

$$\begin{array}{}\text{(51)}& \stackrel{\mathrm{̃}}{\mathbf{M}}{\mathit{Z}}^{\prime }={\stackrel{\mathrm{̃}}{\mathit{R}}}_{\mathbf{M}},\end{array}$$

where $\stackrel{\mathrm{̃}}{\mathbf{M}}$ is an 11×14 matrix, and $({\mathit{Z}}^{\prime }{)}^{\mathrm{T}}=(({\mathit{Y}}^{\prime }{)}^{\mathrm{T}},q^{\prime },{\mathit{\alpha }}^{\prime },{\mathit{\eta }}^{\prime })$, where $({\mathit{Y}}^{\prime }{)}^{\mathrm{T}}=({\mathit{\alpha }}^{\mathrm{v}}{}^{\prime },q^{\mathrm{l}}{}^{\prime },q^{\mathrm{f}}{}^{\prime },q^{\mathrm{v}}{}^{\prime },{\mathit{\alpha }}^{\mathrm{d}}{}^{\prime },{\mathit{\eta }}^{\mathrm{d}}{}^{\prime },{\mathit{\eta }}^{\mathrm{v}}{}^{\prime },{\mathit{\eta }}^{\mathrm{l}}{}^{\prime },{\mathit{\eta }}^{\mathrm{f}}{}^{\prime },{\stackrel{\mathrm{̃}}{\mathit{\lambda }}}^{\mathrm{l}}{}^{\prime },{\stackrel{\mathrm{̃}}{\mathit{\lambda }}}^{\mathrm{f}}{}^{\prime })$. The different ordering of the rows and columns of $\stackrel{\mathrm{̃}}{\mathbf{M}}$ compared to Bowen and Thuburn (2022a) allows a better exploitation of the sparsity in the Gaussian elimination step discussed in Sect. 2.10.

The linearized versions of the p-level thermodynamic equations (Eqs. 18–20) are

$$\begin{array}{}\text{(52)}& {\displaystyle }& {\displaystyle }{\stackrel{\mathrm{‾}}{q^{\mathrm{v}}{}^{\prime }}}^p{\mathit{\alpha }}^{\mathrm{v}\left(p\right)}+{\stackrel{\mathrm{‾}}{q^{\mathrm{v}}}}^p{\mathit{\alpha }}^{\mathrm{v}\left(p\right)}{}^{\prime }+{\stackrel{\mathrm{‾}}{q^{\mathrm{l}}{}^{\prime }}}^p{\mathit{\alpha }}^{\mathrm{l}}+{\stackrel{\mathrm{‾}}{q^{\mathrm{f}}{}^{\prime }}}^p{\mathit{\alpha }}^{\mathrm{f}}+{\displaystyle \frac{{\mathit{\rho }}^{\prime }}{{\mathit{\rho }}_{*}^{\mathrm{2}}}}=\mathrm{0},\text{(53)}& {\displaystyle }& {\displaystyle }{\stackrel{\mathrm{‾}}{q^{\mathrm{v}}{}^{\prime }}}^p{\mathit{\alpha }}^{\mathrm{v}\left(p\right)}+{\stackrel{\mathrm{‾}}{q^{\mathrm{v}}}}^p{\mathit{\alpha }}^{\mathrm{v}\left(p\right)}{}^{\prime }+{\stackrel{\mathrm{‾}}{q{}^{\prime }}}^p{\mathit{\alpha }}^{\mathrm{d}\left(p\right)}-{\stackrel{\mathrm{‾}}{q^{\mathrm{d}}}}^p{\mathit{\alpha }}^{\mathrm{d}\left(p\right)}{}^{\prime }=\mathrm{0},\text{(54)}& {\displaystyle }& {\displaystyle }{p}^{\prime }+e_{\mathit{\alpha }\mathit{\alpha }}^{\mathrm{d}}{\mathit{\alpha }}^{\mathrm{d}\left(p\right){}^{\prime }}+e_{\mathit{\alpha }\mathit{\eta }}^{\mathrm{d}}{\stackrel{\mathrm{‾}}{{\mathit{\eta }}^{\mathrm{d}}{}^{\prime }}}^p+e_{\mathit{\alpha }\mathit{\alpha }}^{\mathrm{v}}{\mathit{\alpha }}^{\mathrm{v}\left(p\right){}^{\prime }}+e_{\mathit{\alpha }\mathit{\eta }}^{\mathrm{v}}{\stackrel{\mathrm{‾}}{{\mathit{\eta }}^{\mathrm{v}}{}^{\prime }}}^p=\mathrm{0}.\end{array}$$

2.10 Derivation of the Helmholtz problem

In order to derive a Helmholtz problem, it will be useful to express all other thermodynamic increments at w levels (the components of Y^′) in terms of q^′, α^′, and η^′. This is done by carrying out a numerical Gaussian elimination on $\stackrel{\mathrm{̃}}{\mathbf{M}}$ to leave

$$\begin{array}{}\text{(55)}& \mathbf{M}{\mathit{Z}}^{\prime }={\mathit{R}}_{\mathbf{M}},\end{array}$$

where M is of the form

$$\begin{array}{}\text{(56)}& \mathbf{M}=\left(\mathbf{I}{\mathit{C}}_{\mathrm{1}}{\mathit{C}}_{\mathrm{2}}{\mathit{C}}_{\mathrm{3}}\right).\end{array}$$

I is the 11×11 identity matrix, and C₁, C₂, and C₃ are columns of (generally) nonzero entries. For efficiency, the Gaussian elimination exploits the known sparsity pattern of the matrix $\stackrel{\mathrm{̃}}{\mathbf{M}}$, and, since the sparsity pattern is the same for all grid points, the elimination can be done without conditional code and can be implemented with the innermost loop over model levels.

The entries M_ij of the eliminated matrix M and the entries R_Mi of the eliminated right-hand side R_M are used in building the coefficients and right-hand side of the Helmholtz problem. A subset of them (rows 2, 3, and 7–11 of C₁, C₂, C₃, and R_M, corresponding to the equations for (${\mathit{\lambda }}^{\mathrm{l}}{}^{\prime },{\mathit{\lambda }}^{\mathrm{f}}{}^{\prime },q^{\mathrm{l}}{}^{\prime },q^{\mathrm{f}}{}^{\prime },{\mathit{\eta }}^{\mathrm{v}}{}^{\prime },{\mathit{\eta }}^{\mathrm{l}}{}^{\prime },{\mathit{\eta }}^{\mathrm{f}}{}^{\prime }$)) are saved for use in back substitution.

Next, we need an equation relating p^′, α^′, and w^′ at w levels and an analogous equation relating p^′, ρ^′, and w^′ at p levels. With the aid of Eq. (55), the pressure perturbation at w levels is given by

$$\begin{array}{}\text{(57)}& \begin{array}{rl}{p}^{\left(w\right)}{}^{\prime }& =-e_{\mathit{\alpha }\mathit{\alpha }}^{\mathrm{d}}{\mathit{\alpha }}^{\mathrm{d}}{}^{\prime }-e_{\mathit{\alpha }\mathit{\eta }}^{\mathrm{d}}{\mathit{\eta }}^{\mathrm{d}}{}^{\prime }-e_{\mathit{\alpha }\mathit{\alpha }}^{\mathrm{v}}{\mathit{\alpha }}^{\mathrm{v}}{}^{\prime }-e_{\mathit{\alpha }\mathit{\eta }}^{\mathrm{v}}{\mathit{\eta }}^{\mathrm{v}}{}^{\prime }\\ & ={\displaystyle \frac{\partial p}{\partial q}}{q}^{\prime }+{\displaystyle \frac{\partial p}{\partial \mathit{\alpha }}}{\mathit{\alpha }}^{\prime }+{\displaystyle \frac{\partial p}{\partial \mathit{\eta }}}{\mathit{\eta }}^{\prime }+R_{\mathrm{TDP}},\end{array}\end{array}$$

where

$$\begin{array}{}\text{(58)}& {\displaystyle }& {\displaystyle \frac{\partial p}{\partial q}}=e_{\mathit{\alpha }\mathit{\alpha }}^{\mathrm{d}}{M}_{\mathrm{5}\mathrm{12}}+e_{\mathit{\alpha }\mathit{\eta }}^{\mathrm{d}}{M}_{\mathrm{6}\mathrm{12}}+e_{\mathit{\alpha }\mathit{\alpha }}^{\mathrm{v}}{M}_{\mathrm{1}\mathrm{12}}+e_{\mathit{\alpha }\mathit{\eta }}^{\mathrm{v}}{M}_{\mathrm{7}\mathrm{12}},\text{(59)}& {\displaystyle }& {\displaystyle \frac{\partial p}{\partial \mathit{\alpha }}}=e_{\mathit{\alpha }\mathit{\alpha }}^{\mathrm{d}}{M}_{\mathrm{5}\mathrm{13}}+e_{\mathit{\alpha }\mathit{\eta }}^{\mathrm{d}}{M}_{\mathrm{6}\mathrm{13}}+e_{\mathit{\alpha }\mathit{\alpha }}^{\mathrm{v}}{M}_{\mathrm{1}\mathrm{13}}+e_{\mathit{\alpha }\mathit{\eta }}^{\mathrm{v}}{M}_{\mathrm{7}\mathrm{13}},\text{(60)}& {\displaystyle }& {\displaystyle \frac{\partial p}{\partial \mathit{\eta }}}=e_{\mathit{\alpha }\mathit{\alpha }}^{\mathrm{d}}{M}_{\mathrm{5}\mathrm{14}}+e_{\mathit{\alpha }\mathit{\eta }}^{\mathrm{d}}{M}_{\mathrm{6}\mathrm{14}}+e_{\mathit{\alpha }\mathit{\alpha }}^{\mathrm{v}}{M}_{\mathrm{1}\mathrm{14}}+e_{\mathit{\alpha }\mathit{\eta }}^{\mathrm{v}}{M}_{\mathrm{7}\mathrm{14}},\text{(61)}& {\displaystyle }& {\displaystyle }{R}_{\mathrm{TDP}}=e_{\mathit{\alpha }\mathit{\alpha }}^{\mathrm{d}}{R}_{\mathbf{M}\mathrm{5}}+e_{\mathit{\alpha }\mathit{\eta }}^{\mathrm{d}}{R}_{\mathbf{M}\mathrm{6}}+e_{\mathit{\alpha }\mathit{\alpha }}^{\mathrm{v}}{R}_{\mathbf{M}\mathrm{1}}+e_{\mathit{\alpha }\mathit{\eta }}^{\mathrm{v}}{R}_{\mathbf{M}\mathrm{7}}.\end{array}$$

Defining the sound speed c by

$$\begin{array}{}\text{(62)}& {\displaystyle \frac{\partial p}{\partial \mathit{\alpha }}}=-{\mathit{\rho }}_{*}^{\mathrm{2}}{c}^{\mathrm{2}}\end{array}$$

and using Eqs. (36) and (37) to eliminate η^′ and q^′, Eq. (49) becomes

$$\begin{array}{}\text{(63)}& {\displaystyle \frac{{\stackrel{\mathrm{‾}}{p^{\prime }}}^w}{c^{\mathrm{2}}}}+{\mathit{\rho }}_{*}^{\mathrm{2}}{\mathit{\alpha }}^{\prime }+a\mathrm{\Delta }t{\mathit{\rho }}_{*}{w}^{\prime }{\displaystyle \frac{N^{\mathrm{2}}}{\mathrm{g}}}=R_{\mathrm{buoy}},\end{array}$$

where

$$\begin{array}{}\text{(64)}& {\mathit{\rho }}_{*}{\displaystyle \frac{N^{\mathrm{2}}}{g}}={\displaystyle \frac{\mathrm{1}}{c^{\mathrm{2}}}}\left\{{\displaystyle \frac{\partial p}{\partial q}}{\displaystyle \frac{\partial q}{\partial z}}+{\displaystyle \frac{\partial p}{\partial \mathit{\eta }}}{\displaystyle \frac{\partial \mathit{\eta }}{\partial z}}\right\}\end{array}$$

and

$$\begin{array}{}\text{(65)}& R_{\mathrm{buoy}}={\displaystyle \frac{\mathrm{1}}{c^{\mathrm{2}}}}\left\{R_{\stackrel{\mathrm{‾}}{p}}+R_{\mathrm{TDP}}+{\displaystyle \frac{\partial p}{\partial q}}{R}_q+{\displaystyle \frac{\partial p}{\partial \mathit{\eta }}}{R}_{\mathit{\eta }}\right\}.\end{array}$$

Deriving the analogous equation on p levels is a little more subtle because p^′ depends both on ρ^′ and, via q^i′ and η^i′, on the w-level α^′ averaged to p levels. Setting q^i′ and η^i′ to zero in Eqs. (52)–(54) shows that

$$\begin{array}{}\text{(66)}& {\left.{\displaystyle \frac{\partial p}{\partial \mathit{\rho }}}\right|}_{{\stackrel{\mathrm{‾}}{\mathit{\alpha }}}^p}={\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{*}^{\mathrm{2}}}}\left({\displaystyle \frac{e_{\mathit{\alpha }\mathit{\alpha }}^{\mathrm{d}}}{{\stackrel{\mathrm{‾}}{q^{\mathrm{d}}}}^p}}+{\displaystyle \frac{e_{\mathit{\alpha }\mathit{\alpha }}^{\mathrm{v}}}{{\stackrel{\mathrm{‾}}{q^{\mathrm{v}}}}^p}}\right)\equiv {\widehat{c}}^{\mathrm{2}},\end{array}$$

where $\widehat{c}$ is a reduced sound speed (typically very close to and slightly larger than c). Hence

$$\begin{array}{}\text{(67)}& {\left.{\displaystyle \frac{\partial p}{\partial {\stackrel{\mathrm{‾}}{\mathit{\alpha }}}^p}}\right|}_{\mathit{\rho }}=-{\mathit{\rho }}_{*}^{\mathrm{2}}\left(c^{\mathrm{2}}-{\widehat{c}}^{\mathrm{2}}\right),\end{array}$$

while $\partial p/\partial q$ and $\partial p/\partial \mathit{\eta }$ are given sufficiently accurately by their w-level values averaged to p levels. Proceeding in this way, Eq. (54) becomes

$$\begin{array}{}\text{(68)}& p{}^{\prime }-{\widehat{c}}^{\mathrm{2}}\mathit{\rho }{}^{\prime }+{\mathit{\rho }}_{*}^{\mathrm{2}}\left(c^{\mathrm{2}}-{\widehat{c}}^{\mathrm{2}}\right){\stackrel{\mathrm{‾}}{\mathit{\alpha }{}^{\prime }}}^p-{\displaystyle \frac{\partial p}{\partial q}}{\stackrel{\mathrm{‾}}{q{}^{\prime }}}^p-{\displaystyle \frac{\partial p}{\partial \mathit{\eta }}}{\stackrel{\mathrm{‾}}{\mathit{\eta }{}^{\prime }}}^p={\stackrel{\mathrm{‾}}{R_{\mathrm{TDP}}}}^p.\end{array}$$

Using Eq. (63) to eliminate α^′ and Eqs. (36) and (37) to eliminate η^′ and q^′ leaves

$$\begin{array}{}\text{(69)}& \begin{array}{rl}{\displaystyle \frac{p{}^{\prime }}{{\widehat{c}}^{\mathrm{2}}}}& -\mathit{\rho }{}^{\prime }+\left({\displaystyle \frac{\mathrm{1}}{c^{\mathrm{2}}}}-{\displaystyle \frac{\mathrm{1}}{{\widehat{c}}^{\mathrm{2}}}}\right){\stackrel{\mathrm{‾}}{{\stackrel{\mathrm{‾}}{p{}^{\prime }}}^w}}^p+a\mathrm{\Delta }t{\stackrel{\mathrm{‾}}{{\mathit{\rho }}_{*}{w}^{\prime }{\displaystyle \frac{N^{\mathrm{2}}}{g}}}}^p\\ & ={\stackrel{\mathrm{‾}}{R_{\mathrm{buoy}}}}^p-{\displaystyle \frac{{\stackrel{\mathrm{‾}}{R_{\stackrel{\mathrm{‾}}{p}}}}^p}{{\widehat{c}}^{\mathrm{2}}}}.\end{array}\end{array}$$

Next, use Eq. (63) to eliminate α^′ from Eq. (35) and combine the terms involving p^′ into a single operator:

$$\begin{array}{}\text{(70)}& {\mathit{\rho }}_{*}{w}^{\prime }+{\mathcal{D}}_{\mathrm{1}}\left(p^{\prime }\right)={\displaystyle \frac{{\mathit{\rho }}_{*}{R}_w+a\mathrm{\Delta }tg{R}_{\mathrm{buoy}}}{\mathrm{1}+a^{\mathrm{2}}\mathrm{\Delta }{t}^{\mathrm{2}}{N}^{\mathrm{2}}}}\equiv {\mathit{\rho }}_{*}{R}_{wDp},\end{array}$$

where

$$\begin{array}{}\text{(71)}& {\mathcal{D}}_{\mathrm{1}}\left(p^{\prime }\right)\equiv {\displaystyle \frac{a\mathrm{\Delta }t}{\mathrm{1}+a^{\mathrm{2}}\mathrm{\Delta }{t}^{\mathrm{2}}{N}^{\mathrm{2}}}}\left({\displaystyle \frac{\partial p^{\prime }}{\partial z}}+{\displaystyle \frac{\mathrm{g}}{c^{\mathrm{2}}}}{\stackrel{\mathrm{‾}}{p^{\prime }}}^w\right).\end{array}$$

Also, use Eq. (69) to eliminate ρ^′ from Eq. (33) and combine the terms involving w^′ into a single operator:

$$\begin{array}{}\text{(72)}& \begin{array}{rl}{\displaystyle \frac{p{}^{\prime }}{{\widehat{c}}^{\mathrm{2}}}}& +\left({\displaystyle \frac{\mathrm{1}}{c^{\mathrm{2}}}}-{\displaystyle \frac{\mathrm{1}}{{\widehat{c}}^{\mathrm{2}}}}\right){\stackrel{\mathrm{‾}}{{\stackrel{\mathrm{‾}}{p{}^{\prime }}}^w}}^p+a\mathrm{\Delta }t{\mathrm{\nabla }}_{\mathbf{H}}\cdot \left({\mathit{\rho }}_{*}\mathit{v}{}^{\prime }\right)+{\mathcal{D}}_{\mathrm{2}}\left({\mathit{\rho }}_{*}w{}^{\prime }\right)\\ & ={\stackrel{\mathrm{‾}}{R_{\mathrm{buoy}}}}^p-{\displaystyle \frac{{\stackrel{\mathrm{‾}}{R_{\stackrel{\mathrm{‾}}{p}}}}^p}{{\widehat{c}}^{\mathrm{2}}}}+R_{\mathit{\rho }},\end{array}\end{array}$$

where

$$\begin{array}{}\text{(73)}& {\mathcal{D}}_{\mathrm{2}}\left({\mathit{\rho }}_{*}{w}^{\prime }\right)\equiv a\mathrm{\Delta }t\left({\displaystyle \frac{\partial }{\partial z}}{\mathit{\rho }}_{*}{w}^{\prime }+{\displaystyle \frac{N^{\mathrm{2}}}{g}}{\stackrel{\mathrm{‾}}{{\mathit{\rho }}_{*}{w}^{\prime }}}^p\right).\end{array}$$

Finally, eliminate u^′, v^′, and w^′ using Eqs. (34) and (70) to obtain the Helmholtz problem

$$\begin{array}{}\text{(74)}& \begin{array}{rl}{\displaystyle \frac{p{}^{\prime }}{{\widehat{c}}^{\mathrm{2}}}}& +\left({\displaystyle \frac{\mathrm{1}}{c^{\mathrm{2}}}}-{\displaystyle \frac{\mathrm{1}}{{\widehat{c}}^{\mathrm{2}}}}\right){\stackrel{\mathrm{‾}}{{\stackrel{\mathrm{‾}}{p{}^{\prime }}}^w}}^p-a^{\mathrm{2}}\mathrm{\Delta }{t}^{\mathrm{2}}{\mathrm{\nabla }}_{\mathbf{H}}\cdot {\mathrm{\nabla }}_{\mathbf{H}}p{}^{\prime }-{\mathcal{D}}_{\mathrm{2}}{\mathcal{D}}_{\mathrm{1}}\left(p{}^{\prime }\right)\\ & ={\stackrel{\mathrm{‾}}{R_{\mathrm{buoy}}}}^p-{\displaystyle \frac{{\stackrel{\mathrm{‾}}{R_{\stackrel{\mathrm{‾}}{p}}}}^p}{{\widehat{c}}^{\mathrm{2}}}}+R_{\mathit{\rho }}-a\mathrm{\Delta }t{\mathrm{\nabla }}_{\mathbf{H}}\cdot {\mathit{R}}_{\mathit{v}}-{\mathcal{D}}_{\mathrm{2}}\left({\mathit{\rho }}_{*}{R}_{wDp}\right).\end{array}\end{array}$$

The form of the Helmholtz problem is slightly unusual in the appearance of the first two terms rather than a single term $p^{\prime }/c^{\mathrm{2}}$. As noted above, the ${\stackrel{\mathrm{‾}}{{\stackrel{\mathrm{‾}}{p^{\prime }}}^w}}^p$ term appears because of the dependence of p^′ on ${\stackrel{\mathrm{‾}}{{\mathit{\alpha }}^{\prime }}}^p$ and of α^′ on ${\stackrel{\mathrm{‾}}{p^{\prime }}}^w$. Nevertheless, this slightly different form does not affect the stencil of the discrete Helmholtz operator or the difficulty of solving the Helmholtz equation numerically.

The vertical part of the Helmholtz operator must be modified near the top and bottom boundaries to impose $w^{\prime }=\mathrm{0}$ there. This modification amounts to

• omitting the contributions from levels $\mathrm{1}/\mathrm{2}$ and $N_z+\mathrm{1}/\mathrm{2}$ in 𝒟₂(ρ_*R_wDp) on the right-hand side of Eq. (74) and

• omitting the contributions 𝒟₁(p^′) from levels $\mathrm{1}/\mathrm{2}$ and $N_z+\mathrm{1}/\mathrm{2}$ in 𝒟₂𝒟₁(p^′).

The switching of the Lagrange multipliers used to enforce non-negativity of q^l and q^f can significantly change the linearization of the thermodynamics, particularly the value of N². Therefore, the matrix M and the coefficients of the Helmholtz problem are rebuilt at every solver iteration.

2.11 Helmholtz solver and back substitution

The Helmholtz problem is solved using a horizontal multigrid method. Each smoother iteration uses a Jacobi method in the horizontal with a tridiagonal direct vertical solve. Key parameters of the multigrid solver – the depth of the V cycles, the number of smoother iterations on the coarsest grid, and the number of V cycles – are automatically chosen to ensure that the pressure increments are sufficiently accurate (Appendix A3).

Having found the solution for p^′, increments to other variables are found through back substitution:

• u^′, v^′, and w^′ are found using Eqs. (34) and (70);

• velocity increments are used to compute mass flux increments F^′ and hence ρ^′ using Eq. (33);

• velocity and mass flux increments are used to update all transported terms and hence increment q and η; and

• α^′ is computed using Eq. (63), and then (${\mathit{\lambda }}^{\mathrm{l}}{}^{\prime },{\mathit{\lambda }}^{\mathrm{f}}{}^{\prime },q^{\mathrm{l}}{}^{\prime },q^{\mathrm{f}}{}^{\prime },{\mathit{\eta }}^{\mathrm{v}}{}^{\prime },{\mathit{\eta }}^{\mathrm{l}}{}^{\prime },{\mathit{\eta }}^{\mathrm{f}}{}^{\prime }$) are computed from q^′, η^′, and α^′ using the coefficients saved after the Gaussian elimination Eq. (55).

Before applying the increments (${\mathit{\lambda }}^{\mathrm{l}}{}^{\prime },{\mathit{\lambda }}^{\mathrm{f}}{}^{\prime },q^{\mathrm{l}}{}^{\prime },q^{\mathrm{f}}{}^{\prime },{\mathit{\eta }}^{\mathrm{v}}{}^{\prime },{\mathit{\eta }}^{\mathrm{l}}{}^{\prime },{\mathit{\eta }}^{\mathrm{f}}{}^{\prime }$), a check is made to see if they would imply breaking any of the constraints q^l≥0, q^f≥0, λ^l≤0, and λ^f≤0. If they would, then, at that grid point, a partial increment is made to these seven variables such that the updated variables sit at the constraint boundary, and the corresponding δⁱ is switched from 0 to 1 or vice versa. All other variables receive their full increments.

After the completion of the solver iterations the model's pressure field is diagnosed using Eqs. (18)–(20). This ensures that the pressure is consistent with the other model thermodynamic fields and enables bit-reproducible restarts without the need to save the pressure field. During early stages of model development, Eqs. (18)–(20) were also used to recompute p at the end of each solver iteration. However, this formulation is adversely impacted by the switching of constraints, as follows. The total entropy always receives its full increment η^′. When constraint switching leads to partial increments η^v′, η^l′, and η^f′, the dry entropy must take up the slack. Consequently, using Eqs. (18)–(20) to update the pressure can result in local pressure changes much larger than the p^′ returned by the Helmholtz problem, leading to large local residuals in the momentum equation. Therefore, in the present formulation, p is updated by adding the increment p^′ returned by the Helmholtz problem at each solver iteration. In the presence of constraint switching, this alternative way of updating the pressure greatly reduces the residuals in the momentum equation but leads to a larger residual $R_{\stackrel{\mathrm{‾}}{p}}$ in Eq. (17). Nevertheless, the overall effect is beneficial, accelerating solver convergence.

3 Test cases

3.1 Generating initial data

The model includes general-purpose routines for constructing a horizontally uniform initial state in discrete hydrostatic balance given specified vertical profiles of a temperature variable, a humidity variable, and the horizontal wind components. The temperature variable may be potential temperature, liquid water potential temperature, virtual potential temperature, or temperature itself. The humidity variable may be total specific humidity or relative humidity.

PTerodaC³TILES version 1.0 has routines in place, which can be selected via namelist options, to generate initial data and forcing terms for the following test cases: ATEX (Stevens et al., 2001), ARM (Brown et al., 2002), BOMEX (Siebesma et al., 2003), BUBBLE (a 3D version of Bryan and Fritsch, 2002), CBL (a dry convective boundary layer, Sullivan and Patton, 2011), DYCOMS (Stevens et al., 2005), and a NEUTRAL boundary layer with shear. The user can easily implement new test cases by using existing routines as templates.

PTerodaC³TILES version 1.0 can also generate initial data for the LBA case (Grabowski et al., 2006). However, the current version cannot represent the microphysics and precipitation necessary to simulate this case successfully, and the forcing terms are not yet implemented.

3.2 Output and diagnostics

A variety of diagnostics are computed online and may be selected by the user for output. These include

• time series of global diagnostics such as total mass and total water, along with accumulated source terms;

• column diagnostics of horizontal means of model fields and derived quantities such as turbulent kinetic energy, vertical fluxes, and cloud fraction;

• two-dimensional slices in the x–y, x–z, and y–z planes of the main model fields plus some derived quantities like cloud top height;

• diagnostics of quasi-Newton solver convergence; and

• diagnostics of quantities potentially related to model stability.

Full details are given in the user manual, available from Zenodo; see the “Code availability” section.

4 Verification

All components of the model have been thoroughly tested during development. This section highlights some aspects that require particularly careful checking or that can appear to work despite errors in formulation or coding.

4.1 Advection

The advection routines have been tested with specified velocity fields, independently of the rest of the model, to verify their overall accuracy, including correct behaviour at vertical and lateral boundaries and, where relevant, to confirm their conservation, consistency, and boundedness properties. It is particularly important that the w-level dual mass budget is satisfied (Sect. 2.2) so that w-level scalars are correctly advected. Conservation, consistency, and boundedness have also been verified in full model simulations to ensure that they are maintained by both the full cascade advection and the cheap advection updates.

4.2 Divergence consistency

It has been verified that the divergent velocity increments discussed in Sect. 2.6 correctly compensate for any discrepancy between the departure volumes computed by SLICE and the trajectory-average divergence. When these compensating increments are not included the model is found to become unstable, as expected theoretically.

4.3 Thermodynamics and linearized thermodynamics

The correctness and consistency of the formulation of thermodynamics in terms of internal energy potentials was verified in standalone code before incorporation in PTerodaC³TILES. The use of switched constraints and partial increments within the quasi-Newton solver allows q^l and q^f to be zero but not negative, as intended. For the PTerodaC³TILES implementation, the correctness of the various test case initial states, constructed to be in hydrostatic balance and thermodynamic equilibrium, and the broadly correct evolution of thermodynamic profiles and clouds provide further verification of the thermodynamics.

The correct linearization of the thermodynamics encoded in the matrix $\stackrel{\mathrm{̃}}{\mathbf{M}}$ is critical to the success of the formulation but is fiddly and susceptible to errors. The linearization has been verified by comparing the actual changes in the left-hand sides of Eqs. (6)–(13) with the changes predicted by the linearization when individual thermodynamic variables are perturbed. This testing considered base states with and without liquid and frozen water to ensure that all cases were covered.

4.4 Solver convergence

Good solver convergence is also critical for the success of the formulation. A useful rule of thumb is that the residuals should decrease by roughly an order of magnitude per solver iteration so that only a small number of iterations is needed per time step. However, ensuring correct performance of the solver is far from trivial, as it depends on correct formulation and implementation of the linearization, derivation of the Helmholtz problem, solution of the Helmholtz problem, and back substitution. Therefore, the ability of the solver to correct known errors was directly verified, as follows.

A known solution for the model state at step n+1 is required. This may be a known steady state such as a horizontally uniform state in hydrostatic balance or a more complex three-dimensional state obtained by taking a sufficiently large number of solver iterations. When a small perturbation is made to this state, the solver should be able to fix that perturbation almost completely in one iteration. Testing included perturbing each model state variable in turn and considered perturbations with different spatial structures, including globally uniform, horizontally uniform vertically localized, horizontally localized vertically uniform, and localized at a single point in the domain interior or at a top or bottom boundary.

Figure 5 shows the maximum residuals in the q and u equations and the maximum pressure increment versus iteration number at hour 9 of the ARM case when both the dynamics and moist thermodynamics are very active. For the purpose of calculating these diagnostics, the model was restarted for a single time step, and 12 solver iterations were taken rather than the default 3. There is a very large drop in the maximum q residual between the first and second iterations associated with the cheap transport updates (Sect. 2.6). There is a rather small reduction in the pressure increment between the first and second iterations, probably related to constraint switching in the thermodynamics. Otherwise, all three quantities steadily decrease by nearly an order of magnitude per iteration. Note that the residuals and pressure increments calculated at the fourth iteration are a measure of the errors due to incomplete solver convergence that remain upon completion of the third iteration.

Figure 5 Maximum absolute value within the model domain of (a) q equation residual (dimensionless), (b) u equation residual (m s⁻¹), and (c) p^′ (Pa) versus iteration number. These diagnostics were computed for a single time step, restarting from hour 9 of the ARM case and taking 12 solver iterations.

As discussed in Sect. 2.11, when the thermodynamic state switches between absence and presence of condensate, a subset of thermodynamic variables local to the switch receive only partial increments. It is important to understand the extent to which such partial increments adversely affect the solver convergence. The DYCOMS case is especially challenging in this regard. If the cloud deck begins to break up, then holes form in the cloud (e.g. Fig. 12). If the horizontal advective Courant number is greater than 1, then a very large number of grid cells per step switch between absence and presence of liquid water. Moreover, the extremely large gradients in humidity and entropy at cloud top, which get folded into the cloud holes, make the thermodynamic state at the cloud top and cloud hole edges very sensitive to the solver advection updates, so that switching can occur at the second and subsequent solver iterations.

An initial attempt to run the DYCOMS case at $\mathrm{128}\times \mathrm{128}\times \mathrm{300}$ resolution ($\mathrm{\Delta }x=\mathrm{\Delta }y=\mathrm{25}$ m, Δz=5 m) with N_ℓ=3 solver iterations resulted in failure of the model after about 32 min. Diagnostics indicated that the switching had not completely settled down after three solver iterations. Nevertheless, the number of switches decreases by 2 orders of magnitude per iteration, and increasing the number of iterations to N_ℓ=4 is sufficient for virtually all of the switching to settle down (Table 1), allowing the DYCOMS case to run successfully. Although constraint switching can slow solver convergence, no cases have been encountered in which constraint switching prevents solver convergence.

Table 1 Representative number of grid points at which switches occurred at different solver iterations. The numbers were monitored over 10 time steps 1 h into the DYCOMS case. The total number of model grid points was $\mathrm{128}\times \mathrm{128}\times \mathrm{300}$.

4.5 Stability

The semi-implicit semi-Lagrangian formulation of PTerodaC³TILES should be stable for large acoustic, gravity-wave, and advective Courant numbers. Diagnostics optionally output by the model confirm that the model does run stably with advective Courant numbers greater than 1 and with very large acoustic Courant numbers. For example, in the ARM case presented in Sect. 5.2 the horizontal and vertical acoustic Courant numbers remain around 53 and 87, respectively. The horizontal advective Courant number in the x direction varies in the range 1.5 to 2, while the vertical advective Courant number peaks at about 3. Both the gravity-wave and convective Courant numbers $\sqrt{\max \left(N^{\mathrm{2}}\right)}\mathrm{\Delta }t$ and $\sqrt{\max (-N^{\mathrm{2}})}\mathrm{\Delta }t$, respectively, peak at around 0.6. For the DYCOMS case presented in Sect. 5.3, the horizontal and vertical acoustic Courant numbers remain around 49 and 344, respectively. The horizontal advective Courant numbers in both the x and y directions are greater than 1, while the vertical advective Courant number peaks at over 3. The gravity-wave Courant number is about 0.9 while the convective Courant number varies between 0.8 and 1.4.

Because the ILES formulation does not include an eddy-diffusion-based subgrid scheme, the diffusive Courant number is not relevant for model stability.

In all cases the model does become unstable if the time step is chosen too large. Further careful study is needed to understand exactly what limits the model stability. However, based on the theoretical properties of the methods and experience to date, likely candidates for instability mechanisms include the following.

1. When the deformational Courant number, or, more precisely, the maximum magnitude of the eigenvalues of bΔt∇u, approaches 1, the semi-Lagrangian trajectory calculations become increasingly inaccurate. This might lead to badly distorted departure cells such that the simple small-amplitude advecting velocity correction cannot restore divergence consistency (Sect. 2.6), or there might be other feedbacks via advection that amplify errors. Diagnostics of the components of Δt∇u for the cases presented in Sect. 5 suggest that the model is close to that regime in the surface layer where vertical shear is strongest.

2. If a²N²Δt² approaches −1, then the factor 1+a²N²Δt² in the denominator of Eqs. (70) and (71) approaches zero and those coefficients blow up. This is a known limitation of the semi-implicit treatment of gravity waves when the stratification becomes unstable (Davies et al., 2005); it can be mitigated by bounding the value of N² that appears in those coefficients. In PTerodaC³TILES no such bound is applied.

3. Any factors that inhibit convergence of the quasi-Newton solver can mean that the theoretical stability of a semi-implicit, semi-Lagrangian scheme is not attained. This could occur, for example,

   • if the change in state over one time step is too large, hence too nonlinear, to be captured by the quasi-Newton linearization;

   • if insufficient iterations are taken for constraint switching to settle down (as discussed in Sect. 4.4); and

   • if terms omitted from the linearization, such as $u\partial \mathit{\eta }/\partial x$, become important.

   (Coriolis terms are omitted from the linearization Eqs. (34) and (35), but $\left|\mathbf{\Omega }\right|\mathrm{\Delta }t$ is unlikely to be large enough, in practice, for this to matter.)

It would be a valuable addition to the model to include mitigation strategies, such as adaptive time stepping, if the mechanisms limiting the model stability could be better understood and quantified.

4.6 Conservation

PTerodaC³TILES is designed to have closed budgets for mass, water, and entropy. This property requires conservative advection of these quantities by both SLICE and the cheap transport updates (Sect. 2.6), as well as careful formulation of surface and interior sources (Sect. 2.7), all taking into account the w-level dual mass budget (Sect. 2.2).

Figure 6 shows time series of various budget quantities for the ARM case. Figure 6a shows that the change in total mass (blue curve) agrees with the accumulated source of total mass, in this case due to the surface moisture flux (blue symbols), while the dry mass is exactly conserved. Figure 6b shows that the change in total water mass (black curve) agrees with the accumulated source of water (black symbols), even while some water condenses into liquid (red curve). Figure 6d shows that the change in total entropy (black curve) agrees with the accumulated source of entropy (black symbols). Figure 6c confirms that the momentum budget is not closed because the semi-Lagrangian advection of momentum is not conservative. Figure 6e shows that the change in total energy (black curve), which averages 432 W m⁻², is slightly less than the accumulated total energy source (black symbols). The difference, shown in Fig. 6f, is due to numerical dissipation of kinetic energy and mixing of water and entropy. This energy loss averages about 12 W m⁻².

Figure 6 Global changes and accumulated sources versus time for key budget quantities for the ARM case. (a) Total mass and dry mass (kg m⁻²). (b) Total water; 100× global liquid water and frozen water (kg m⁻²). (c) u momentum and v momentum (kg m s⁻¹ m⁻²). (d) Entropy (J K⁻¹ m⁻²). (e) Energy (J m⁻²). (f) Numerical energy change, i.e. the actual energy change minus the change due to explicit sources (J m⁻²).

Dissipation of kinetic energy and mixing of constituents and heat should result in a source of entropy that exactly compensates for this energy loss; however, this entropy source is currently neglected. See further discussion in Sect. 6. Because the current formulation assumes thermodynamic equilibrium in each grid cell, entropy sources due to departures from equilibrium (the J^TP term in Eq. 4) are zero.

4.7 Bit reproducibility

PTerodaC³TILES produces bit-reproducible results when restarted from a checkpoint file³. This property is invaluable for development and testing, as well as when rerunning sections of a simulation to obtain additional diagnostics. PTerodaC³TILES also produces bit-identical answers whether run with or without OpenMP shared memory parallelism, verifying the correctness of the parallel implementation.

5 Evaluation

This section presents results from some standard LES test cases to demonstrate the performance of PTerodaC³TILES version 1.0⁴. The same spatial resolution is used as in the original intercomparison articles defining the test case specifications. Although relatively coarse by current-day standards, this facilitates comparison with those published results and helps to highlight any limitations of PTerodaC³TILES version 1.0 that might be less conspicuous at finer resolution.

All results shown in this section used SLICE advection with piecewise constant remapping for density and parabolic spline remapping with a limiter and the Charney–Phillips grid correction for water and entropy. Semi-Lagrangian advection of velocity components used no limiter and the freeslip option for extrapolation near the bottom boundary.

5.1 BOMEX

The BOMEX test case (Siebesma et al., 2003) is based on observations made during the Barbados Oceanographic and Meteorological Experiment. It simulates a scenario of shallow cumulus over the ocean in which large-scale forcing, radiation, and turbulent and convective fluxes maintain a quasi-steady balance.

As in Siebesma et al. (2003), PTerodaC³TILES used a $\mathrm{64}\times \mathrm{64}\times \mathrm{75}$ grid with $\mathrm{\Delta }x=\mathrm{\Delta }y=\mathrm{100}$ m, Δz=40 m. The time step was Δt=10 s and the simulation was run for 6 h.

Figure 7 shows time series of three key quantities: total cloud cover, liquid water path (LWP), and turbulent kinetic energy (TKE). All three time series agree broadly with Fig. 2 of Siebesma et al. (2003). After an initial spin-up during the first hour or so, the total cloud cover and LWP fluctuate but have little trend, while the TKE continues to grow slowly. The PTerodaC³TILES cloud cover is slightly lower and the TKE slightly larger than the ensemble means in Siebesma et al. (2003) but within the typical inter-model spread.

Figure 7 Time series of (a) total cloud cover, (b) liquid water path, and (c) turbulent kinetic energy for the BOMEX case.

Figure 8 shows several horizontally averaged profiles from the BOMEX case. Figure 8b–d agree well with the corresponding figures from Siebesma et al. (2003) (their Figs. 6, 3d, and 4a). Figure 8a, e, and f also broadly agree with the corresponding figures from Siebesma et al. (2003) (their Figs. 3c, 4e, and 5a). However, an excessively strong shear layer has formed between model levels 1 and 2, consistent with the idea that the ILES approach poorly represents vertical subgrid transports near a horizontal boundary. This excessively strong shear layer leads to noise in the lowest two to three levels in the profiles of u, momentum flux, and TKE.

Figure 8 Vertical profiles of horizontally averaged quantities from the BOMEX case. (a) u and v (m s⁻¹). (b) Cloud fraction (dimensionless). (c) Liquid water q^l (dimensionless). (d) Vertical eddy flux of total water $\langle w^{\prime }{q}^{\prime }\rangle$ (m s⁻¹). (e) Vertical eddy fluxes of momentum $\langle w^{\prime }{u}^{\prime }\rangle$ and $\langle w^{\prime }{v}^{\prime }\rangle$ (m² s⁻²). (f) TKE (m² s⁻²). Panels (a)–(c) are averaged over the last hour of simulation; panels (d)–(f) are averaged over the last 3 h.

On a closely related point, numerical experimentation revealed that the momentum budget and the boundary layer u profile are strongly sensitive to the details of the semi-Lagrangian interpolation scheme for velocity components near the bottom boundary (Sect. 2.6). For example, switching on the limiter for velocity advection or using the noslip$\bigsqcup \bigsqcup$ bottom boundary extrapolation option resulted in a large spurious numerical source of eastward momentum and a significant shift to the right of the boundary layer u profile.

5.2 ARM

The ARM test case (Brown et al., 2002) simulates the diurnal evolution of shallow convection over land, starting from an initially stable and cloud-free boundary layer. The number, size, and depth of clouds evolve during the day, and clouds often overshoot their level of neutral buoyancy.

As in Brown et al. (2002), PTerodaC³TILES used a $\mathrm{96}\times \mathrm{96}\times \mathrm{110}$ grid with $\mathrm{\Delta }x=\mathrm{\Delta }y=\mathrm{66.7}$ m, Δz=40 m. The time step was Δt=10 s and the simulation was run for 14.5 h. Large-scale forcing terms were omitted, since they have only a small effect on the simulation (Brown et al., 2002).

Figure 9 shows a time–height plot of cloud fraction and TKE. The evolution of the cloud fraction, which peaks at a little over 0.15 around 6 or 7 h, as well as the height of cloud base and cloud top, agrees well with Fig. 5 of Brown et al. (2002). The intensity of TKE in the boundary layer grows and then decays in concert with the strength of the surface heat flux, and there is a clear signature of the formation of gravity waves in and above the cloud layer.

Figure 9 TKE (shading) and cloud fraction (contours) versus time and height for the ARM case. The contour values are 0.0001, 0.05, 0.1, and 0.15. The peak TKE value is 1.74 m² s⁻².

Figure 10 shows profiles of u and w variances at 3 and 9 h. The peak in w variance after 3 h is somewhat smaller than most ensemble members in Brown et al. (2002) (their Fig. 6), but otherwise these profiles agree well with Brown et al. (2002), including the secondary peak in w variance in the cloud layer at 9 h and the small peak in u variance near the boundary layer top at 3 h, which disappears after clouds form (see discussion in Brown et al., 2002, Sect. 3b).

Figure 10 Vertical profiles of (a) w variance and (b) u variance for the ARM case at 3 and 9 h. The profiles are also averaged in time between minus and plus 30 min of the nominal diagnostic time.

Figure 11 shows snapshots of cloud top height at selected times. The behaviour is consistent with that documented for the ARM case, with many small and shallow clouds at early times, gradually growing in depth, evolving towards fewer, larger, and deeper clouds with smaller total cloud cover at later times. The horizontal motion of the higher cloud tops is almost exactly in the x direction (the direction of the background geostrophic wind), while near cloud base there is a small component of motion in the y direction, leading to a distinct characteristic tilt to the clouds.

Figure 11 Snapshots of cloud top height (grey) for the ARM case at 6, 8, and 10 h. For this purpose, grid cells with $q^{\mathrm{l}}+q^{\mathrm{f}}>{\mathrm{10}}^{-\mathrm{5}}$ are defined to be cloudy (though q^f≡0 for ARM). The green background indicates no cloud in that column, and the contour interval is 100 m; refer to Fig. 9 for the lowest cloud base and highest cloud top at these times.

5.3 DYCOMS

The DYCOMS test case (Stevens et al., 2005) simulates a very different boundary layer regime from BOMEX and ARM: a nocturnal stratocumulus cloud layer over the ocean. The cloud layer is capped by a very strong and sharp inversion, with q decreasing by a factor of 6 and a jump in potential temperature of more than 9 K over 5 m, the recommended vertical grid spacing. Observations suggest that the cloud cover should be maintained close to 100 %. However, LES often fail to maintain the cloud cover, since excessive mixing across the inversion can lead to evaporative cooling driving cloud-free downdraughts.

Figure 12 Snapshots of cloud top height (pale grey) for the DYCOMS case at 1, 2, and 3 h. The blue background indicates no cloud in that column, and the contour interval is 50 m; refer to Fig. 14 for the lowest cloud base and highest cloud top at these times.

The DYCOMS case is expected to be particularly testing for the PTerodaC³TILES formulation, for several reasons. First, there are no physical parameters that can be adjusted to control the strength of mixing near the inversion, since the mixing is entirely associated with the numerics. There is, however, some sensitivity to the choice of numerical options and parameters, as discussed below. Second, good convergence of the semi-implicit iterative solver depends on having a sufficiently good linearization. The presence of near discontinuities in the distributions of total entropy and total specific humidity mean that the linearized advection terms Eqs. (41) and (42) are necessarily less accurate. The solver convergence is also affected by the switching of thermodynamic constraints (Sects. 2.11 and 4.4). Once cloud-free downdraughts form (Fig. 12) and are advected across the grid with horizontal Courant numbers greater than 1, condensate appears or disappears at large numbers of grid points every step. Moreover, the large gradients in entropy and specific humidity mean that transport increments resulting from small velocity increments at one solver iteration are enough to cause constraint switching at the next solver iteration. Experience to date suggests that more solver iterations per step are required for DYCOMS than for other test cases in order for the constraint switching to settle down (Sect. 4.4).

As in Stevens et al. (2005), PTerodaC³TILES used a $\mathrm{96}\times \mathrm{96}\times \mathrm{300}$ grid with $\mathrm{\Delta }x=\mathrm{\Delta }y=\mathrm{35}\mathrm{m}$, Δz=5m. The time step was Δt=5s and the simulation was run for 4h.

Figure 13 (compare Stevens et al., 2005, Fig. 2) shows that in the PTerodaC³TILES simulation cloud cover gradually falls to about 40 % while the liquid water path falls to around 0.02 kg m⁻² (about 0.035 kg m⁻² for the ensemble mean in Stevens et al., 2005). The TKE settles down at about 200 kg s⁻², about half of the ensemble mean in Stevens et al. (2005) (noting the different units used). The gradual reduction in cloud cover is also clearly illustrated in Fig. 14. The figure confirms that turbulence grows initially both from the surface (due to the surface sensible heat flux) and from cloud top (due to cloud top radiative cooling).

Figure 13 Time series of (a) total cloud cover, (b) liquid water path, and (c) turbulent kinetic energy for the DYCOMS case.

Figure 14 TKE (shading) and cloud fraction (contours) versus time and height for the DYCOMS case. The contour values are 0.3 to 0.9 in steps of 0.1 and 0.99. The peak TKE value is 1.16 m² s⁻².

The PTerodaC³TILES simulation of cloud in the DYCOMS case appears to be quite sensitive to the choice of numerical options and parameters. For example, when the time step is reduced from 5 to 2 s, the cloud breakup is slowed, leaving 70 % cloud cover after 4 h. On the other hand, when parabolic spline remapping is replaced by piecewise parabolic method remapping for the advection of entropy and water, almost all cloud disappears after 4 h. Further experiments and diagnostics to understand these sensitivities would be valuable.

5.4 Computational cost

While noting the usual caveats that computational costs are sensitive to computing platform, compiler, output files written, and many other factors, it is nevertheless useful to give would-be users an idea of the computational cost of PTerodaC³TILES. Table 2 gives the wall-clock times for the three test cases discussed in this section. The cases were run on a 2022 MacBook Pro with an 8 core M2 chip. The gfortran compiler was used with the −O₃ optimization flag and OpenMP for shared memory parallelism. The parallel speed-up ranged from 574 % for BOMEX to 652 % for DYCOMS.

Table 2 Computational cost of PTerodaC³TILES for three standard test cases on a 2022 MacBook Pro laptop.

The multigrid elliptic solver and the full transport calculation dominate the cost, accounting for approximately 40 % and 35 % of the total cost, respectively, with some variation from case to case. Creating the thermodynamic subsystem matrix and carrying out the Gaussian elimination account for less than 15 % of the total cost.

The cost of the model scales linearly with the number of time steps and hence almost inversely with the time step size Δt. The scaling is not exactly inverse because the multigrid scheme parameters take account of the acoustic Courant number (Appendix A3). The scaling of the cost with spatial resolution is affected by memory usage as well as floating point operations, so it will be strongly system dependent.

6 Conclusions and discussion

It has been demonstrated that a consistent treatment of moist thermodynamics, expressed in terms of internal energy potentials, can be incorporated, without excessive computational expense, in a three-dimensional computational fluid dynamics code suitable for the study of atmospheric boundary layers and shallow convection. In the current implementation the moist thermodynamics is fully incorporated within the dynamical core rather than treated as separate physics source terms.

The iterative solver for the semi-implicit time integration scheme requires a linearization of the thermodynamics, with elimination of unknowns to leave a typical Helmholtz problem for the pressure increment. The moderate sparsity and fixed sparsity pattern of the Jacobian matrix of the (w-level) thermodynamic subsystem are exploited to reduce the cost of the elimination.

The use of internal energy potentials has several advantages over the Gibbs function approach used by Thuburn (2017b) (see Sect. 1). However, one disadvantage of the internal energy potential approach is that it does not permit seamlessly switching to a (quasi-)Boussinesq equation of state since e_α terms (used throughout the algorithm) would become undefined. A nearly Boussinesq fluid could be simulated by making the sound speed extremely large. Alternatively, if a Boussinesq option is a requirement, then the above difficulty could be avoided by using specific enthalpies as the potentials in formulating the thermodynamics (Chris Eldred, personal communication, 2024).

The numerical methods used by PTerodaC³TILES are more typical of those used in global weather and climate models than traditional LES models. The semi-implicit semi-Lagrangian scheme permits time steps significantly larger than are commonly used in traditional LES models at similar resolution. Experience to date, consistent with the theoretical properties of the numerical methods, suggests that the deformational Courant number is most often the factor limiting the maximum stable time step. However, there is no simple, easily monitored stability criterion (computing the eigenvalues of Δt∇u at every grid point and every step would be expensive), so a degree of experimentation is required to find a suitable time step for each simulation. An important caveat is that stability does not imply accuracy. For example, the DYCOMS case (Sect. 5.3) shows significant sensitivity to the size of the time step. Further investigation is required to quantify the extent to which accuracy declines as we push the time step towards the stability limit and which flow regimes, like DYCOMS, are especially sensitive to the time step.

For the BOMEX, ARM, and DYCOMS cases presented here, and other cases not discussed, PTerodaC³TILES can produce plausible simulations of standard LES test cases, comparable to other model results in the literature. These results are encouraging for the ILES approach and, more widely, for the use of global models at sub-kilometre resolution. At the same time, these results have given some useful initial indications of the limitations of the ILES approach, as well as highlighting areas in need of further investigation.

Most notably, ILES produces weak vertical fluxes of momentum and scalars near the bottom boundary. These weak fluxes lead to excessive vertical gradients, which, in turn, result in further errors. For example, in BOMEX the spuriously strong vertical shear exacerbates the conservation errors of the cubic semi-Lagrangian advection of momentum and makes the results very sensitive to the details of the semi-Lagrangian interpolation in the lowest levels. In an ARM simulation with finer vertical resolution near the surface (not shown), the spuriously weak mixing near the surface allows a thin layer of fog to form briefly during the first hour. Some initial experimentation showed benefits in distributing the convergence of the surface momentum flux over several layers near the surface rather than just a single layer. In future work it is planned to investigate this idea further, applied to scalars as well as momentum, and to try to determine how the optimal distribution depends on the flow regime as well as numerical factors such as the grid resolution and anisotropy.

Another broad area in need of further investigation is the sensitivity to numerical methods and parameters. The sensitivity of the DYCOMS case to the time step and the advection remapping scheme has been mentioned already. Other potentially significant factors include the grid isotropy, the use of different advection options such as the use of limiters and the Charney–Phillips grid correction (Eq. 27), and the modifications to momentum interpolation near the bottom boundary.

Taking inspiration from global models that predict potential temperature, PTerodaC³TILES has a closed budget for entropy rather than energy. The entropy source that should be associated with numerical mixing of scalars and dissipation of kinetic energy is neglected, resulting in a small but systematic energy loss (Fig. 6). If desired, the global numerical energy loss could easily be diagnosed and returned as an entropy source to close the energy budget. Initial attempts to diagnose the local numerical energy loss (not shown) suggest that the calculation is subtle and far from trivial and might not even be well defined, in large part because of the Charney–Phillips vertical grid staggering. Further work will be needed to diagnose the local numerical energy loss, if, indeed, it is possible at all. A strong motivation to continue these attempts is that such an estimate of energy dissipation would be a key input to a stochastic parameterization of backscatter (e.g. Mason and Thomson, 1992; Brown et al., 1994), which would be a useful extension to the model's functionality.

Two closely related priorities for future work are the inclusion of a simple microphysics scheme with precipitation and the inclusion of thermodynamic nonequilibrium effects. These developments will enhance the capabilities of PTerodaC³TILES, allowing it to be applied to a wider range of cases. Equally importantly, they will test whether the thermodynamic potential approach can be applied straightforwardly and efficiently to more complex physical processes, beyond the coupling of the equation of state to the dynamical core.

Appendix A

A1 Stretched vertical grid

When a stretched vertical grid is chosen, the height of the model level with index k is given by

$$\begin{array}{}\text{(A1)}& z_k=a_{\left(z\right)}\log \left(\mathrm{1}+d_{\left(z\right)}{b}_{\left(z\right)}^k\right)+c_{\left(z\right)},\end{array}$$

where k is an integer for p levels and an integer plus $\mathrm{1}/\mathrm{2}$ for w levels. The parameters in Eq. (A1) are set so that the grid spacing approaches a stretching factor b_(z) per level for small k and a uniform grid spacing u_(z) for large k, with b_(z) and u_(z) specified by the user when creating initial data. In terms of N_z, the number of p levels, and $D_z=z_{N_z+\mathrm{1}/\mathrm{2}}-z_{\mathrm{1}/\mathrm{2}}$, the domain depth with $z_{\mathrm{1}/\mathrm{2}}=\mathrm{0}$ the height of the bottom boundary, the parameters are given by

$$\begin{array}{}\text{(A2)}& {\displaystyle }& {\displaystyle }{a}_{\left(z\right)}=u_{\left(z\right)}/\log b_{\left(z\right)},\text{(A3)}& {\displaystyle }& {\displaystyle }{d}_{\left(z\right)}={\displaystyle \frac{\mathrm{1}-\exp (D_z/a_{\left(z\right)})}{b_{\left(z\right)}^{\mathrm{1}/\mathrm{2}}\left(b_{\left(z\right)}^{N_z}-\exp \left(D_z/a_{\left(z\right)}\right)\right)}},\text{(A4)}& {\displaystyle }& {\displaystyle }{c}_{\left(z\right)}=D_z-a_{\left(z\right)}\log \left(\mathrm{1}+d_{\left(z\right)}{b}_{\left(z\right)}^{N_z+\mathrm{1}/\mathrm{2}}\right).\end{array}$$

Importantly, Eq. (A1) has the virtue of being invertible, allowing the cell index k to be determined for any height z without the need for expensive searching in the semi-Lagrangian departure-point calculations.

A2 Expressions for internal energy and related quantities

The expressions used for specific internal energy of dry air, water vapour, liquid water, and frozen water are

$$\begin{array}{}\text{(A5)}& {\displaystyle }& {\displaystyle }{e}^{\mathrm{d}}\left({\mathit{\alpha }}^{\mathrm{d}},{\mathit{\eta }}^{\mathrm{d}}\right)=C_{\mathrm{v}}^{\mathrm{d}}{T}_{\mathrm{0}}\exp \left({\displaystyle \frac{{\mathit{\eta }}^{\mathrm{d}}-C_{\mathrm{p}}^{\mathrm{d}}}{C_{\mathrm{v}}^{\mathrm{d}}}}\right){\left({\displaystyle \frac{{\mathit{\alpha }}_{\mathrm{0}}^{\mathrm{d}}}{{\mathit{\alpha }}^{\mathrm{d}}}}\right)}^{R^{\mathrm{d}}/C_{\mathrm{v}}^{\mathrm{d}}},\text{(A6)}& {\displaystyle }& {\displaystyle }{e}^{\mathrm{v}}\left({\mathit{\alpha }}^{\mathrm{v}},{\mathit{\eta }}^{\mathrm{v}}\right)=C_{\mathrm{v}}^{\mathrm{v}}{T}_{\mathrm{0}}\exp \left({\displaystyle \frac{{\mathit{\eta }}^{\mathrm{v}}-C_{\mathrm{p}}^{\mathrm{v}}}{C_{\mathrm{v}}^{\mathrm{v}}}}\right){\left({\displaystyle \frac{{\mathit{\alpha }}_{\mathrm{0}}^{\mathrm{v}}}{{\mathit{\alpha }}^{\mathrm{v}}}}\right)}^{R^{\mathrm{v}}/C_{\mathrm{v}}^{\mathrm{v}}}+L_{\mathrm{00}}^{\mathrm{s}},\end{array}$$

$$\begin{array}{}\text{(A7)}& {\displaystyle }& {\displaystyle }{e}^{\mathrm{l}}\left({\mathit{\eta }}^{\mathrm{l}}\right)=C^{\mathrm{l}}{T}_{\mathrm{0}}\exp \left({\displaystyle \frac{{\mathit{\eta }}^{\mathrm{l}}-C^{\mathrm{l}}}{C^{\mathrm{l}}}}+{\displaystyle \frac{L_{\mathrm{00}}^{\mathrm{v}}-{\mathit{\alpha }}^{\mathrm{l}}{p}_{\mathrm{0}}^{\mathrm{sat}}}{C^{\mathrm{l}}{T}_{\mathrm{0}}}}\right)+L_{\mathrm{00}}^{\mathrm{f}},\text{(A8)}& {\displaystyle }& {\displaystyle }{e}^{\mathrm{f}}\left({\mathit{\eta }}^{\mathrm{f}}\right)=C^{\mathrm{f}}{T}_{\mathrm{0}}\exp \left({\displaystyle \frac{{\mathit{\eta }}^{\mathrm{f}}-C^{\mathrm{f}}}{C^{\mathrm{f}}}}+{\displaystyle \frac{L_{\mathrm{00}}^{\mathrm{s}}-{\mathit{\alpha }}^{\mathrm{f}}{p}_{\mathrm{0}}^{\mathrm{sat}}}{C^{\mathrm{f}}{T}_{\mathrm{0}}}}\right).\end{array}$$

These expressions differ from those used by Bowen and Thuburn (2022a) in that the constant $L_{\mathrm{00}}^{\mathrm{s}}/T_{\mathrm{0}}$ has been subtracted from the specific entropy of all water phases.

For $i=\mathrm{d},\mathrm{v},\mathrm{l},\mathrm{f}$, the species temperature is given by

$$\begin{array}{}\text{(A9)}& T^i=e_{\mathit{\eta }}^i.\end{array}$$

For $i=\mathrm{d},\mathrm{v}$, the species partial pressure is given by

$$\begin{array}{}\text{(A10)}& p^i=-e_{\mathit{\alpha }}^i.\end{array}$$

The Gibbs function for water vapour is given by

$$\begin{array}{}\text{(A11)}& g^{\mathrm{v}}=e^{\mathrm{v}}+{\mathit{\alpha }}^{\mathrm{v}}{p}^{\mathrm{v}}-{\mathit{\eta }}^{\mathrm{v}}{T}^{\mathrm{v}},\end{array}$$

while for $i=\mathrm{l},\mathrm{f}$,

$$\begin{array}{}\text{(A12)}& g^i=e^i+{\mathit{\alpha }}^{i}p-{\mathit{\eta }}^i{T}^i,\end{array}$$

with $p=p^{\mathrm{d}}+p^{\mathrm{v}}$. (The Gibbs function for dry air is not needed.)

A3 Parameters for multigrid scheme

The tuning of parameters in numerical solvers for elliptic problems is often empirical and case-dependent. Here, both the Helmholtz problem and the properties of the multigrid solver are well understood, which allows some key parameters in the multigrid scheme to be set automatically.

Let $C_{\mathrm{ac}}^x=c\mathrm{\Delta }t/\mathrm{\Delta }x$ and $C_{\mathrm{ac}}^y=c\mathrm{\Delta }t/\mathrm{\Delta }y$ be the horizontal acoustic wave Courant numbers in the x and y directions, respectively, appropriate to whichever grid in the hierarchy is under consideration. Let $C_{\mathrm{ac}}=\max (C_{\mathrm{ac}}^x,C_{\mathrm{ac}}^y)$. In the following calculations the temporal off-centring parameter a is effectively set to 1.

Each smoother iteration uses a direct solve in the vertical direction; thus, we must consider the action of the smoother on errors of different horizontal scales. The error amplification factor for a uniform error for one Jacobi smoother iteration is

$$\begin{array}{}\text{(A13)}& A_{\mathit{\epsilon }}=\mathrm{1}-{\displaystyle \frac{\mathit{\mu }}{\mathrm{1}+\mathrm{2}\left({\left(C_{\mathrm{ac}}^x\right)}^{\mathrm{2}}+{\left(C_{\mathrm{ac}}^y\right)}^{\mathrm{2}}\right)}},\end{array}$$

where μ is the under-relaxation parameter, while the error amplification factor for a checkerboard pattern error is

$$\begin{array}{}\text{(A14)}& A_{\mathit{\epsilon }}=\mathrm{1}-\mathit{\mu }-{\displaystyle \frac{\mathrm{2}\mathit{\mu }\left({\left(C_{\mathrm{ac}}^x\right)}^{\mathrm{2}}+{\left(C_{\mathrm{ac}}^y\right)}^{\mathrm{2}}\right)}{\mathrm{1}+\mathrm{2}\left({\left(C_{\mathrm{ac}}^x\right)}^{\mathrm{2}}+{\left(C_{\mathrm{ac}}^y\right)}^{\mathrm{2}}\right)}}.\end{array}$$

1. Under-relaxation parameter. The under-relaxation parameter is set to μ=0.8. This provides an appropriate compromise between the damping of large-scale errors on the coarsest grid, where C_ac≈1 implies A_ε≈0.85, and the damping of grid-scale errors on the finest grid, where C_ac≫1 implies A_ε≈0.6.

2. Depth of V cycles. As the horizontal grid is coarsened, the horizontal acoustic wave Courant number C_ac decreases, the horizontal part of the Helmholtz problem becomes more diagonally dominant, and the smoother iterations damp the error more and more quickly. Once the horizontal acoustic Courant number reaches about 1, there is little to be gained by coarsening the grid further. Thus, the desired V-cycle depth is set to the smallest depth needed to reduce the horizontal acoustic Courant number below 1. The maximum possible V-cycle depth permitted by the number of grid points in the x and y directions is also computed. The actual V-cycle depth is set equal to the minimum of the desired depth and the maximum permitted depth. If the maximum permitted depth is smaller than the desired depth, then a warning message printed and the scheme attempts to compensate by increasing the number of coarsest-grid smoother iterations.

3. Number of coarsest-grid smoother iterations. The rate at which large-scale errors are damped is determined by the damping of uniform errors on the coarsest grid. Using Eq. (A13) to estimate the large-scale error damping rate, the number of smoother iterations on the coarsest grid is chosen so that a uniform error is damped by a factor $\mathrm{1}/\mathrm{10}$ on each V cycle.

4. Number of V cycles. To ensure mass conservation, the back substitution computes density increments from the divergence of mass flux increments, with velocity increments computed from the gradient of pressure increments. As a consequence of the gradient and divergence operations, any grid-scale errors in the pressure increments due to imperfect convergence of the multigrid solver get amplified by a factor equal to the (fine-grid) horizontal acoustic Courant number squared. Thus, to ensure good convergence of the quasi-Newton solver, any grid-scale errors in the pressure increments computed by the multigrid solver must be much smaller, by a factor an order of magnitude smaller than $\mathrm{1}/C_{\mathrm{ac}}^{\mathrm{2}}$, than the pressure increments themselves. Using the estimate from Eq. (A14) that $\left|A_{\mathit{\epsilon }}\right|\approx |\mathrm{1}-\mathrm{2}\mathit{\mu }|$ for large C_ac, we can estimate the total number of fine-grid smoother iterations required and hence the number of V cycles.

A4 Alternative formulation of Monin–Obukhov theory for surface momentum flux

Monin–Obukhov similarity theory gives the flow speed U(z) at height z in terms of the friction velocity U_* (e.g. Stull, 1988):

$$\begin{array}{}\text{(A15)}& U\left(z\right)={\displaystyle \frac{U_{*}}{\mathit{\kappa }}}\left\{\log \left(z/z_{\mathrm{0}}\right)+\mathrm{\Psi }\left(\mathit{\zeta },{\mathit{\zeta }}_{\mathrm{0}}\right)\right\},\end{array}$$

where κ is the von Kármán constant; $\mathit{\zeta }=z/L$ is a non-dimensional height; ${\mathit{\zeta }}_{\mathrm{0}}=z_{\mathrm{0}}/L$, where z₀ is the surface roughness length; $L=-U_{*}^{\mathrm{3}}/F^{\mathrm{b}}$ is the Obukhov length with F^b the surface buoyancy flux; and Ψ is an integrated stability function, discussed below. As written, Eq. (A15) must be solved iteratively to obtain U_*, given U(z) at some z. Moreover, in the stable case (F^b<0) with light winds Monin–Obukhov similarity theory is valid only for sufficiently small z and may be outside its range of validity at the height of the lowest model level; in that case, with commonly used expressions for Ψ, the solution for U_* might not be unique or might not exist (e.g. Fig. A1).

Figure A1 Dimensionless friction velocity ${\widehat{U}}_{*}={\mathit{\zeta }}^{\mathrm{1}/\mathrm{3}}$ versus dimensionless wind speed $\widehat{U}=\mathrm{1}/\widehat{s}$ at height z in the stable regime given by Monin–Obukhov theory with the Businger–Dyer stability function (Dyer, 1974, asterisks), with the Cheng and Brutsaert (2005) stability function (open circles), with the modified Cheng and Brutsaert stability function Eq. (A26) (filled circles), and with the fit Eq. (A29) (thick curve) (a) for $r=z_{\mathrm{0}}/z=\mathrm{0.001}$ and (b) for $r=z_{\mathrm{0}}/z=\mathrm{0.1}$. Note that with the Businger–Dyer stability function there may be zero or two solutions for ${\widehat{U}}_{*}$, while with both the original and modified Cheng and Brutsaert stability functions there may be multiple solutions for ${\widehat{U}}_{*}$.

Here a slight modification of the usual Monin–Obukhov similarity theory is presented that ensures the existence of a unique solution for U_* even when Monin–Obukhov similarity theory is outside its range of validity. (The modification does not, of course, extend that range of validity.) We focus on the case in which F^b is known and the task is to determine U_*. The case in which F^b is also to be determined is more complicated (Bull and Derbyshire, 1990).

Since the expressions used for stability functions are typically derived by curve fitting to observations or simulations, it seems reasonable to attempt to provide an inverse to Eq. (A15) with aid of some curve fitting. The fitting can be done in such a way as to guarantee the existence of a unique solution for U_* satisfying the reasonable requirement $U_{*}\to \mathrm{0}$ as U(z)→0.

It is convenient to define $r=z_{\mathrm{0}}/z={\mathit{\zeta }}_{\mathrm{0}}/\mathit{\zeta }$ and to introduce a non-dimensional inverse flow speed,

$$\begin{array}{}\text{(A16)}& \widehat{s}={\displaystyle \frac{{\left(-F^{\mathrm{b}}z\right)}^{\mathrm{1}/\mathrm{3}}}{\mathit{\kappa }\left(U\left(z_{\mathrm{1}}\right)+\mathit{\epsilon }\right)}},\end{array}$$

where ε is a small safety parameter to prevent infinite $\widehat{s}$ as U(z₁)→0. Then Eq. (A15) becomes

$$\begin{array}{}\text{(A17)}& \widehat{s}={\displaystyle \frac{{\mathit{\zeta }}^{\mathrm{1}/\mathrm{3}}}{\mathit{\{}-\log r+\mathrm{\Psi }(\mathit{\zeta },r\mathit{\zeta }\left)\mathit{\right\}}}}.\end{array}$$

To parameterize the surface momentum flux we need to be able to compute $\mathit{\zeta }=\mathit{\zeta }(\widehat{s},r)$; ζ is a function of the two parameters $\widehat{s}$ and r.

To construct a functional fit for $\mathit{\zeta }(\widehat{s},r)$, proceed as follows. For both the stable case and the unstable case, carry out an asymptotic expansion of Eq. (A17) for the limits of small $\left|\widehat{s}\right|$ and large $\left|\widehat{s}\right|$ and rearrange to obtain the limiting expressions for $\mathit{\zeta }(\widehat{s},r)$. Then seek a simple expression that agrees with the asymptotic expressions in the two limits.

For the unstable case we begin with the integrated stability function given by Benoit (1977):

$$\begin{array}{}\text{(A18)}& \begin{array}{rl}\mathrm{\Psi }(\mathit{\zeta },r\mathit{\zeta })& =\log \left\{{\displaystyle \frac{\left(x_{\mathrm{0}}^{\mathrm{2}}+\mathrm{1}\right){\left(x_{\mathrm{0}}+\mathrm{1}\right)}^{\mathrm{2}}}{\left(x^{\mathrm{2}}+\mathrm{1}\right)(x+\mathrm{1}{)}^{\mathrm{2}}}}\right\}\\ & +\mathrm{2}\left(\arctan \left(x\right)-\arctan \left(x_{\mathrm{0}}\right)\right),\end{array}\end{array}$$

where

$$\begin{array}{}\text{(A19)}& x=(\mathrm{1}-\mathrm{15}\mathit{\zeta }{)}^{\mathrm{1}/\mathrm{4}},x_{\mathrm{0}}=(\mathrm{1}-\mathrm{15}r\mathit{\zeta }{)}^{\mathrm{1}/\mathrm{4}}.\end{array}$$

We find

$$\begin{array}{}\text{(A20)}& {\mathit{\zeta }}^{\mathrm{1}/\mathrm{3}}\sim \widehat{s}\left\{-\log r-{\displaystyle \frac{\mathrm{15}}{\mathrm{4}}}(\mathrm{1}-r)(\log r{)}^{\mathrm{3}}{\widehat{s}}^{\mathrm{3}}\right\}\mathrm{as}\widehat{s}\to \mathrm{0}\end{array}$$

and

$$\begin{array}{}\text{(A21)}& {\mathit{\zeta }}^{\mathrm{1}/\mathrm{3}}\sim \widehat{s}{\left\{-{\displaystyle \frac{\mathrm{1}}{\mathrm{15}}}{\left(\mathrm{4}\left(r^{-\mathrm{1}/\mathrm{4}}-\mathrm{1}\right)\right)}^{\mathrm{4}}{\widehat{s}}^{-\mathrm{3}}\right\}}^{\mathrm{1}/\mathrm{7}}\mathrm{as}\widehat{s}\to \mathrm{\infty }.\end{array}$$

These asymptotic limits are captured by the functional fit

$$\begin{array}{}\text{(A22)}& {\mathit{\zeta }}^{\mathrm{1}/\mathrm{3}}={\left(\stackrel{\mathrm{̃}}{a}+\stackrel{\mathrm{̃}}{b}{\left(\stackrel{\mathrm{̃}}{c}+(-\widehat{s}{)}^{\stackrel{\mathrm{̃}}{p}}\right)}^{\stackrel{\mathrm{̃}}{q}}\right)}^{\stackrel{\mathrm{̃}}{r}},\end{array}$$

with the parameters given by

$$\begin{array}{}\text{(A23)}& {\displaystyle }\stackrel{\mathrm{̃}}{p}=\mathrm{3};\stackrel{\mathrm{̃}}{q}=\mathrm{2}/\mathrm{7};\stackrel{\mathrm{̃}}{r}=\mathrm{3}/\left(\mathrm{7}\stackrel{\mathrm{̃}}{p}\stackrel{\mathrm{̃}}{q}\right);\text{(A24)}& {\displaystyle }\begin{array}{rl}\stackrel{\mathrm{̃}}{A}=-\log r;& \stackrel{\mathrm{̃}}{B}={\displaystyle \frac{\mathrm{15}}{\mathrm{4}}}(\mathrm{1}-r)(\log r{)}^{\mathrm{3}};\\ & \stackrel{\mathrm{̃}}{C}={\left\{{\displaystyle \frac{\mathrm{1}}{\mathrm{15}}}{\left(\mathrm{4}\left(r^{-\mathrm{1}/\mathrm{4}}-\mathrm{1}\right)\right)}^{\mathrm{4}}\right\}}^{\mathrm{1}/\mathrm{7}};\end{array}\text{(A25)}& {\displaystyle }\stackrel{\mathrm{̃}}{a}={\stackrel{\mathrm{̃}}{A}}^{\mathrm{1}/\stackrel{\mathrm{̃}}{r}}-\stackrel{\mathrm{̃}}{b}{\stackrel{\mathrm{̃}}{c}}^{\stackrel{\mathrm{̃}}{q}};\stackrel{\mathrm{̃}}{b}={\stackrel{\mathrm{̃}}{C}}^{\mathrm{1}/\stackrel{\mathrm{̃}}{r}};\stackrel{\mathrm{̃}}{c}={\left({\displaystyle \frac{\stackrel{\mathrm{̃}}{B}{\stackrel{\mathrm{̃}}{A}}^{(\mathrm{1}-\stackrel{\mathrm{̃}}{r})/\stackrel{\mathrm{̃}}{r}}}{\stackrel{\mathrm{̃}}{b}\stackrel{\mathrm{̃}}{q}\stackrel{\mathrm{̃}}{r}}}\right)}^{\mathrm{1}/(\stackrel{\mathrm{̃}}{q}-\mathrm{1})}.\end{array}$$

For the stable case we begin with a modification of the integrated stability function given by Cheng and Brutsaert (2005):

$$\begin{array}{}\text{(A26)}& \mathrm{\Psi }(\mathit{\zeta },r\mathit{\zeta })=a\log \left({\displaystyle \frac{\mathit{\zeta }+{\left(c^{\mathrm{b}}+{\mathit{\zeta }}^{\mathrm{b}}\right)}^{\mathrm{1}/\mathrm{b}}}{{\mathit{\zeta }}_{\mathrm{0}}+{\left(c^{\mathrm{b}}+{\mathit{\zeta }}_{\mathrm{0}}^{\mathrm{b}}\right)}^{\mathrm{1}/\mathrm{b}}}}\right),\end{array}$$

with a=3, b=2.5, and c=0.5. (The original Cheng and Brutsaert (2005) scheme corresponds to a=6.1, b=2.5, and c=1. Figure A1 shows the effect of this modification.) We find

$$\begin{array}{}\text{(A27)}& {\mathit{\zeta }}^{\mathrm{1}/\mathrm{3}}\sim \widehat{s}\left\{-\log r+{\displaystyle \frac{a}{c}}(\mathrm{1}-r)(-\log r{)}^{\mathrm{3}}{\widehat{s}}^{\mathrm{3}}\right\}\mathrm{as}\widehat{s}\to \mathrm{0}\end{array}$$

and

$$\begin{array}{}\text{(A28)}& {\mathit{\zeta }}^{\mathrm{1}/\mathrm{3}}\sim \widehat{s}\mathit{\{}-(\mathrm{1}+a\left)\log r\mathit{\right\}}\mathrm{as}\widehat{s}\to \mathrm{\infty }.\end{array}$$

These asymptotic limits are captured by the functional fit

$$\begin{array}{}\text{(A29)}& {\mathit{\zeta }}^{\mathrm{1}/\mathrm{3}}=\stackrel{\mathrm{̃}}{c}+{\left(\stackrel{\mathrm{̃}}{a}+\stackrel{\mathrm{̃}}{b}\widehat{s}{)}^{\stackrel{\mathrm{̃}}{p}}\right)}^{\stackrel{\mathrm{̃}}{q}},\end{array}$$

with the parameters given by

$$\begin{array}{}\text{(A30)}& {\displaystyle }\stackrel{\mathrm{̃}}{p}=\mathrm{3};\stackrel{\mathrm{̃}}{q}=-\mathrm{3};\text{(A31)}& {\displaystyle }\begin{array}{rl}\stackrel{\mathrm{̃}}{A}=-\log r;& \stackrel{\mathrm{̃}}{B}={\displaystyle \frac{a}{c}}(\mathrm{1}-r)(-\log r{)}^{\mathrm{3}};\\ & \stackrel{\mathrm{̃}}{C}=-(\mathrm{1}+a)\log r;\end{array}\text{(A32)}& {\displaystyle }\stackrel{\mathrm{̃}}{a}=(\stackrel{\mathrm{̃}}{A}-\stackrel{\mathrm{̃}}{c}{)}^{\mathrm{1}/\stackrel{\mathrm{̃}}{q}};\stackrel{\mathrm{̃}}{b}=\stackrel{\mathrm{̃}}{B}{\stackrel{\mathrm{̃}}{a}}^{\mathrm{1}/\stackrel{\mathrm{̃}}{q}}/\stackrel{\mathrm{̃}}{q};\stackrel{\mathrm{̃}}{c}=\stackrel{\mathrm{̃}}{C}.\end{array}$$