The Hydro-ABC model (Version 2.0): a simplified convective-scale model with moist dynamics

The prediction of convection (in terms of position, timing, and strength) is important to achieve for high-resolution weather forecasting. This problem requires not only good convective-scale models, but also data assimilation systems that give initial conditions which neither improperly hinder nor improperly hasten convection in the ensuing forecasts. Solving this problem is difficult and expensive using operational-scale numerical weather prediction systems, and so a simplified model of convective-scale flow is under development (called the “ABC model”). This paper extends the existing ABC model of dry convective-scale flow to include mixing ratios of vapour and condensate phases of water. The revised model is called “Hydro-ABC”.

Hydro-ABC includes transport of the vapour and condensate mixing ratios within a dynamical core, and it transitions between these two phases via a micro-physics scheme. A saturated mixing ratio is derived from model quantities, which helps determine whether evaporation or condensation happens. Latent heat is exchanged with the buoyancy variable (ABC's potential-temperature-like variable) in such a way to conserve total energy, where total energy is the sum of dry energy and latent heat. The model equations are designed to conserve the domain-total mass, water, and energy.

An example numerical model integration is performed and analysed, which shows the development of a realistic looking anvil cloud and excitation of inertio-gravity and acoustic modes over a wide range of frequencies. This behaviour means that Hydro-ABC is a sufficiently challenging model which will allow experimentation with innovative data assimilation strategies in future work. An ensemble of Hydro-ABC integrations is performed in order to study the possible forecast error covariance statistics (knowledge of which is necessary for data assimilation). These show patterns that are dependent on the presence of convective activity (at any model's vertical column), thus giving a taste of flow-dependent error statistics. Candidate indicators/harbingers of convection are also evaluated (namely relative humidity, hydrostatic imbalance, horizontal divergence, convective available potential energy, convective inhibition, vertical wind, and the condensate mixing ratio), some of which appear to be reliable diagnostics concerning the presence of convection. These diagnostics will be useful in the selection of the relevant forecast error covariance statistics when data assimilation for Hydro-ABC is developed.

1 Introduction

Numerical models used for modern numerical weather prediction (NWP) services have evolved rapidly in sophistication and scale over recent decades. One of the most notable advances concerns the reduction in size of the models' grid boxes. Using smaller grid sizes not only allows for finer flows to be resolved, but also permits new flow regimes that could not be represented with coarser models (Clark et al., 2016; Yano et al., 2018). These flow regimes are more non-linear and inherently less predictable than larger-scale flows (Hohenegger and Schär, 2007; Leung et al., 2019; Lorenz, 1969), and thus there is an increased need for more observations to constantly correct for fast-growing errors in order to produce useful forecasts (e.g. Banos et al., 2022). This is not just a dynamical effect (for instance due to larger Rossby or Reynolds numbers), but also due to the role of phase transitions of water (namely evaporation and condensation, thus absorbing or releasing latent heat; Errico et al., 2007), which can vary over short length scales, especially where moist processes occur (Montmerle and Berre, 2010).

By resolving moist processes that lead to convective motion (partially achieved with kilometre-scale grid lengths), some aspects of numerical models have become simpler, as much of the fine-scale transport of water and latent heat can be handled explicitly, thus removing the need for a convection scheme. Conversely, the necessary techniques needed to realistically estimate the initial conditions (namely data assimilation, DA) have become more complicated owing to the breakdown in many of the assumptions that are made in traditional DA methods, like linearity of the models and observation operators, Gaussianity of the background and observation errors, and the homogeneity and quasi-static nature of the background error statistics (Bannister et al., 2020). These (and other) problems contribute to the difficulty of assimilating data (especially radar data) in such models (Fabry and Meunier, 2020).

High-resolution models of the atmosphere and their DA counterparts remain very expensive to run, requiring costly supercomputers. It is therefore desirable to study the convective-scale DA problem in simplified (or “toy”) models that are capable of being run at less cost on smaller computers. One of the first such simplified models was presented by Würsch and Craig (2014). Their model is based on the non-rotating 1D shallow-water equations, comprising an equation each for fluid velocity, u; fluid depth, h; and rainwater, r. A further variable, the geopotential, ϕ, is tied to h in a non-trivial way. When the total height, H+h (where H is the topography height), is less than a specified level, the system is considered sub-saturated and ϕ and h are related via $\mathit{\varphi }=g(H+h)$ (where g is the acceleration due to gravity). When h is above the specified level, a cloud is assumed to form and an analogue of a “low-pressure” region is formed by redefining $\mathit{\varphi }=gH+{\mathit{\varphi }}_{\mathrm{c}}$ (where ϕ_c has a small value). This step encourages mass convergence and a form of convection to be initiated. Rain is then formed when h is above a second specified level in regions of mass convergence. The model exhibits convection-like behaviour, showing the sporadic appearance of clouds. Such a model and its developments (Kent et al., 2017) provide useful first steps to providing a cheap and effective model on which to base developments in convective-scale data assimilation.

The ABC model (Petrie et al., 2017) was separately developed as a low-cost model of convective-scale behaviour also to be used as a basis to investigate convective DA methods. The ABC model is based on the three-dimensional Euler equations but is reduced to two dimensions (longitude and height). It is therefore more complex than the shallow-water-based model of Würsch and Craig (2014), most notably in that it allows a height dependence. A variational DA capability was later developed for the ABC model (Bannister, 2020). This has since been used to show that enforcing geostrophic and hydrostatic balances as part of the background error covariance matrix (a key ingredient of the traditional variational DA problem) affects the DA's ability to analyse observations (Bannister, 2021). That study found that these balances, in the mid-latitude cases, are an advantage when analysing the larger scales but are a disadvantage at small scales. The model has also been used in a tropical setting, where the DA method has been expanded to include a hybrid (ensemble-variational) capability (Lee et al., 2022). Until now however, the ABC model has lacked a water variable, meaning that studies have been limited to dry dynamics only, where the model shows little evidence of sporadic convective behaviour. In this paper we show how moist processes can be included in the ABC model and how such behaviour emerges, providing a more appropriate setting for DA schemes to be tested and developed.

2 Dry ABC to Hydro-ABC

We seek modifications to the original ABC equations (Petrie et al., 2017), which are as simple as possible and require minimal change in the underlying numerical advection scheme. We refer to this moist model as “Hydro-ABC”. The ABC model comprises five variables on a two-dimensional (longitude–height) grid: zonal wind u, meridional wind v, vertical wind w, scaled density perturbation ${\stackrel{\mathrm{̃}}{\mathit{\rho }}}^{\prime }$ (a pressure-like variable), and buoyancy perturbation b^′ (a potential-temperature-like variable). Hydro-ABC contains two new variables, namely water vapour q and condensate q_c mixing ratios, the latter of which may also be thought of as a cloud variable. The model allows exchange between these two species according to the local thermodynamic conditions, with an associated latent heat exchange (see below). In order to keep the model as simple as possible, no rain is currently represented. The continuous dry equations support mass and total energy conservation, as shown in Petrie et al. (2017). Hydro-ABC must maintain similar conservation properties but with the addition of conservation of total water, q+q_c, and a total energy that includes a contribution from latent heat.

2.1 The Hydro-ABC model equations

The continuous Hydro-ABC equations are as follows:

$$\begin{array}{}\text{(1a)}& {\displaystyle \frac{\partial u}{\partial t}}+B\mathit{u}\cdot \mathrm{\nabla }u+C{\displaystyle \frac{\partial {\stackrel{\mathrm{̃}}{\mathit{\rho }}}^{\prime }}{\partial x}}-fv=\mathrm{0},\text{(1b)}& {\displaystyle \frac{\partial v}{\partial t}}+B\mathit{u}\cdot \mathrm{\nabla }v+fu=\mathrm{0},\text{(1c)}& {\displaystyle \frac{\partial w}{\partial t}}+B\mathit{u}\cdot \mathrm{\nabla }w+C{\displaystyle \frac{\partial {\stackrel{\mathrm{̃}}{\mathit{\rho }}}^{\prime }}{\partial z}}-b^{\prime }=\mathrm{0},\text{(1d)}& {\displaystyle \frac{\partial {\stackrel{\mathrm{̃}}{\mathit{\rho }}}^{\prime }}{\partial t}}+B\mathrm{\nabla }\cdot \left(\stackrel{\mathrm{̃}}{\mathit{\rho }}\mathit{u}\right)=\mathrm{0},\text{(1e)}& {\displaystyle \frac{\partial b^{\prime }}{\partial t}}+B\mathit{u}\cdot \mathrm{\nabla }{b}^{\prime }+A^{\mathrm{2}}w={\mathcal{S}}_{{\mathrm{b}}^{\prime }},\text{(1f)}& {\displaystyle }{p}^{\prime }=C{\mathit{\rho }}_{\mathrm{0}}{\stackrel{\mathrm{̃}}{\mathit{\rho }}}^{\prime },\text{(1g)}& {\displaystyle \frac{\partial \stackrel{\mathrm{̃}}{\mathit{\rho }}q}{\partial t}}+B\mathrm{\nabla }\cdot \left(\stackrel{\mathrm{̃}}{\mathit{\rho }}q\mathit{u}\right)=\stackrel{\mathrm{̃}}{\mathit{\rho }}(\text{Ev}-\text{Co}),\text{(1h)}& {\displaystyle \frac{\partial \stackrel{\mathrm{̃}}{\mathit{\rho }}{q}_{\mathrm{c}}}{\partial t}}+B\mathrm{\nabla }\cdot \left(\stackrel{\mathrm{̃}}{\mathit{\rho }}{q}_{\mathrm{c}}\mathit{u}\right)=-\stackrel{\mathrm{̃}}{\mathit{\rho }}(\text{Ev}-\text{Co}).\end{array}$$

Here A is the pure gravity wave frequency, B is the modulation of the advection and mass divergence terms, C is the inverse compressibility coefficient, t is time, x and z are positions in the longitude and height directions respectively, $\mathit{u}=(u,v,w)$, f is the Coriolis parameter, ρ₀ is a reference density, $\stackrel{\mathrm{̃}}{\mathit{\rho }}$ is the scaled density ($\stackrel{\mathrm{̃}}{\mathit{\rho }}=({\mathit{\rho }}_{\mathrm{0}}+{\mathit{\rho }}^{\prime })/{\mathit{\rho }}_{\mathrm{0}}$), p^′ is the pressure perturbation (diagnosed purely from ${\stackrel{\mathrm{̃}}{\mathit{\rho }}}^{\prime }$), ${\mathcal{S}}_{{\mathrm{b}}^{\prime }}$ is the buoyancy source term (related to latent heat exchange, which is activated when there is a change in water phase), Ev is the evaporation rate, Co is the condensation rate, and $\mathrm{\nabla }=\left(\begin{array}{cc}\partial /\partial x& \partial /\partial z\end{array}\right)$ is the two-dimensional vector derivative. The key differences between Hydro-ABC and the dry ABC system are ${\mathcal{S}}_{{\mathrm{b}}^{\prime }}$ in Eq. (1e) and the additional equations, Eqs. (1g) and (1h). How ${\mathcal{S}}_{{\mathrm{b}}^{\prime }}$ and Ev−Co are parametrised in Hydro-ABC as described below.

2.2 Conservation properties of the continuous equations

Petrie et al. (2017) (their Appendix B) use Gauss' divergence theorem, the model's vertical boundary conditions (their Table 2), and periodic lateral boundary conditions to show that the dry equations conserve total mass, $\int \int \mathrm{d}x\mathrm{d}z\mathit{\rho }$, and total dry energy, $\int \int \mathrm{d}x\mathrm{d}z{E}_{\mathrm{dry}}$, where the integration is performed over the whole domain (E_dry is defined in Sect. 2.2.2). Here we show that Hydro-ABC's Eqs. (1a–h) also conserve total water and total moist energy.

2.2.1 Conservation of total water

Total water at a particular location, q_t, is the sum of the vapour and condensate mixing ratios, $q_{\mathrm{t}}=q+q_{\mathrm{c}}$. Adding Eqs. (1g) and (1h) gives the governing equation for q_t: $\partial \left(\stackrel{\mathrm{̃}}{\mathit{\rho }}{q}_{\mathrm{t}}\right)/\partial t+B\mathrm{\nabla }\cdot \left(\stackrel{\mathrm{̃}}{\mathit{\rho }}{q}_{\mathrm{t}}\mathit{u}\right)=\mathrm{0}$. Multiplying this equation by the constant reference density ρ₀, recognising that $\mathit{\rho }={\mathit{\rho }}_{\mathrm{0}}\stackrel{\mathrm{̃}}{\mathit{\rho }}$, and integrating over the domain give the following equation for the total amount of water in the domain:

$$\begin{array}{}\text{(2)}& {\displaystyle \frac{\partial }{\partial t}}\int \int \mathit{\rho }{q}_{\mathrm{t}}\mathrm{d}x\mathrm{d}z+B\int \int \mathrm{\nabla }\cdot \left(\mathit{\rho }{q}_{\mathrm{t}}\mathit{u}\right)\mathrm{d}x\mathrm{d}z=\mathrm{0}.\end{array}$$

By following the same procedure as that in Petrie et al. (2017) (their Appendix B), it is straightforward to show that the second term is zero, thus demonstrating that the total amount of water in the domain is constant in time.

2.2.2 Conservation of total moist energy

From Petrie et al. (2017) the kinetic energy is $E_{\mathrm{k}}=\stackrel{\mathrm{̃}}{\mathit{\rho }}(u^{\mathrm{2}}+v^{\mathrm{2}}+w^{\mathrm{2}})/\mathrm{2}$, the buoyancy energy is $E_{\mathrm{b}}=\stackrel{\mathrm{̃}}{\mathit{\rho }}{b^{\prime }}^{\mathrm{2}}/\left(\mathrm{2}{A}^{\mathrm{2}}\right)$, and the elastic energy is $E_{\mathrm{e}}=C{{\stackrel{\mathrm{̃}}{\mathit{\rho }}}^{\prime }}^{\mathrm{2}}/\left(\mathrm{2}B\right)$. These make up the total dry energy, $E_{\mathrm{dry}}=E_{\mathrm{k}}+E_{\mathrm{b}}+E_{\mathrm{e}}$. In Hydro-ABC there is an additional latent heat energy which we define as

$$\begin{array}{}\text{(3)}& E_{\mathrm{l}}=\stackrel{\mathrm{̃}}{\mathit{\rho }}{L}_{\mathrm{v}}q,\end{array}$$

where L_v is the specific latent heat of vaporisation of water. In Hydro-ABC, the total moist energy is defined as the sum of the dry energy and latent heat, $E_{\mathrm{moist}}=E_{\mathrm{dry}}+E_{\mathrm{l}}=E_{\mathrm{k}}+E_{\mathrm{b}}+E_{\mathrm{e}}+E_{\mathrm{l}}$.

In its continuous form, Hydro-ABC is designed to conserve $\int \int \mathrm{d}x\mathrm{d}z{E}_{\mathrm{moist}}$ (see Appendix A). We, however, solve Hydro-ABC's Eq. (1) over a time step Δt in two parts: a “dynamical core” in which ${\mathcal{S}}_{{\mathrm{b}}^{\prime }}$, Ev, and Co are each zero, followed by a micro-physics step. A similar time-step splitting is done in other models, e.g. WRF (Weather Research and Forecasting model; Skamarock et al., 2021). The dynamical core is equivalent to integrating the dry equations (with extra advection equations for q and q_c). The evolution of the total energy is represented by the following sum:

$$\begin{array}{}\text{(4)}& \begin{array}{rl}& {\displaystyle \frac{\partial E_{\mathrm{moist}}}{\partial t}}=\\ & \underset{{\left(\partial E_{\mathrm{moist}}/\partial t\right)}^{\mathrm{dc}}}{\underbrace{{\left({\displaystyle \frac{\partial E_{\mathrm{dry}}}{\partial t}}\right)}^{\mathrm{dc}}+{\left({\displaystyle \frac{\partial E_{\mathrm{l}}}{\partial t}}\right)}^{\mathrm{dc}}}}+\underset{{\left(\partial E_{\mathrm{moist}}/\partial t\right)}^{\mathrm{mp}}}{\underbrace{{\left({\displaystyle \frac{\partial E_{\mathrm{l}}}{\partial t}}\right)}^{\mathrm{mp}}+{\left({\displaystyle \frac{\partial E_{\mathrm{b}}}{\partial t}}\right)}^{\mathrm{mp}}}},\end{array}\end{array}$$

where the large brackets indicate that the changes are performed over the dynamical core (dc) or the micro-physics (mp).

The first term represents pointwise changes in the dry energy by the dynamical core; this is the part that is shown in Petrie et al. (2017) to integrate to zero over the model's domain. The second term represents pointwise changes in the latent heat energy by the dynamical core, which is driven by transport of q. Although q (and hence E_l) can change at a particular point, the dynamical core does not change the phase of water, so this term also integrates to zero over the model's domain. The third and fourth terms are each non-zero, but the micro-physics scheme is designed so that their pointwise sum is zero (see Sect. 2.3). Hence integrating $\partial E_{\mathrm{moist}}/\partial t$ over the entire domain is zero. Note that the reason why it is the buoyancy energy that appears in the last term in the micro-physics step in Eq. (4) (and not kinetic or elastic energies) is explained below. In Sect. 3.2 we examine how well such energy conservation is approximated in the time–space discretised equations.

2.3 The parametrisations of ${\mathcal{S}}_{{\mathrm{b}}^{\prime }}$, Ev, and Co

As explained above, Hydro-ABC solves Eq. (1) over a time step Δt in two stages, by first assuming ${\mathcal{S}}_{{\mathrm{b}}^{\prime }}$ and Ev−Co are each zero (the dynamical core) and then adding their effects in a separate step (micro-physics, described below). The pointwise conservation of energy for the micro-physics step, Δ^mpE_moist=0 (where ${\mathrm{\Delta }}^{\mathrm{mp}}•\approx {\left(\partial •/\partial t\right)}^{\mathrm{mp}}\mathrm{\Delta }t$), is the basis for our micro-physics parametrisation, which is now described. The purpose of the micro-physics step is to determine values of ${\mathcal{S}}_{{\mathrm{b}}^{\prime }}$ and Ev−Co. There are no general closed-form expressions for these quantities, so some equations need to be solved implicitly, which are established in this section.

The relationship between the micro-physical change in q and Ev−Co is as follows: ${\mathrm{\Delta }}^{\mathrm{mp}}q=\left(\text{Ev}-\text{Co}\right)\mathrm{\Delta }t$. Further, the relationship between the change in E_l and the change in q is ${\mathrm{\Delta }}^{\mathrm{mp}}{E}_{\mathrm{l}}=\stackrel{\mathrm{̃}}{\mathit{\rho }}{L}_{\mathrm{v}}{\mathrm{\Delta }}^{\mathrm{mp}}q$. The above constraint that Δ^mpE_moist=0 may be put in the following way: ${\mathrm{\Delta }}^{\mathrm{mp}}{E}_{\mathrm{moist}}={\mathrm{\Delta }}^{\mathrm{mp}}{E}_{\mathrm{k}}+{\mathrm{\Delta }}^{\mathrm{mp}}{E}_{\mathrm{b}}+{\mathrm{\Delta }}^{\mathrm{mp}}{E}_{\mathrm{e}}+{\mathrm{\Delta }}^{\mathrm{mp}}{E}_{\mathrm{l}}=\mathrm{0}$. To facilitate closure, we assume that buoyancy energy is the only form of the dry energy that is exchanged with the latent heat. This is a reasonable assumption because in ABC, buoyancy is a potential-temperature-like variable (hence only E_b appears in the micro-physics terms of Eq. 4). The last equation then simplifies to ${\mathrm{\Delta }}^{\mathrm{mp}}{E}_{\mathrm{b}}+{\mathrm{\Delta }}^{\mathrm{mp}}{E}_{\mathrm{l}}=\mathrm{0}$. Combining the above formulae with the form of E_b (which is given in the opening paragraph of Sect. 2.2.2) yields

$$\begin{array}{}\text{(5)}& {\displaystyle \frac{\stackrel{\mathrm{̃}}{\mathit{\rho }}}{\mathrm{2}{A}^{\mathrm{2}}}}{\mathrm{\Delta }}^{\mathrm{mp}}\left({b^{\prime }}^{\mathrm{2}}\right)=-\stackrel{\mathrm{̃}}{\mathit{\rho }}{L}_{\mathrm{v}}\left(\text{Ev}-\text{Co}\right)\mathrm{\Delta }t\text{or}\text{(6)}& {\displaystyle }{\mathrm{\Delta }}^{\mathrm{mp}}\left({b^{\prime }}^{\mathrm{2}}\right)=-\mathrm{2}{A}^{\mathrm{2}}{L}_{\mathrm{v}}\left(\text{Ev}-\text{Co}\right)\mathrm{\Delta }t.\end{array}$$

Note that ${\mathrm{\Delta }}^{\mathrm{mp}}{b}^{\prime }={\mathcal{S}}_{{\mathrm{b}}^{\prime }}\mathrm{\Delta }t$. Equation (6) will now be used with further rules (below) to determine how our micro-physics scheme will modify q, q_c, and b^′ (which will reveal how ${\mathcal{S}}_{{\mathrm{b}}^{\prime }}$ and Ev−Co are found). The description of the scheme may be read in combination with the pseudo-code presented in Appendix B.

Figure 1 Graphical representation of Eq. (7), showing the change in q, q_c, and b^′ over a micro-physics step. Water vapour mixing ratio is shown in red, condensate is in blue, and buoyancy is in green. The water vapour mixing ratio according to saturation is the dotted black line. The zero line for all quantities spans all panels. Case (1) concerns condensation of water vapour due to super-saturation and case (2) concerns evaporation of condensate into vapour where there are three possibilities. Case (3) is not shown, as it represents no micro-physical changes to the fields.

The task now is to determine (Ev−Co)Δt, which is the change in q (or minus the change in q_c) over the micro-physics step. The following descriptions are made with reference to steps in the pseudo-code in Appendix B (in the form of Ap-B, step x). Consider the following change in vapour over the micro-physics step, which depends on whether q is super- or sub-saturated after the dynamical core is run for one time step:

$$\begin{array}{ll}\text{(7)}& {\displaystyle }& {\displaystyle }\left(\text{Ev}-\text{Co}\right)\mathrm{\Delta }t={\displaystyle }& {\displaystyle }\left\{\begin{array}{ll}-\left(q-q_{\mathrm{s}}\right)& \left\{\begin{array}{c}\text{condensation}\\ \text{(super-saturated):}\\ \text{when }q\ge q_{\mathrm{s}}\text{ and}\\ \left(b^{\prime }>\mathrm{0}\text{ or }\right.\\ \left.\left[b^{\prime }<\mathrm{0}\text{ and }\stackrel{\mathrm{̃}}{\mathit{\rho }}{L}_{\mathrm{v}}\left|q-q_{\mathrm{s}}\right|/E_{\mathrm{b}}>\mathrm{\Gamma }\right]\right),\end{array}\right.\\ min\left(\left(q_{\mathrm{s}}-q\right)\left(\mathrm{1}-e^{-\mathrm{\Delta }t/\mathit{\tau }}\right),\right.\\ \left.q_{\mathrm{c}},{b^{\prime }}^{\mathrm{2}}/\left(\mathrm{2}{A}^{\mathrm{2}}{L}_{\mathrm{v}}\right)\right)& \left\{\begin{array}{c}\text{evaporation}\\ \text{(sub-saturated):}\\ \text{when }q<q_{\mathrm{s}}\\ \text{ and }{b}^{\prime }>\mathrm{0},\end{array}\right.\\ \mathrm{0}& \text{otherwise},\end{array}\right.\end{array}$$

where q_s is the saturated mixing ratio (this is a function of the model state; see Sect. 2.4), Γ is a threshold for convection (see below), τ is the prescribed evaporation timescale, and the “min” function selects the smallest of its three arguments. This equation is represented in Fig. 1. The three cases of this equation are now discussed.

• 1.

  The first case in Eq. (7) concerns condensation of water vapour due to super-saturation (Fig. 1a and Ap-B, step 3). In this case q will instantaneously decrease by the amount q−q_s (red line) and q_c will increase by the same amount (blue line). According to Eq. (6), associated with this step is an increase in ${b^{\prime }}^{\mathrm{2}}$ by the amount ${\mathrm{\Delta }}^{\mathrm{mp}}\left({b^{\prime }}^{\mathrm{2}}\right)=\mathrm{2}{A}^{\mathrm{2}}{L}_{\mathrm{v}}(q-q_{\mathrm{s}})$. This is due to the release of latent heat, which is converted to buoyancy energy. Writing ${\mathrm{\Delta }}^{\mathrm{mp}}\left({b^{\prime }}^{\mathrm{2}}\right)$ as ${b^{\prime }}_{\mathrm{mp}}^{\mathrm{2}}-{b^{\prime }}^{\mathrm{2}}$ (where $b_{\mathrm{mp}}^{\prime }$ is the value of buoyancy after micro-physics and b^′ is the value before), the equation for ${\mathrm{\Delta }}^{\mathrm{mp}}\left({b^{\prime }}^{\mathrm{2}}\right)$ translates to

  $$\begin{array}{}\text{(8)}& b_{\mathrm{mp}}^{\prime }=+\sqrt{{b^{\prime }}^{\mathrm{2}}+\mathrm{2}{A}^{\mathrm{2}}{L}_{\mathrm{v}}(q-q_{\mathrm{s}})}\end{array}$$

  (green line). The positive square root is assumed because b^′ is associated with potential temperature, which would rise after latent heat is released.

  • a.

    Note that q_s is shown in Sect. 2.4 to be a function of ${\stackrel{\mathrm{̃}}{\mathit{\rho }}}^{\prime }$ and b^′. Changing b^′ as prescribed by Eq. (8) will therefore change q_s. Equation (8) therefore represents a non-linear equation. This is solved by finding the root of the function $f\left(q_{\mathrm{mp}}\right)=q_{\mathrm{mp}}-q_{\mathrm{s}}({\stackrel{\mathrm{̃}}{\mathit{\rho }}}^{\prime },b_{\mathrm{mp}}^{\prime }(q_{\mathrm{mp}}\left)\right)$ given that $b_{\mathrm{mp}}^{\prime }$ is a function of q_mp, i.e. $b_{\mathrm{mp}}^{\prime }=+\sqrt{{b^{\prime }}^{\mathrm{2}}+\mathrm{2}{A}^{\mathrm{2}}{L}_{\mathrm{v}}(q-q_{\mathrm{mp}})}$, and assuming that ${\stackrel{\mathrm{̃}}{\mathit{\rho }}}^{\prime }$ is constant over the micro-physics step. When $b^{\prime }\ge \mathrm{0}$, this is done using a Newton–Raphson procedure while assuming ${\stackrel{\mathrm{̃}}{\mathit{\rho }}}^{\prime }$ is constant. The value of (Ev−Co)Δt for the first branch of Eq. (7) is then $-\left(q-q_{\mathrm{s}}({\stackrel{\mathrm{̃}}{\mathit{\rho }}}^{\prime },b_{\mathrm{mp}}^{\prime }(q_{\mathrm{mp}}\left)\right)\right)$. See Ap-B, step 3a.

  • b.

    Step a, however, becomes potentially unphysical if b^′ is large and negative and q−q_s is small, as this would flip the sign of b^′ to be positive, resulting is an unrealistically large increment in b^′. For this reason, included in the first case is the additional condition that $\stackrel{\mathrm{̃}}{\mathit{\rho }}{L}_{\mathrm{v}}\left|q-q_{\mathrm{s}}\right|/E_{\mathrm{b}}>\mathrm{\Gamma }$, where Γ≫1. This is the condition (which must be met only when $b^{\prime }<\mathrm{0}$) where, for condensation to happen, the release of latent heat must be much larger than the amount of buoyancy energy. The value of Γ=10 is chosen for the results in this paper. Also, for simplicity, when $b^{\prime }<\mathrm{0}$, the amount of latent heat released is based on the value of q_s at the pre-micro-physics value of b^′ and not on the value of q_s at the increased value of b^′ due to the latent heat release. This means that the non-linear equation in (a) above is not used in this case. See Ap-B, step 3b. The cause of this issue seems to be in the square (of b^′) and square root in Eq. (8), which forces the sign of b^′ to become invisible, followed by the need to choose the sign of the solution.

• 2.

  The second case in Eq. (7) concerns evaporation of condensate into vapour (see Ap-B, step 3c), where q is assumed to relax exponentially upwards towards q_s with timescale τ (Fig. 1b, red line) and q_c relaxes downwards (blue line). According to Eq. (6), associated with this step is a change in ${b^{\prime }}^{\mathrm{2}}$ by the amount ${\mathrm{\Delta }}^{\mathrm{mp}}\left({b^{\prime }}^{\mathrm{2}}\right)=-\mathrm{2}{A}^{\mathrm{2}}{L}_{\mathrm{v}}\left(\text{Ev}-\text{Co}\right)\mathrm{\Delta }t$ (where $\text{Ev}-\text{Co}>\mathrm{0}$). This translates to a new value $b_{\mathrm{mp}}^{\prime }$ according to the following:

  $$\begin{array}{}\text{(9)}& b_{\mathrm{mp}}^{\prime }=+\sqrt{{b^{\prime }}^{\mathrm{2}}-\mathrm{2}{A}^{\mathrm{2}}{L}_{\mathrm{v}}\left(\text{Ev}-\text{Co}\right)\mathrm{\Delta }t},\end{array}$$

  where the positive square root is taken, which corresponds to a reduction in b^′ due to evaporation (green line). This process is restricted to when $b^{\prime }>\mathrm{0}$ because, in the case of $b^{\prime }<\mathrm{0}$, Eq. (9) dictates that b^′ will increase (and will instantly become positive). Even if one chooses to take the negative square root in cases when $b^{\prime }<\mathrm{0}$, Eq. (9) will yield an increased b^′ (i.e. $b_{\mathrm{mp}}^{\prime }-b^{\prime }>\mathrm{0}$), which is unphysical.
  This process assumes that there is an adequate amount of condensate to evaporate and an adequate amount of buoyancy energy to convert to latent heat during evaporation (this corresponds to the first argument of the min function in Eq. 7 and to Ap-B, step 4a). The evaporation process is limited should either of these requirements not be met, which corresponds to the second and third arguments of the min function. When the second argument is smallest, evaporation is limited by available q_c (panel c, where the micro-physics step evaporates all available q_c, and see Ap-B, step 4b), and when the third argument is smallest, evaporation is limited by available buoyancy energy (panel d, where the micro-physics step converts all available E_b, and see Ap-B, step 4c). As evaporation is much slower than condensation, (Ev−Co)Δt (case 2 in Eq. 7) is based on the value of q_s before evaporation lowers the value of b^′ (and hence lowers q_s). This avoids the need to solve the non-linear equations mentioned previously (but overall non-linearity is not avoided due to the presence of the “switches” in Eq. 7).

• 3.

  The third case is invoked for the residual cases (i.e. for super-saturation but when the threshold for condensation is not met or for sub-saturation but when $b^{\prime }<\mathrm{0}$). In this case there are no micro-physics-related changes to the fields.

The change in buoyancy due to micro-physics, $b_{\mathrm{mp}}^{\prime }-b^{\prime }$, yields the value ${\mathcal{S}}_{{\mathrm{b}}^{\prime }}\mathrm{\Delta }t$ in Eq. (1e), and the change in vapour, q_mp−q, gives the value (Ev−Co)Δt in Eqs. (1g) and (1h).

2.4 The saturated vapour mixing ratio

The saturated mixing ratio, q_s, specifies the upper limit for water vapour. In standard atmospheric models it is a function of pressure and temperature. We use the following formula:

$$\begin{array}{}\text{(10)}& q_{\mathrm{s}}=\mathrm{1000}\mathrm{g}{\mathrm{kg}}^{-\mathrm{1}}{\displaystyle \frac{\mathrm{380}\mathrm{Pa}}{p}}\exp \left({\displaystyle \frac{\mathrm{17.3}(T-\mathrm{273.2}\mathrm{K})}{T-\mathrm{35.9}\mathrm{K}}}\right),\end{array}$$

which is based on Eq. (9.8) of Pielke (2002) (adapted so that pressure is in pascals and the mixing ratio is in grams per kilogram). In Eq. (10), p and T (and hence q_s) have an implicit dependence on x, z, and t. Neither p nor T is an ABC variable, and so we now take steps to produce physically reasonable formulae for p and T from ABC variables.

For pressure, a reference density profile, ρ₀(z), is required. In the ABC model (Petrie et al., 2017), ρ₀ is a constant, but for the purposes of computing q_s we let this depend on z according to ${\mathit{\rho }}_{\mathrm{0}}\left(z\right)={\mathit{\rho }}_{\mathrm{00}}\exp (-z/H)$, where ρ₀₀ is the surface air density (1.225 kg m⁻³) and H is the scale height (9 km). For pressure, $p=p_{\mathrm{0}}+p^{\prime }$, where p₀ (the reference pressure) is in hydrostatic balance with ρ₀, $\mathrm{d}{p}_{\mathrm{0}}\left(z\right)/\mathrm{d}z=-{\mathit{\rho }}_{\mathrm{0}}\left(z\right)g$, where g is the acceleration due to gravity. The solution to this equation, given ρ₀(z), is $p_{\mathrm{0}}\left(z\right)=p_{\mathrm{00}}-H{\mathit{\rho }}_{\mathrm{00}}g\left(\mathrm{1}-\exp (-z/H)\right)$, where p₀₀ is the reference surface pressure. The relationship between surface reference pressure and surface reference density is p₀₀=Hρ₀₀g, which is derived from the previous equation with the condition $p_{\mathrm{0}}(z\to \mathrm{\infty })\to \mathrm{0}$. ABC's equation of state, Eq. (1f), relates the ABC variable ${\stackrel{\mathrm{̃}}{\mathit{\rho }}}^{\prime }$ to p^′ and, using the above form of p_o(z), gives

$$\begin{array}{}\text{(11)}& \begin{array}{rl}p(x,z,t)& =p_{\mathrm{0}}\left(z\right)+p^{\prime }(x,z,t)=p_{\mathrm{00}}\\ & -H{\mathit{\rho }}_{\mathrm{00}}g\left(\mathrm{1}-\exp (-z/H)\right)\\ & +C{\mathit{\rho }}_{\mathrm{00}}\exp (-z/H){\stackrel{\mathrm{̃}}{\mathit{\rho }}}^{\prime }(x,z,t).\end{array}\end{array}$$

In Petrie et al. (2017), b^′ is related to a potential temperature increment θ^′ via ${\mathit{\theta }}^{\prime }=({\mathit{\theta }}_R/g)b^{\prime }$, where θ_R is a constant (value 273 K).¹ Hence we focus now on potential temperature. The potential temperature in Petrie et al. (2017) is defined as $\mathit{\theta }(x,z,t)={\mathit{\theta }}_R+{\mathit{\theta }}_{\mathrm{0}}\left(z\right)+{\mathit{\theta }}^{\prime }(x,z,t)$. To derive a reference profile, θ₀(z), we use the definition of the Brunt–Väisälä frequency $N^{\mathrm{2}}=(g/{\mathit{\theta }}_R)\mathrm{d}{\mathit{\theta }}_{\mathrm{0}}\left(z\right)/\mathrm{d}z$ (see e.g. Holton and Hakim, 2013). The parameter A in Eq. (1e) takes the place of N, and so we solve the equation $A^{\mathrm{2}}=(g/{\mathit{\theta }}_R)\mathrm{d}{\mathit{\theta }}_{\mathrm{0}}\left(z\right)/\mathrm{d}z$. The solution to this equation is ${\mathit{\theta }}_{\mathrm{0}}\left(z\right)={\mathit{\theta }}_{\mathrm{00}}-{\mathit{\theta }}_R+(A^{\mathrm{2}}{\mathit{\theta }}_R/g)z$, where θ₀₀ is the surface potential temperature (300 K). The potential temperature is therefore

$$\begin{array}{}\text{(12)}& \begin{array}{rl}\mathit{\theta }(x,z,t)& ={\mathit{\theta }}_R+{\mathit{\theta }}_{\mathrm{0}}\left(z\right)+{\mathit{\theta }}^{\prime }(x,z,t)\\ & ={\mathit{\theta }}_{\mathrm{00}}+{\displaystyle \frac{{\mathit{\theta }}_R}{g}}\left(A^{\mathrm{2}}z+b^{\prime }(x,z,t)\right).\end{array}\end{array}$$

The standard relationship for T in terms of θ and p is $T(x,z,t)=\mathit{\theta }(x,z,t){\left(p(x,z,t)/p_{\mathrm{00}}\right)}^{\mathit{\kappa }}$ (Holton and Hakim, 2013), where κ is the constant $R/c_{\mathrm{p}}$ (where R is the gas constant for dry air and c_p is the specific heat capacity at constant pressure; the value of κ is 0.286). Using Eqs. (12) and (11) gives

$$\begin{array}{}\text{(13)}& \begin{array}{rl}& T(x,z,t)=\left[{\mathit{\theta }}_{\mathrm{00}}+{\displaystyle \frac{{\mathit{\theta }}_R}{g}}\left(A^{\mathrm{2}}z+b^{\prime }(x,z,t)\right)\right]\\ & {\left[{\displaystyle \frac{p_{\mathrm{00}}-H{\mathit{\rho }}_{\mathrm{00}}g\left(\mathrm{1}-\exp (-z/H)\right)+C{\mathit{\rho }}_{\mathrm{00}}\exp (-z/H){\stackrel{\mathrm{̃}}{\mathit{\rho }}}^{\prime }(x,z,t)}{p_{\mathrm{00}}}}\right]}^{\mathit{\kappa }}.\end{array}\end{array}$$

Equations (11) and (13) can now be substituted into Eq. (10) to yield a physically reasonable function for the saturated vapour mixing ratio in terms of ABC variables b^′ and ${\stackrel{\mathrm{̃}}{\mathit{\rho }}}^{\prime }$.

2.5 Numerical details

The ABC model stores variables on an Arakawa C grid in the horizontal and a Charney–Phillips grid in the vertical (see Petrie et al., 2017, their Fig. 1). The two extra variables in Hydro-ABC, q and q_c, are each stored at the same location as $\stackrel{\mathrm{̃}}{\mathit{\rho }}$ and ${\stackrel{\mathrm{̃}}{\mathit{\rho }}}^{\prime }$. This choice makes sense from a numerical perspective as both q and q_c are multiplied by $\stackrel{\mathrm{̃}}{\mathit{\rho }}$ in Eqs. (1g) and (1h), and most of the finite differences that approximate the derivatives in these equations are evaluated at the correct locations on the grid.

In Hydro-ABC, Eqs. (1g) and (1h) are solved as advection/source equations rather than as the continuity equations shown. By employing the product rule for differentiation, Eq. (1g) becomes $\stackrel{\mathrm{̃}}{\mathit{\rho }}\partial q/\partial t+q\partial \stackrel{\mathrm{̃}}{\mathit{\rho }}/\partial t+Bq\mathrm{\nabla }\cdot \left(\stackrel{\mathrm{̃}}{\mathit{\rho }}\mathit{u}\right)+B\stackrel{\mathrm{̃}}{\mathit{\rho }}\mathit{u}\cdot \mathrm{\nabla }q=\stackrel{\mathrm{̃}}{\mathit{\rho }}(\text{Ev}-\text{Co})$. The second and third terms sum to zero due to the mass continuity (Eq. 1d), leaving the advection/source equation for q as follows:

$$\begin{array}{}\text{(14)}& {\displaystyle \frac{\partial q}{\partial t}}+B\mathit{u}\cdot \mathrm{\nabla }q=\text{Ev}-\text{Co}.\end{array}$$

A similar equation is also found for q_c.

The dynamical core's scheme is based on Cullen and Davies (1991) and has two parts: an adjustment stage and an advection stage. The adjustment stage operates over two sub-time steps and deals with the momentum and thermodynamic equations (omitting the advective terms) followed by treatment of the continuity equation. The advection stage advects the fields using the winds found in the adjustment stage averaged over the two sub-time steps. Upwind gradients are used for the advection. For more details see Sect. 3.3 of Petrie et al. (2017). The equations for q and q_c (with $\text{Ev}-\text{Co}=\mathrm{0}$) are each solved by means of the same advection scheme as for the other variables in ABC. Ev−Co is subsequently found with the micro-physics parametrisation of Sect. 2.3 using values of variables after the dynamical core is called.

3 Example run of the model

In this section we show example behaviour of Hydro-ABC. The initial conditions for the dry variables (u, v, w, ${\stackrel{\mathrm{̃}}{\mathit{\rho }}}^{\prime }$, and b^′) are set according to the prescription in Petrie et al. (2017) (their Sect. 5.2). This involves extracting u and v from a longitude–height slice output from a 1.5 km Unified Model (UM) run, which has the same grid staggering as the ABC model, and then adjusting u and v to have smooth periodic boundary conditions in x. Then ${\stackrel{\mathrm{̃}}{\mathit{\rho }}}^{\prime }$ is set using geostrophic balance, b^′ is set using hydrostatic balance, and w is set by imposing a 3D divergence-free condition. The initial water vapour is set as shown in the top panel of Fig. 2 (set so the maximum relative humidity is 0.98), and the initial q_c=0.

3.1 Reference run

The rest of Fig. 2 (rows 2, 3, and 4) shows the evolution of w (left column) and q, q_c (right column) over the first 6 h of the integration using the reference parameters indicated in the caption. Before the first hour of the simulation is complete there are no obvious signs of convection (not shown). It is only at t=60 min that localised upward motion starts to develop at ground level, most strongly at x=148.5 km but also at 156 km (also not shown). By 120 min, these have developed into striking convective plumes, which transport vapour upwards. This coincides with the position of the water vapour maximum in the initial conditions. In this region there is significant conversion of vapour to condensate, thus releasing latent heat, which is converted to buoyancy energy, thus inducing strong upward motion. The magnitude of the plume's vertical motion is very high (peaking at ∼200 m s⁻¹ at the strongest convecting point), but in the ABC model this speed is moderated by the small value of B (0.05), so it is equivalent to a peak vertical mass transport of 10 m s⁻¹, which is realistic. The vertical grid spacing is ∼250 m, making the vertical Courant number 0.08 for the unmodulated motion and 0.004 allowing for the B modulation. These values are well within the Courant–Friedrichs–Lewy condition for the Courant number to be less than unity.

Figure 2 Initial conditions of the example reference run of Hydro-ABC ($A=\mathrm{0.01}{\mathrm{s}}^{-\mathrm{1}}$, B=0.05, $C=\mathrm{2}\times {\mathrm{10}}^{\mathrm{4}}{\mathrm{m}}^{\mathrm{2}}{\mathrm{s}}^{-\mathrm{2}}$, $f={\mathrm{10}}^{-\mathrm{4}}{\mathrm{s}}^{-\mathrm{1}}$, Δt=0.1 s, τ=1000 s, $L_{\mathrm{v}}=\mathrm{2500}\mathrm{J}{\mathrm{g}}^{-\mathrm{1}}$, and Γ=10). The initial conditions of the dry variables are set as described in the text, the initial values of q are shown in the top panel (set so the maximum relative humidity is 98 %), and the initial q_c=0. Snapshots of the subsequent integration of Hydro-ABC are shown in the remaining panels, zoomed into the region indicated by the red box in the top panel for 120, 240, and 360 min. w is shown in the left column (m s⁻¹), and q (shading) and q_c (contours) are shown in the right column. The contours are the black (1 g kg⁻¹), orange (6 g kg⁻¹), and red (12 g kg⁻¹) lines. Three zones are defined, the “convection zone” (columns close to the convection), the “moist zone” (columns of non-zero moisture in the initial conditions), and the “dry zone” (zero moisture in the initial conditions).

By 240 min, an anvil cloud of vapour and condensate has established itself between heights of ∼5.5 and 8.5 km, which spreads out laterally in both directions. The zonal winds at these levels peak at around ±100 m s⁻¹ (not shown), but again once B moderation is taken into account; this is equivalent to a horizontal mass transport of ±5 m s⁻¹. At this time, the second smaller plume has been extinguished by downdraughts from the main plume, but smaller convective plumes of 11 to 13 km have developed on either side of the main plume. Between 240 and 360 min, the smaller plumes are themselves impeded by downwelling air on each side of the main plume. By 360 min the condensate anvils have continued to develop laterally, and the central plume has started the process of weakening.

3.2 Numerical effects on the total energy evolution

Hydro-ABC in its continuous form, as described by Eqs. (1a–h) and the algorithm in Sect. 2.3, will conserve total moist energy, $\int \int \mathrm{d}x\mathrm{d}z{E}_{\mathrm{moist}}\left(t\right)=\mathrm{constant}$. Discretisation of the equations in space and time will however lead to loss of energy conservation. Five numerical experiments – based on the reference parameters but of different Δt values from 0.1 to 1.0 s – are performed to analyse how in practice the total moist energy deviates from conservation. Figure 3 shows the evolution of relative total moist energy (relative to the total moist energy at t=0), where the labels are the Δt values of the model runs. All runs show fluctuations in this quantity, which become visible on the plot's scale when convection starts just after t=60 min. As expected Δt=0.1 s deviates from conservation the least, and Δt=1.0 s deviates the most. The slight increase in energy in the Δt=0.1 and 0.2 s integrations is small but unexpected. Our assumption is that the Δt=0.1 s setting represents an acceptable deviation from conservation, especially given the fast-moving nature of the features that are being simulated. Δt=0.1 s is the setting that has been used for the result in Fig. 2 and for the remainder of this paper.

Figure 3 Evolution of the relative total energy of runs of Hydro-ABC for different model integration time steps as indicated (in  s). The relative total energy is $\left[\int \int \mathrm{d}x\mathrm{d}z{E}_{\mathrm{moist}}\left(t\right)\right]/\left[\int \int \mathrm{d}x\mathrm{d}z{E}_{\mathrm{moist}}\left(\mathrm{0}\right)\right]$, where $E_{\mathrm{moist}}=E_{\mathrm{k}}+E_{\mathrm{b}}+E_{\mathrm{e}}+E_{\mathrm{l}}$ (Sect. 2.2.2).

3.3 Excitation of waves by convection

The ABC model is capable of supporting Rossby-like, inertio-gravity, and acoustic waves. The continuous equations can support all horizontal wavenumbers allowed by the periodic boundary conditions and all frequencies. In practice though, the maximum horizontal wavenumber and frequency are set by Δx (the grid box length) and Δt (the time step) respectively. In runs of the dry model, high wavenumbers and frequencies have a very small weight, but their presence is still important in order to represent scales that would otherwise anomalously dissipate energy. Allowing moist convection, as facilitated by Hydro-ABC, potentially excites appreciable weight to modes of high wavenumber and frequency.

Figure 4 shows w for level 30 (near the top of the anvil cloud) of two model runs: one for the dry model (left) and another for Hydro-ABC (right). The top panels are Hovmöller plots of w for the longitude–time domain. Both versions of the model show westward- and eastward-propagating waves. The bottom panels of Fig. 4 show the spectral density of w at the same height. These plots are simply the magnitude of the bi-Fourier transform of the top panels. This reveals the dispersion relations of the dominant waves in the system. By analogy with waves in systems that represent more closely the real atmosphere, each feature in these plots represents a wave with a different vertical wavenumber (Salby, 1996, their Fig. 14.9) and the lines separate into two groups. The first group is the lower-frequency inertio-gravity waves, whose frequencies increase with horizontal wavenumber and then level off at $A/\left(\mathrm{2}\mathit{\pi }\right)=\mathrm{1.59}\times {\mathrm{10}}^{-\mathrm{3}}{\mathrm{s}}^{-\mathrm{1}}$. This is the wave frequency of the highest-frequency inertio-gravity waves in the high wavenumber limit (the horizontal dashed lines), which is consistent with classical wave theory, e.g. Gill (1982) or Salby (1996). The second group is the higher-frequency acoustic waves, whose frequencies increase with horizontal wavenumber. The quantity $\sqrt{\mathrm{BC}}$ is the pure sound wave speed (the gradient of the sloping dashed lines). This lines separate the inertio-gravity waves and the acoustic waves and is occupied by Lamb waves in the real atmosphere (there is no visible evidence of Lamb waves in these (Hydro)-ABC runs, but we will still call this line the “Lamb line”). In the (Hydro)-ABC runs some of the high-wavenumber acoustic waves have frequencies lower than the Lamb line (visible especially in the top-right section of the bottom-right panel). This is a numerical effect where the frequencies of high-wavenumber waves are too low, leading to numerical dispersion (see e.g. Chap. 2 of Durran, 1999).

Comparing the dry ABC and Hydro-ABC spectra, the shapes of the dispersion curves look similar. This indicates that adding the adiabatic moisture scheme does not noticeably modify the underlying dispersion relations of the dynamical core. Instead, the striking effects that moist convection in Hydro-ABC causes are (i) the excitation of stronger inertio-gravity and acoustic modes across the spectrum and (ii) broadening of their frequencies. Effect (i) is visible at all wavenumbers and frequencies, especially for inertio-gravity waves for wavenumbers less than ∼70. Effect (ii) is likely due to the gradual onset of convection, where the degree of broadening is inversely proportional to the timescale of the onset of convection. In the top-right panel of Fig. 4, there is a gradual onset of the strength of convection, which leads to a convolution of the spectra with this “low-frequency” onset.

Figure 4 Vertical wind in the longitude–time (a, b) and wavenumber–frequency (c, d) domains for level 30 in the model (7.7 km height). The left panels are for dry ABC, and the right panels are for Hydro-ABC. The initial conditions of u, v, w, ${\stackrel{\mathrm{̃}}{\mathit{\rho }}}^{\prime }$, and b^′ are identical for ABC and Hydro-ABC, and the additional initial conditions for q and q_c are the same as those used in Fig. 2. In (c) and (d), the horizontal dashed lines are each at frequency $A/\left(\mathrm{2}\mathit{\pi }\right)$ (the pure gravity wave frequency), and the sloped dashed lines each have gradient $\mathrm{d}f/\mathrm{d}n=L^{-\mathrm{1}}\sqrt{\mathrm{BC}}$ (the pure sound wave speed), where f is frequency (the y axes), n is the wavenumber index (the x axes), and L is the length of the domain, 540×10³ m.

4 Convection-dependent covariances and identification of areas prone to convection

It is well known that latent heat release can lead to convective instabilities. Such rapidly changing conditions present challenges for DA, not least because the background error covariances can become highly flow-dependent (e.g. Lorenc, 2007; Montmerle and Berre, 2010; Sun et al., 2014; Bannister et al., 2020). In this section we show how vertical error statistics of model variables depend on the background conditions for hydrodynamically stable and unstable regions. An important issue in convective-scale DA is the identification of convective regions so that the (otherwise static) background error covariance matrix can be modified with the relevant vertical covariance statistics should the model be used with data assimilation.

4.1 How forecast statistics are affected by convectively prone regions

In order to give a taste of moisture-affected covariances and correlations between different variables and positions, as well as in different regimes in Hydro-ABC, we generate an ensemble of Hydro-ABC forecasts. Longitude–height slices from a 93-member ensemble of convective-scale UM forecasts (Bannister et al., 2017) are each processed in the same way as described in Sect. 3 to give a 93-member ensemble of Hydro-ABC forecasts (run out to 360 min). Each ensemble member is prepared with the same q and q_c fields (as in Fig. 2), but the relative humidity ($q/q_{\mathrm{s}}$, the field that matters for convection) is different between members because q_s is a function of b^′ and ${\stackrel{\mathrm{̃}}{\mathit{\rho }}}^{\prime }$ (see Eq. 10, Sect. 2.4, where p and T depend on ${\stackrel{\mathrm{̃}}{\mathit{\rho }}}^{\prime }$ and b^′ as shown in Eqs. 11 and 13). The example Hydro-ABC run in Sect. 3 is one of these members, and convection happens around the same place in all but one members, allowing the vertical statistics (covariances between different vertical levels) to be based on slightly different versions of the convective situation, which is centred at x=148.5 km as in Fig. 2.

The covariances of and between different model variables will have different units and magnitudes. In order to simplify comparison of covariances, we first normalise each model variable by its (constant) maximum standard deviation, where the maximum is calculated over all positions and times. This will result in covariances (referred to as normalised covariances) that are bounded between −1 and 1. These values are not correlations though as correlations are based on normalisation by the local standard deviations (local in space and time), although we do later also show correlations of different quantities.

Figure 5 Evolution of a selection of normalised vertical covariances in the convection zone (x=148.5 km) at times between 40 and 360 min. The rows are for w–w, u–w, q–w, and q_c–w covariances. In each plot of p₁–p₂ covariances, the y axis is for p₁ and the x axis is for p₂, and the surface is at the bottom left of each square. Red (blue) indicates positive (negative) covariances according to the logarithmic scale shown.

Figure 5 shows the normalised covariances for a selection of quantities (see figure caption) at different times in the convection zone (x=148.5 km). Near the start of the integrations (t=40 min), all covariances shown are too small to be detected on the scale used. At 80 min, convection is just starting, which appears in a region of non-zero covariances within 4.5 km of the surface. There are positive covariances between w and itself at different levels, and with q (and more weakly with q_c), and there are split negative–positive covariances with u. These results show that larger w is associated with greater vertical advection of vapour and the formation of a shear zone in zonal wind, with westward wind in the lower part and eastward wind in the upper part of this 4.5 km layer. By 120 min the non-zero covariance region has grown to within 9–10 km of the surface, including growth of the wind shear zone. An increase in the vertical length scales in q errors is also seen in the WRF-based study of Michel et al. (2011) and in the AROME-based study of Montmerle and Berre (2010), in each case when the underlying conditions are raining (which is assumed to be comparable to the convective situation here where the condensate would fall as rain if precipitation processes were included). Additionally, q_c has developed strong positive covariances with w in this region, indicating that condensate has formed, with larger amounts with increasing w. As indicated in Fig. 2 (contours in right panels), the condensate concentrations increase with height up to ∼9 km. This suggests that the region of positive q_c–w covariances just mentioned is not caused by advection (the covariances would have a negative sign if they were) but instead by the process of more condensate being produced in association with stronger w. Note that in most parts of the cloud layer, q and q_c are both positively correlated with w, which may appear counter-intuitive. In the convection zone, convergence and condensation processes are happening. The convergence can increase q while the condensation will decrease q but increase q_c. As long as the increase in q caused by convergence is larger than the loss of q by condensation, there is a net increase in q, in which case both q and q_c could have a positive correlation with the vertical wind. In fact statistics show that often q and q_c are mutually positively correlated at points inside the cloud (not shown).

At 160 min a shallow surface layer of negative q_c–w covariances has developed, which may be due to a circulation being induced by the convection, which is drawing in air with lower q_c from the surroundings. There is an additional region of weak, but non-zero, w–w and u–w correlations above the top of the convection, the latter indicating a further (but weak and shallow) shear region at the higher levels, which continues to 200 min. From 240 min the covariances develop more complicated, but weaker, patterns, including the development of a further possible u shear layer near the surface. There is also an influence of the convection in vertical covariances at locations other than at x=148.5 km, including at locations within and beyond the moist “bubble” set-up in the initial conditions (covariances not shown, but correlations are shown below).

Figure 6 Evolution of a selection of vertical correlations at times between 40 and 360 min. The top set of three rows is for ${\stackrel{\mathrm{̃}}{\mathit{\rho }}}^{\prime }$–${\stackrel{\mathrm{̃}}{\mathit{\rho }}}^{\prime }$ correlations, and the bottom set is for b^′–${\stackrel{\mathrm{̃}}{\mathit{\rho }}}^{\prime }$ correlations. For each set of three rows, the rows are in the convection zone (x=148.5 km – the same position as used for Fig. 5), the moist zone (286.5 km), and the dry zone (421.5 km) respectively. In each plot of p₁–p₂ correlations, the y axis is for p₁, the x axis is for p₂, and the surface is at the bottom left of each square. Red (blue) indicates positive (negative) correlations. The circular and triangular symbols represent points where consistency with hydrostatic balance is tested. These are placed at a selection of heights in the top three rows where there is a local maximum or minimum in ${\stackrel{\mathrm{̃}}{\mathit{\rho }}}^{\prime }$–${\stackrel{\mathrm{̃}}{\mathit{\rho }}}^{\prime }$ correlation with respect to the y axis. Each symbol is also plotted in the corresponding bottom three rows to see if it associated with a “small” b^′–${\stackrel{\mathrm{̃}}{\mathit{\rho }}}^{\prime }$ correlation. Circles are associated with small b^′–${\stackrel{\mathrm{̃}}{\mathit{\rho }}}^{\prime }$ correlations (≤0.3, associated with hydrostatic balance), and triangles are given otherwise. Some symbols are labelled with letters, which are referenced in the text.

Figure 6 shows vertical correlations between a different selection of quantities to those in Fig. 5. The top three rows show auto-correlations of ${\stackrel{\mathrm{̃}}{\mathit{\rho }}}^{\prime }$ for the convection, moist, and dry zones respectively (see Fig. 2 for the definition of the zones where the columns are selected and Fig. 6 caption for their lateral positions). At t=40 min the different zones show similar patterns. The model variable ${\stackrel{\mathrm{̃}}{\mathit{\rho }}}^{\prime }$ appears to comprise two layers of correlations (one below 6 km and the other above), where in each, ${\stackrel{\mathrm{̃}}{\mathit{\rho }}}^{\prime }$ is strongly correlated vertically but only weakly (and negatively) across each layer. In the convection zone there are some small-scale vertical oscillations, which are not present in the other zones shown and so are likely to be due to early signs of convection. At 80 min in the convection zone, the lower layer level itself splits into two negatively correlated sub-layers, which are associated with the establishment of upper- and lower-pressure changes formed in connection with the u–w covariances in Fig. 5 at this time (the wind shearing). The covariances in the moist and dry zones have not changed qualitatively very much at this time. At 120 and 160 min in the convection zone the sub-layering has disappeared (probably because the shear zone has moved upwards), and the vertical length scales above 6 km have dramatically reduced. This shortening of the vertical length scales is also seen in the moist and dry zones but with a delay of tens of minutes. From 200 min onwards the vertical covariances in the convection zone do not change very much in structure, but the other zones develop strong large-scale (vertically) correlation structures presumably due to outward propagation of disturbances from the convection zone.

The bottom three rows of Fig. 6 show vertical correlations between b^′ at different heights (y axis) and ${\stackrel{\mathrm{̃}}{\mathit{\rho }}}^{\prime }$ at different heights (x axis). These are included in combination with the ${\stackrel{\mathrm{̃}}{\mathit{\rho }}}^{\prime }$ auto-correlations to look for signatures of hydrostatic balance (HB). The HB equation is found from Eq. (1c) after the Lagrangian derivative terms are neglected, resulting in $C\partial {\stackrel{\mathrm{̃}}{\mathit{\rho }}}^{\prime }/\partial z=b^{\prime }$. This relates b^′ to vertical gradients in ${\stackrel{\mathrm{̃}}{\mathit{\rho }}}^{\prime }$. Thus the property of HB is assumed when the vertical derivative of the ${\stackrel{\mathrm{̃}}{\mathit{\rho }}}^{\prime }\text{–}{\stackrel{\mathrm{̃}}{\mathit{\rho }}}^{\prime }$ correlation function at a point in the top three rows of Fig. 6 is zero and is matched by the $b^{\prime }\text{–}{\stackrel{\mathrm{̃}}{\mathit{\rho }}}^{\prime }$ correlation also being zero at the corresponding point in the bottom three panels. This property is seen at many points in Fig. 6, and a sample of points is highlighted in the figure and below.

In the top three rows, symbols are placed at a selection of points where the vertical correlation is stationary ($\partial \left({\stackrel{\mathrm{̃}}{\mathit{\rho }}}^{\prime }\text{–}{\stackrel{\mathrm{̃}}{\mathit{\rho }}}^{\prime }\text{ correlation}\right)/$ $\partial z\left({\stackrel{\mathrm{̃}}{\mathit{\rho }}}^{\prime }\right)=\mathrm{0}$). The corresponding symbols are placed at the same positions in the bottom three panels. A symbol is a circle if it has a small (proxy for zero) b^′–${\stackrel{\mathrm{̃}}{\mathit{\rho }}}^{\prime }$ correlation magnitude (≤0.3), which is associated with HB, or a triangle otherwise. At 80 min for instance in the moist zone, there are four stationary points highlighted (A, B, C, and D). These correspond to near-zero values of b^′–${\stackrel{\mathrm{̃}}{\mathit{\rho }}}^{\prime }$ correlations (a, b, c, and d). There are only a few points in the moist and dry zones in the early stages of the evolution that do not satisfy this condition (triangles). At 160 min in the convection zone, there are two stationary points highlighted (E and F). These do not correspond to near-zero values of b^′–${\stackrel{\mathrm{̃}}{\mathit{\rho }}}^{\prime }$ correlations (e and f), suggesting a clear breakdown of HB. In the first half of the forecast period, there are more non-hydrostatic signatures in the convection zone than in the other zones. Later in the convection zone, the hydrostatic signatures return as the system stabilises. In the second half of the forecast period, the moist and then dry zones experience more non-hydrostatic signatures as a result of the outward propagation of disturbances from the convection zone. At 320 min in the dry zone for instance, there are two stationary points highlighted (G and H). These do not correspond to near-zero values of b^′–${\stackrel{\mathrm{̃}}{\mathit{\rho }}}^{\prime }$ correlations (g and h), suggesting a propagated breakdown of HB. Bannister et al. (2011) and Vetra-Carvalho et al. (2012) similarly found evidence of non-hydrostatic patterns in the covariance patterns in rainy regions (compared to dry regions) of forecast errors in the convective-scale UM. These statistics encourage the development of data assimilation systems that do not always constrain the background error covariances to be hydrostatic when convection is happening. Hydrostatic imbalance is further mentioned in Sect. 4.2.

4.2 Possible indicators and harbingers of convection

Since there are significant differences in the covariance and correlation structures for convecting and non-convecting regions, convective-scale data assimilation systems should gain advantages by being aware of the convection “status” of any geographical region and thus should use appropriate background error statistics. This is done automatically in ensemble-based systems (e.g. ensemble Kalman filters or ensemble-variational methods; Bannister, 2017) and partially in hybrid systems (where the static background error covariances of traditional variational methods are averaged with ensemble-derived covariances; Bannister, 2017). It is well known though that the use of ensemble-based covariances leads to rank deficiency and under-sampling (Houtekamer and Mitchell, 1998; van Leeuwen, 1999; Hamill et al., 2001; Houtekamer and Mitchell, 2005; Ehrendorfer, 2007; Zhang and Zhang, 2012; Houtekamer and Zhang, 2016), and so any progress that allows the (otherwise) static background error covariance matrix to gain flow dependence is still valuable. This is the strategy behind a number of studies, including Montmerle and Berre (2010), Ménétrier and Montmerle (2011), Michel et al. (2011), Montmerle (2012), and Yang et al. (2022), although these studies used the strength of precipitation or the presence of fog (rather than convection) as a criterion for selection of the most appropriate covariance statistics. Here, we explore ways to determine locations of convecting and non-convecting flows. Our strategy is to determine these regions ideally before convection has fully developed so that a data assimilation system that uses the right covariances with new observations does not lead to an analysis state that will impede the development of convection in the ensuing forecasts by imposing HB wrongly.

Seven quantities are considered here to help in the discrimination. (i) Convective available potential energy (CAPE) is a vertically integrated quantity, which can be used to help determine if an atmospheric column is “primed” for convection. (ii) Convective inhibition (CIN) is also a vertically integrated quantity, which can be used to help determine if a column is “resistant” to convection. CAPE and CIN are derived from Hydro-ABC quantities in Appendix C. (iii) RH at any horizontal location is defined as the maximum value of $q/q_{\mathrm{s}}$ in the vertical column, where q_s is the saturated specific humidity defined in Sect. 2.4. Large values of RH are expected during and just before moist convection. The column maximum (iv) vertical wind, w^max, and (v) condensate, $q_{\mathrm{c}}^{\max }$, are important signatures of active convection but not necessarily of pre-convective conditions. (vi) Hydrostatic imbalance (HI) at any horizontal location is defined as the largest absolute value of $\left(C\partial {\stackrel{\mathrm{̃}}{\mathit{\rho }}}^{\prime }/\partial z-b^{\prime }\right)/$ $\left(\mathrm{rms}\left(C\partial {\stackrel{\mathrm{̃}}{\mathit{\rho }}}^{\prime }/\partial z\right)+\mathrm{rms}\left(b^{\prime }\right)\right)$ in the column, where rms is the root mean square calculated over the whole domain. A value of zero represents perfect hydrostatic balance, so convection is associated with non-zero HI. (vii) The maximum value of $\partial u/\partial x$ in a column (horizontal divergence, HD – recall there is no y dependence in ABC) is the last quantity considered.

Figure 7 shows how these quantities vary with horizontal position in the domain at t=10 min, which is before convection has developed. CAPE (red, top panel) increases towards the convection zone but is not a smooth function of position and drops slightly at the future convection point. CAPE increases abruptly at ∼0 km and drops at ∼285 km, even though these points are within the region where q>0 in the initial conditions. CIN (blue, top) also changes abruptly at these positions but in the opposite way: CIN is small where CAPE is large and vice versa. RH (green, top) in this case is a smooth function of position, peaking at unity around the convection zone. While CAPE, CIN, and RH give physically reasonable indications of where convection is possible, they do not pinpoint where convection develops (namely in the convection zone).

w^max (red, bottom), $q_{\mathrm{c}}^{\max }$ (magenta, bottom), and HI (green, bottom) are fields with a finer scale than CAPE, CIN, and RH, and they peak at the precise location where convection develops. The HD (blue, bottom) is also fine-scale but does not peak at the correct location of the future convection. The quantities w^max, $q_{\mathrm{c}}^{\max }$, and HI are therefore promising harbingers of convection, which may be useful for selecting vertical background error statistics that are characteristic of convection (rather than quiescent conditions) for data assimilation. This is the case especially when the data assimilation method is otherwise non-flow-dependent, as in traditional variational schemes. Further work will be required to test the reliability of these indicators.

Figure 7 The geographical variation in various indicators/harbingers of convection at t=10 min. Plotted in (a) are convective available potential energy (red line), convective inhibition (blue), and relative humidity (green). Plotted in (b) are column maximum values of vertical wind (red), the condensate mixing ratio (magenta), hydrostatic imbalance (green), and horizontal divergence (blue). The vertical black lines mark locations that are within the convection, moist, and dry zones (see Fig. 2).

5 Conclusions

Described in this paper is how the simplified ABC model of convective-scale flow has been extended to include moist processes, namely evaporation and condensation, yielding a revised model called Hydro-ABC. This is the next step in the aim to provide a low-cost framework to study convective-scale data assimilation (DA) strategies. The following new variables, terms, and parameters have been introduced.

• The model's vapour mixing ratio field is q (the gaseous phase of water), Eq. (1g).

• The model's condensate mixing ratio field is q_c (the solid and/or liquid phases of water), Eq. (1h).

• The net rate of the source of vapour (or equivalently the net rate of the sink of condensate) is Ev−Co (the evaporation minus the condensation rates) for each grid point. Only the difference needs to be known and not the values of Ev and Co separately.

• Latent heat is released (absorbed) for net condensation (evaporation). The latent heating rate is ${\mathcal{S}}_{{\mathrm{b}}^{\prime }}$ and is a source term for buoyancy (Eq. 1e) (latent heat is assumed to be exchanged directly only with buoyancy energy). ${\mathcal{S}}_{{\mathrm{b}}^{\prime }}$ and Ev−Co are not determined explicitly but are the result of potentially solving a non-linear system of equations as set out in Sect. 2.3. The latent heat of vaporisation, L_v (set to 2500 J g⁻¹), is introduced and relates changes in vapour to changes in latent heat energy, Eq. (3).

• Condensation happens when situations are super-saturated (q>q_s), in which case q drops immediately (practically over a model time step) to a new q_s value, and b^′ is increased from its previous value due to latent heat release. This is the reason for the non-linear equations as mentioned above; see Eq. (8). When $b^{\prime }<\mathrm{0}$, condensation is not allowed when the ratio of latent heat release to buoyancy energy is greater than Γ (set to 10). This is to avoid excessively large jumps in buoyancy.

• Evaporation happens when situations are sub-saturated, when there is a finite amount of condensate to evaporate, and when b^′ is positive. The last condition is imposed as otherwise b^′ will increase (i.e. become less negative) when evaporation happens, which is unphysical. Evaporation happens at a slower rate (timescale τ, set to 1000 s) than condensation; see Eq. (7).

• The saturated mixing ratio, q_s, sets a maximum q that can be supported. In real systems q_s is a function of temperature and pressure. Since these quantities do not exist in the ABC model, analogous quantities are derived based on an assumed reference density profile, hydrostatic balance, and the correspondence between buoyancy and potential temperature. A formula is derived for q_s in terms of ABC variables b^′ and ${\stackrel{\mathrm{̃}}{\mathit{\rho }}}^{\prime }$ as set out in Sect. 2.4. The conditions used to relate b^′ and ${\stackrel{\mathrm{̃}}{\mathit{\rho }}}^{\prime }$ to q_s do not need to be exactly relevant to the ABC system but are used only as a basis for a sensible relationship.

Hydro-ABC is constructed to satisfy certain conservation properties over the whole domain. These are for total water (q+q_c, Sect. 2.2.1) and total moist energy (defined as dry energy plus latent heat energy, Sect. 2.2.2). There are no precipitation processes in Hydro-ABC; this is to keep the system as simple as possible while allowing highly non-linear processes.

An example run of Hydro-ABC shows a strong convective event driven by latent heat release. This forms a realistic anvil cloud (Sect. 3.1) and the excitation of inertio-gravity and acoustic modes, higher in frequency than those normally generated by the dry ABC model (Sect. 3.3).

A very important concern in DA applied to convective-scale flows involves how the background error covariance statistics are affected by the flow situation. This is an issue in real forecasting and DA systems where attempts have been made in the literature to use statistics that depend on the presence of background properties like rain and fog. This study gives an example of this by demonstrating how flow-dependent statistics (found from a 93-member ensemble) can depend on the presence of convection. It is found that there can be significant differences between the covariances in convecting and non-convecting regions, including the breakdown of hydrostatic balance. Many operational systems base their static background error covariance model, or the training data used to calibrate it, on hydrostatic balance (Berre, 2000; Bannister, 2008; Gustafsson et al., 2018; Bannister et al., 2020), and so this result may urge changes to this practice.

Given that the covariances depend on convection, it is necessary to be able to diagnose reliably the presence of convection (or impending convection), especially if a suitable ensemble is not available and one needs to instead prescribe covariances for the DA. Such an indicator can be applied to the background state in a DA application in order to decide the most relevant set of error covariances. Using climatologically derived (static) error covariances, or covariances that support hydrostatic flow, may well impede the progress of convection prediction. A number of convection diagnostics have been proposed here as an indicator or harbinger of convection, namely relative humidity, hydrostatic imbalance, horizontal divergence, convective available potential energy, convective inhibition, vertical wind, and the condensate mixing ratio (Sect. 4.2). It is found in the example shown that all of these quantities are sensitive to convection, but hydrostatic imbalance, vertical wind, and the condensate mixing ratio seem to be the best indicators of pinpointing impending convection. Further work though is required to decide which quantity or combination of quantities is best and which threshold values should be used. This will be part of the next stage of this work, namely to adapt the dry ABC–DA system to assimilate water-related information. Another interesting avenue that may be explored is to add a precipitation process (i.e. an additional fast downward transport of condensate), which could give rise to cold pools.

Appendix A: Continuous total energy equations of Hydro-ABC

The formulation of Hydro-ABC is based on a time-step split into a dynamical core and a micro-physics parametrisation. Hydro-ABC can also be represented as continuous equations, i.e. as Eqs. (1a–h) but with Eq. (1e) replaced with

$$\begin{array}{}\text{(A1)}& {\displaystyle \frac{\partial b^{\prime }}{\partial t}}+B\mathit{u}\cdot \mathrm{\nabla }{b}^{\prime }+A^{\mathrm{2}}w=-{\displaystyle \frac{A^{\mathrm{2}}{L}_{\mathrm{v}}}{b^{\prime }}}(\text{Ev}-\text{Co}),\end{array}$$

which is an explicit form of the buoyancy equation. This is consistent with the basis of our micro-physics parametrisation, Eq. (6), where ${\mathrm{\Delta }}^{\mathrm{mp}}\left({b^{\prime }}^{\mathrm{2}}\right)=\frac{\partial }{\partial t}\left({b^{\prime }}^{\mathrm{2}}\right)\mathrm{\Delta }t=\mathrm{2}{b}^{\prime }\frac{\partial }{\partial t}\left(b^{\prime }\right)\mathrm{\Delta }t=\mathrm{2}{b}^{\prime }{\mathrm{\Delta }}^{\mathrm{mp}}{b}^{\prime }$. Plugging this into Eq. (6) gives

$${\displaystyle }{\mathrm{\Delta }}^{\mathrm{mp}}{b}^{\prime }=-{\displaystyle \frac{A^{\mathrm{2}}{L}_{\mathrm{v}}}{b^{\prime }}}\left(\text{Ev}-\text{Co}\right)\mathrm{\Delta }t,$$

i.e. the right-hand side of Eq. (A1). In order to show that this leads to an energy-conserving system of equations, multiply Eq. (A1) by $\stackrel{\mathrm{̃}}{\mathit{\rho }}{b}^{\prime }/A^{\mathrm{2}}$:

$$\begin{array}{c}{\displaystyle \frac{\stackrel{\mathrm{̃}}{\mathit{\rho }}{b}^{\prime }}{A^{\mathrm{2}}}}{\displaystyle \frac{\partial b^{\prime }}{\partial t}}+B{\displaystyle \frac{\stackrel{\mathrm{̃}}{\mathit{\rho }}{b}^{\prime }}{A^{\mathrm{2}}}}\mathit{u}\cdot \mathrm{\nabla }{b}^{\prime }+\stackrel{\mathrm{̃}}{\mathit{\rho }}{b}^{\prime }w=-\stackrel{\mathrm{̃}}{\mathit{\rho }}{L}_{\mathrm{v}}(\text{Ev}-\text{Co})\\ {\displaystyle }\begin{array}{rl}& \stackrel{\mathrm{̃}}{\mathit{\rho }}{\displaystyle \frac{\partial \left[{b^{\prime }}^{\mathrm{2}}/\left(\mathrm{2}{A}^{\mathrm{2}}\right)\right]}{\partial t}}+B\stackrel{\mathrm{̃}}{\mathit{\rho }}\mathit{u}\cdot \mathrm{\nabla }\left[{b^{\prime }}^{\mathrm{2}}/\left(\mathrm{2}{A}^{\mathrm{2}}\right)\right]+\stackrel{\mathrm{̃}}{\mathit{\rho }}{b}^{\prime }w\\ & =-\stackrel{\mathrm{̃}}{\mathit{\rho }}{L}_{\mathrm{v}}(\text{Ev}-\text{Co}).\end{array}\end{array}$$

Then use the following identity (Petrie et al., 2017, their Eq. 20):

$$\begin{array}{}\text{(A2)}& \stackrel{\mathrm{̃}}{\mathit{\rho }}\left({\displaystyle \frac{\partial \mathit{\gamma }}{\partial t}}+B\mathit{u}\cdot \mathrm{\nabla }\mathit{\gamma }\right)={\displaystyle \frac{\partial \left(\stackrel{\mathrm{̃}}{\mathit{\rho }}\mathit{\gamma }\right)}{\partial t}}+B\mathrm{\nabla }\cdot \left(\stackrel{\mathrm{̃}}{\mathit{\rho }}\mathit{\gamma }\mathit{u}\right),\end{array}$$

with $\mathit{\gamma }={b^{\prime }}^{\mathrm{2}}/\left(\mathrm{2}{A}^{\mathrm{2}}\right)$. That is,

$$\begin{array}{ll}{\displaystyle }& {\displaystyle \frac{\partial \left[\stackrel{\mathrm{̃}}{\mathit{\rho }}{b^{\prime }}^{\mathrm{2}}/\left(\mathrm{2}{A}^{\mathrm{2}}\right)\right]}{\partial t}}+B\mathrm{\nabla }\cdot \left[\stackrel{\mathrm{̃}}{\mathit{\rho }}{b^{\prime }}^{\mathrm{2}}/\left(\mathrm{2}{A}^{\mathrm{2}}\right)\mathit{u}\right]+\stackrel{\mathrm{̃}}{\mathit{\rho }}{b}^{\prime }w\\ {\displaystyle }& {\displaystyle }=-\stackrel{\mathrm{̃}}{\mathit{\rho }}{L}_{\mathrm{v}}(\text{Ev}-\text{Co}).\end{array}$$

Making the substitution $E_{\mathrm{b}}=\stackrel{\mathrm{̃}}{\mathit{\rho }}{b^{\prime }}^{\mathrm{2}}/\left(\mathrm{2}{A}^{\mathrm{2}}\right)$ (see Sect. 2.2.2), gives

$$\begin{array}{}\text{(A3)}& {\displaystyle \frac{\partial E_{\mathrm{b}}}{\partial t}}+B\mathrm{\nabla }\cdot \left[E_{\mathrm{b}}\mathit{u}\right]+\stackrel{\mathrm{̃}}{\mathit{\rho }}{b}^{\prime }w=-\stackrel{\mathrm{̃}}{\mathit{\rho }}{L}_{\mathrm{v}}(\text{Ev}-\text{Co}).\end{array}$$

This is the continuous form of the equation for buoyancy energy. Next multiply Eq. (1g) by L_v,

$$\begin{array}{}\text{(A4)}& {\displaystyle \frac{\partial \stackrel{\mathrm{̃}}{\mathit{\rho }}{L}_{\mathrm{v}}q}{\partial t}}+B\mathrm{\nabla }\cdot \left(\stackrel{\mathrm{̃}}{\mathit{\rho }}{L}_{\mathrm{v}}q\mathit{u}\right)=\stackrel{\mathrm{̃}}{\mathit{\rho }}{L}_{\mathrm{v}}(\text{Ev}-\text{Co}),\end{array}$$

and make the substitution from Eq. (3),

$$\begin{array}{}\text{(A5)}& {\displaystyle \frac{\partial E_{\mathrm{l}}}{\partial t}}+B\mathrm{\nabla }\cdot \left(E_{\mathrm{l}}\mathit{u}\right)=\stackrel{\mathrm{̃}}{\mathit{\rho }}{L}_{\mathrm{v}}(\text{Ev}-\text{Co}).\end{array}$$

This is the continuous form of the equation for latent heat energy. The equations for kinetic and elastic energies (E_k and E_e respectively) are derived in a similar way to Eq. (A3) (see Petrie et al., 2017, their Sect. 2.3):

$$\begin{array}{}\text{(A6)}& {\displaystyle \frac{\partial E_{\mathrm{k}}}{\partial t}}+B\mathrm{\nabla }\cdot \left(E_{\mathrm{k}}\mathit{u}\right)-\stackrel{\mathrm{̃}}{\mathit{\rho }}{b}^{\prime }w+C\stackrel{\mathrm{̃}}{\mathit{\rho }}\mathit{u}\cdot \mathrm{\nabla }{\stackrel{\mathrm{̃}}{\mathit{\rho }}}^{\prime }=\mathrm{0},\text{(A7)}& {\displaystyle \frac{\partial E_{\mathrm{e}}}{\partial t}}+C{\stackrel{\mathrm{̃}}{\mathit{\rho }}}^{\prime }\mathrm{\nabla }\cdot \left(\stackrel{\mathrm{̃}}{\mathit{\rho }}\mathit{u}\right)=\mathrm{0}.\end{array}$$

Adding Eqs. (A3), (A5), (A6), and (A7) leads to the equation for the total moist energy, $E_{\mathrm{moist}}=E_{\mathrm{k}}+E_{\mathrm{b}}+E_{\mathrm{e}}+E_{\mathrm{l}}$:

$${\displaystyle \frac{\partial E_{\mathrm{moist}}}{\partial t}}+B\mathrm{\nabla }\cdot \left[\left(E_{\mathrm{b}}+E_{\mathrm{k}}+E_{\mathrm{l}}\right)\mathit{u}\right]+C\mathrm{\nabla }\cdot \left({\stackrel{\mathrm{̃}}{\mathit{\rho }}}^{\prime }\stackrel{\mathrm{̃}}{\mathit{\rho }}\mathit{u}\right)=\mathrm{0}.$$

Integrating this quantity over the whole domain, using the divergence theorem, and exploiting the horizontal periodic boundary conditions, as well as the zero vertical wind at the top and bottom boundaries (see Petrie et al., 2017, their Sect. 3.2), the following is found:

$$\begin{array}{}\text{(A8)}& {\displaystyle \frac{\partial }{\partial t}}\left(\int \int \mathrm{d}x\mathrm{d}z{E}_{\mathrm{moist}}\right)=\mathrm{0}.\end{array}$$

This proves that the integral of total moist energy is conserved under the governing equations, Eqs. (1) and (A1).

The issue with Eq. (A1) is the $\mathrm{1}/b^{\prime }$ factor on the right-hand side. This issue is avoided by using the micro-physics parametrisation step represented by Eq. (6).

Appendix B: Pseudo-code for the micro-physics parametrisation

Annotated pseudo-code describing the dynamics/micro-physics parametrisation is given in this part of the Appendix.

• 1.

  Propagate the Hydro-ABC fields u, v, w, ${\stackrel{\mathrm{̃}}{\mathit{\rho }}}^{\prime }$, b^′, q, and q_c (described by Eq. 1) by one time step.

• 2.

  Calculate the pre-micro-physics saturated specific humidity field, q_s(x,z) (Eq. 10). This requires the corresponding pressure, p(x,z) (Eq. 11), and temperature, T(x,z) (Eq. 13), which are related to model variables. How the micro-physics behaves depends on Eq. (7).

• 3.

  If $q(x,z)>q_{\mathrm{s}}(x,z)$, condensation may happen (see points 1a and 1b in Sect. 2.3). The following are performed for each x and z.

  • a.

    If $b^{\prime }\ge \mathrm{0}$, allow condensation to happen by solving the following non-linear simultaneous equations,

    $$\begin{array}{c}{\displaystyle }{b}_{\mathrm{mp}}^{\prime }=+\sqrt{{b^{\prime }}^{\mathrm{2}}+\mathrm{2}{A}^{\mathrm{2}}{L}_{\mathrm{v}}\left(q-q_{\mathrm{mp}}\right)}\text{and}\\ {\displaystyle }{q}_{\mathrm{mp}}=q_{\mathrm{s}}\left(b_{\mathrm{mp}}^{\prime }\right),\end{array}$$

    and then setting ${q_{\mathrm{c}}}_{\mathrm{mp}}=q_{\mathrm{c}}+\left(q-q_{\mathrm{mp}}\right).$

  • b.

    Else if $\stackrel{\mathrm{̃}}{\mathit{\rho }}{L}_{\mathrm{v}}\left|q-q_{\mathrm{s}}\right|/E_{\mathrm{b}}>\mathrm{\Gamma }$, then the latent heat energy released by condensation will be much greater than the current buoyancy energy. This will allow condensation to happen by setting

    $$\begin{array}{c}{\displaystyle }{b}_{\mathrm{mp}}^{\prime }=+\sqrt{{b^{\prime }}^{\mathrm{2}}+\mathrm{2}{A}^{\mathrm{2}}{L}_{\mathrm{v}}\left(q-q_{\mathrm{s}}\right)},\\ {\displaystyle }{q}_{\mathrm{mp}}=q_{\mathrm{s}},\\ {\displaystyle }{q_{\mathrm{c}}}_{\mathrm{mp}}=q_{\mathrm{c}}+\left(q-q_{\mathrm{mp}}\right).\end{array}$$
  • c.

    Else do nothing (to avoid unrealistic situations).

• 4.

  Else if $b^{\prime }>\mathrm{0}$, evaporation may happen (see point 2 in Sect. 2.3 and Eq. 7). The following are performed for each x and z. Calculate three positive semi-definite quantities ($\left(q_{\mathrm{s}}-q\right)\left(\mathrm{1}-e^{-\mathrm{\Delta }t/\mathit{\tau }}\right)$, q_c, and ${b^{\prime }}^{\mathrm{2}}/\left(\mathrm{2}{A}^{\mathrm{2}}{L}_{\mathrm{v}}\right)$) and find which has the smallest value.

  • a.

    If $(q_{\mathrm{s}}-q)\left(\mathrm{1}-e^{-\mathrm{\Delta }t/\mathit{\tau }}\right)$ is the smallest, then there is enough condensate and buoyancy energy in the grid box to allow this amount of evaporation to occur by setting

    $$\begin{array}{c}{\displaystyle }{b}_{\mathrm{mp}}^{\prime }=+\sqrt{{b^{\prime }}^{\mathrm{2}}-\mathrm{2}{A}^{\mathrm{2}}{L}_{\mathrm{v}}\left(q_{\mathrm{s}}-q\right)\left(\mathrm{1}-e^{-\mathrm{\Delta }t/\mathit{\tau }}\right)},\\ {\displaystyle }{q_{\mathrm{c}}}_{\mathrm{mp}}=q_{\mathrm{c}}-\left(q_{\mathrm{s}}-q\right)\left(\mathrm{1}-e^{-\mathrm{\Delta }t/\mathit{\tau }}\right),\\ {\displaystyle }{q}_{\mathrm{mp}}=q+\left(q_{\mathrm{s}}-q\right)\left(\mathrm{1}-e^{-\mathrm{\Delta }t/\mathit{\tau }}\right).\end{array}$$
  • b.

    If q_c is the smallest, then evaporation is limited only by availability of condensate (there is enough buoyancy energy available to allow all of this condensate in the grid box to evaporate) by transferring all condensate to vapour, leaving

    $$\begin{array}{c}{\displaystyle }{b}_{\mathrm{mp}}^{\prime }=+\sqrt{{b^{\prime }}^{\mathrm{2}}-\mathrm{2}{A}^{\mathrm{2}}{L}_{\mathrm{v}}{q}_{\mathrm{c}}},\\ {\displaystyle }{q_{\mathrm{c}}}_{\mathrm{mp}}=\mathrm{0},\\ {\displaystyle }{q}_{\mathrm{mp}}=q+q_{\mathrm{c}}.\end{array}$$
  • c.

    Else if ${b^{\prime }}^{\mathrm{2}}/\left(\mathrm{2}{A}^{\mathrm{2}}{L}_{\mathrm{v}}\right)$ is the smallest, then evaporation is limited only by availability of buoyancy energy, $E_{\mathrm{b}}={b^{\prime }}^{\mathrm{2}}/\left(\mathrm{2}{A}^{\mathrm{2}}\stackrel{\mathrm{̃}}{\mathit{\rho }}\right)$ (there is enough condensate available), by transferring all buoyancy energy to latent heat, leaving

    $$\begin{array}{c}{\displaystyle }{b}_{\mathrm{mp}}^{\prime }=\mathrm{0},\\ {\displaystyle }{q_{\mathrm{c}}}_{\mathrm{mp}}=q_{\mathrm{c}}-{b^{\prime }}^{\mathrm{2}}/\left(\mathrm{2}{A}^{\mathrm{2}}{L}_{\mathrm{v}}\right),\\ {\displaystyle }{q}_{\mathrm{mp}}=q+{b^{\prime }}^{\mathrm{2}}/\left(\mathrm{2}{A}^{\mathrm{2}}{L}_{\mathrm{v}}\right).\end{array}$$
• 5.

  Overwrite the model's variables, b^′, q, and q_c, with their post-micro-physics values, $b_{\mathrm{mp}}^{\prime }$, q_mp, and ${q_{\mathrm{c}}}_{\mathrm{mc}}$ respectively.

Appendix C: Computation of the convective available potential energy (CAPE)

In the real atmosphere, CAPE is proportional to the vertical integral of the difference between a moist parcel's temperature, $\widehat{T}$, and the environmental temperature, T_env (e.g. Sect. 7.4.1 of Salby, 1996; $\mathrm{CAPE}=g{\int }_{\mathrm{LFC}}^{\mathrm{LNB}}\left(\widehat{T}-T_{\mathrm{env}}\right)/T_{\mathrm{env}}\mathrm{d}z$, where LFC is the level of free convection and LNB is the level of neutral buoyancy – see below). CAPE represents the amount of energy that can be converted to kinetic energy by convective instability. A positive value indicates energy is available for convection and is thus a measure of instability. The above equation for CAPE is relevant to the standard equations of motion. In Hydro-ABC there is no explicit temperature variable, but we can derive a physically reasonable CAPE-like quantity by analysing the energetics of an air parcel using the ABC equations.

Consider an air parcel, whose properties are indicated with a “hat”. The Lagrangian derivative of $\widehat{w}$ is

$$\begin{array}{}\text{(C1)}& {\displaystyle \frac{\partial \widehat{w}}{\partial t}}+B\widehat{\mathit{u}}\cdot \mathrm{\nabla }\widehat{w}={\displaystyle \frac{D\widehat{w}}{\mathrm{D}t}}={\displaystyle \frac{D^{\mathrm{2}}\widehat{z}}{\mathrm{D}{t}^{\mathrm{2}}}},\end{array}$$

where $\widehat{z}$ is the vertical position of the parcel. Equation (C1) is the vertical acceleration of the parcel. We can substitute for $\partial \widehat{w}/\partial t+B\widehat{\mathit{u}}\cdot \mathrm{\nabla }\widehat{w}$ using the vertical momentum Eq. (1c). The parcel equation is then

$$\begin{array}{}\text{(C2)}& {\displaystyle \frac{D^{\mathrm{2}}\widehat{z}}{\mathrm{D}{t}^{\mathrm{2}}}}=-C{\displaystyle \frac{\partial \widehat{\stackrel{\mathrm{̃}}{{\mathit{\rho }}^{\prime }}}}{\partial z}}+\widehat{b^{\prime }}.\end{array}$$

The environmental profile (the mean state around the convecting column, found by averaging the model fields over 101 points horizontally centred on the profile of interest) is found to be in hydrostatic equilibrium (not shown):

$$\begin{array}{}\text{(C3)}& -C{\displaystyle \frac{\partial {\stackrel{\mathrm{̃}}{\mathit{\rho }}}_{\mathrm{env}}^{\prime }}{\partial z}}+b_{\mathrm{env}}^{\prime }=\mathrm{0}.\end{array}$$

Taking the difference between Eqs. (C2) and (C3) gives

$$\begin{array}{}\text{(C4)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{D^{\mathrm{2}}\widehat{z}}{\mathrm{D}{t}^{\mathrm{2}}}}& =-C{\displaystyle \frac{\partial \left(\widehat{\stackrel{\mathrm{̃}}{{\mathit{\rho }}^{\prime }}}-{\stackrel{\mathrm{̃}}{\mathit{\rho }}}_{\mathrm{env}}^{\prime }\right)}{\partial z}}+\left(\widehat{b^{\prime }}-b_{\mathrm{env}}^{\prime }\right)\\ & \approx \widehat{b^{\prime }}-b_{\mathrm{env}}^{\prime },\end{array}\end{array}$$

where in the second line we have assumed mechanical equilibrium (that $\widehat{\stackrel{\mathrm{̃}}{{\mathit{\rho }}^{\prime }}}={\stackrel{\mathrm{̃}}{\mathit{\rho }}}_{\mathrm{env}}^{\prime }$). This gives a formula for the parcel's acceleration (or specific force), which is used in the CAPE formula below.

The parcel's buoyancy profile, $\widehat{b^{\prime }}\left(z_i\right)$, is found by raising a hypothetical parcel from the surface. At the surface the parcel's buoyancy and vapour mixing ratio take the model's values, $\widehat{b^{\prime }}\left(z_{\mathrm{0}}\right)=b^{\prime }\left(z_{\mathrm{0}}\right)$ and $\widehat{q}\left(z_{\mathrm{0}}\right)=q\left(z_{\mathrm{0}}\right)$ respectively. Start with level i=1. By raising the parcel to the next level, the total buoyancy is assumed to be conserved: $\widehat{b}\left(z_i\right)=\widehat{b}\left(z_{i-\mathrm{1}}\right)$. The total buoyancy is the sum of the reference state and the perturbation, namely $\widehat{b}\left(z_i\right)={\widehat{b}}_{\mathrm{0}}\left(z_i\right)+\widehat{b^{\prime }}\left(z_i\right)$. By definition (Petrie et al., 2017), A² is the rate of increase in reference buoyancy with height, namely ${\widehat{b}}_{\mathrm{0}}\left(z_{i+\mathrm{1}}\right)-{\widehat{b}}_{\mathrm{0}}\left(z_i\right)=A^{\mathrm{2}}(z_{i+\mathrm{1}}-z_i)$ in discrete form. Putting the last three equations together gives the parcel's perturbation value at z_i:

$$\begin{array}{}\text{(C5)}& \widehat{b^{\prime }}\left(z_i\right)=\widehat{b^{\prime }}\left(z_{i-\mathrm{1}}\right)-A^{\mathrm{2}}(z_i-z_{i-\mathrm{1}}).\end{array}$$

The parcel's mixing ratio is also assumed to be conserved: $\widehat{q}\left(z_i\right)=\widehat{q}\left(z_{i-\mathrm{1}}\right)$.

Should the parcel become super-saturated, the same non-linear equations as used by Hydro-ABC's micro-physics scheme to adjust $\widehat{b^{\prime }}\left(z_i\right)$ and $\widehat{q}\left(z_i\right)$ towards saturation are used here. This calculates revised values of $\widehat{b^{\prime }}\left(z_i\right)$ and $\widehat{q}\left(z_i\right)$ (call these ${\widehat{b^{\prime }}}_{\mathrm{mp}}\left(z_i\right)$ and ${\widehat{q}}_{\mathrm{mp}}\left(z_i\right)$ respectively). The equations are given in point 1a of Sect. 2.3, where a hat should be added to those equations to make them relevant here. Once ${\widehat{b^{\prime }}}_{\mathrm{mp}}\left(z_i\right)$ and ${\widehat{q}}_{\mathrm{mp}}\left(z_i\right)$ have been found, $\widehat{b^{\prime }}\left(z_i\right)$ and $\widehat{q}\left(z_i\right)$ are then overwritten respectively by these adjusted values. This process is repeated by replacing $z_i\to z_{i+\mathrm{1}}$ until the top of the model is reached.

Once this procedure is complete, two special levels are found. Starting from the bottom of the domain, the first level encountered that has the property $\widehat{b^{\prime }}\ge b_{\mathrm{env}}^{\prime }$ is called the level of free convection, z_LFC, and the first level above that where $\widehat{b^{\prime }}\le b_{\mathrm{env}}^{\prime }$ is called the level of neutral buoyancy, z_LNB. CAPE is then computed from the integral:

$$\begin{array}{}\text{(C6)}& \mathrm{CAPE}=\underset{z_{\mathrm{LFC}}}{\overset{z_{\mathrm{LNB}}}{\int }}\left(\widehat{b^{\prime }}\left(z\right)-b_{\mathrm{env}}^{\prime }\left(z\right)\right)\mathrm{d}z,\end{array}$$

which has units of joules per kilogram (J kg⁻¹). The related quantity convective inhibition (CIN) can also be found from the parcel's profile:

$$\begin{array}{}\text{(C7)}& \mathrm{CIN}=-\underset{\mathrm{0}}{\overset{z_{\mathrm{LFC}}}{\int }}\left(\widehat{b^{\prime }}\left(z\right)-b_{\mathrm{env}}^{\prime }\left(z\right)\right)\mathrm{d}z.\end{array}$$

A large CAPE is associated with large positive differences, $\widehat{b^{\prime }}\left(z\right)-b_{\mathrm{env}}^{\prime }\left(z\right)$, between z_LFC and z_LNB, which are associated with buoyant air (air locally warmer than its surroundings). Similarly a large value of CIN is associated with the reverse – cooler air at a given location compared to its surroundings.

Should a value of z_LFC not be found (owing to $\widehat{b^{\prime }}<b_{\mathrm{env}}^{\prime }$ for all levels), CAPE would be zero and CIN would be found by integrating over the entire profile. Should a value of z_LNB not be found (owing to $\widehat{b^{\prime }}>b_{\mathrm{env}}^{\prime }$ for all levels above z_LFC), then CAPE would be found by integrating to the top of the model.

Example profiles of the parcel and environmental buoyancy are shown in Fig. C1 for three positions in the domain at t=10 min. In the convection zone at 148.5 km (left panel), $\widehat{b^{\prime }}>b_{\mathrm{env}}^{\prime }$ for all levels, so CAPE has a large value, consistent with the strong convection there. Two further profiles are shown for the moist zone. One, at 271.5 km (middle panel), is where CAPE remains large but where a value of z_LFC is found. There is therefore a non-zero CIN associated with this position, but the moist instability still has influence. The other, at 286.5 km (right panel), is just past the position where CAPE drops away (see Fig. 7). Here $\widehat{b^{\prime }}<b_{\mathrm{env}}^{\prime }$ for all levels, so there is no CAPE but a large value of CIN. This is the situation in the dry zone too, where a very similar plot is found, for instance at 420 km (not shown).

Figure C1 Example profiles of a parcel's computed profile ($\widehat{b^{\prime }}$, solid line) and the environmental profile ($b_{\mathrm{env}}^{\prime }$, dashed line) for three horizontal positions in the domain (calculated at t=10 min). The horizontal positions are in the convection zone (x=148.5 km) and two points within the moist zone (271.5, and 286.5 km). The area of the hatched red region is CAPE, and the area of the hatched turquoise region is CIN.