the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
The HydroABC model (Version 2.0): a simplified convectivescale model with moist dynamics
Jiangshan Zhu
Ross Noel Bannister
The prediction of convection (in terms of position, timing, and strength) is important to achieve for highresolution weather forecasting. This problem requires not only good convectivescale 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 operationalscale numerical weather prediction systems, and so a simplified model of convectivescale flow is under development (called the “ABC model”). This paper extends the existing ABC model of dry convectivescale flow to include mixing ratios of vapour and condensate phases of water. The revised model is called “HydroABC”.
HydroABC includes transport of the vapour and condensate mixing ratios within a dynamical core, and it transitions between these two phases via a microphysics 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 potentialtemperaturelike 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 domaintotal 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 inertiogravity and acoustic modes over a wide range of frequencies. This behaviour means that HydroABC is a sufficiently challenging model which will allow experimentation with innovative data assimilation strategies in future work. An ensemble of HydroABC 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 flowdependent 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 HydroABC is developed.
 Article
(7274 KB)  Fulltext XML
 BibTeX
 EndNote
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 nonlinear and inherently less predictable than largerscale 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 fastgrowing 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 kilometrescale grid lengths), some aspects of numerical models have become simpler, as much of the finescale 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 quasistatic 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).
Highresolution models of the atmosphere and their DA counterparts remain very expensive to run, requiring costly supercomputers. It is therefore desirable to study the convectivescale 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 nonrotating 1D shallowwater 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 nontrivial way. When the total height, H+h (where H is the topography height), is less than a specified level, the system is considered subsaturated 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 “lowpressure” 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 convectionlike 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 convectivescale data assimilation.
The ABC model (Petrie et al., 2017) was separately developed as a lowcost model of convectivescale behaviour also to be used as a basis to investigate convective DA methods. The ABC model is based on the threedimensional Euler equations but is reduced to two dimensions (longitude and height). It is therefore more complex than the shallowwaterbased 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 midlatitude 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 (ensemblevariational) 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.
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 “HydroABC”. The ABC model comprises five variables on a twodimensional (longitude–height) grid: zonal wind u, meridional wind v, vertical wind w, scaled density perturbation ${\stackrel{\mathrm{\u0303}}{\mathit{\rho}}}^{\prime}$ (a pressurelike variable), and buoyancy perturbation b^{′} (a potentialtemperaturelike variable). HydroABC 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). HydroABC 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 HydroABC model equations
The continuous HydroABC equations are as follows:
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, ρ_{0} is a reference density, $\stackrel{\mathrm{\u0303}}{\mathit{\rho}}$ is the scaled density ($\stackrel{\mathrm{\u0303}}{\mathit{\rho}}=({\mathit{\rho}}_{\mathrm{0}}+{\mathit{\rho}}^{\prime})/{\mathit{\rho}}_{\mathrm{0}}$), p^{′} is the pressure perturbation (diagnosed purely from ${\stackrel{\mathrm{\u0303}}{\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 twodimensional vector derivative. The key differences between HydroABC 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 HydroABC 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\phantom{\rule{0.33em}{0ex}}\mathit{\rho}$, and total dry energy, $\int \int \mathrm{d}x\mathrm{d}z\phantom{\rule{0.125em}{0ex}}{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 HydroABC'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{\u0303}}{\mathit{\rho}}{q}_{\mathrm{t}}\right)/\partial t+B\mathrm{\nabla}\cdot \left(\stackrel{\mathrm{\u0303}}{\mathit{\rho}}{q}_{\mathrm{t}}\mathit{u}\right)=\mathrm{0}$. Multiplying this equation by the constant reference density ρ_{0}, recognising that $\mathit{\rho}={\mathit{\rho}}_{\mathrm{0}}\stackrel{\mathrm{\u0303}}{\mathit{\rho}}$, and integrating over the domain give the following equation for the total amount of water in the domain:
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{\u0303}}{\mathit{\rho}}({u}^{\mathrm{2}}+{v}^{\mathrm{2}}+{w}^{\mathrm{2}})/\mathrm{2}$, the buoyancy energy is ${E}_{\mathrm{b}}=\stackrel{\mathrm{\u0303}}{\mathit{\rho}}{{b}^{\prime}}^{\mathrm{2}}/\left(\mathrm{2}{A}^{\mathrm{2}}\right)$, and the elastic energy is ${E}_{\mathrm{e}}=C{{\stackrel{\mathrm{\u0303}}{\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 HydroABC there is an additional latent heat energy which we define as
where L_{v} is the specific latent heat of vaporisation of water. In HydroABC, 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, HydroABC is designed to conserve $\int \int \mathrm{d}x\mathrm{d}z\phantom{\rule{0.125em}{0ex}}{E}_{\mathrm{moist}}$ (see Appendix A). We, however, solve HydroABC'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 microphysics step. A similar timestep 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:
where the large brackets indicate that the changes are performed over the dynamical core (dc) or the microphysics (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 nonzero, but the microphysics 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 microphysics 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, HydroABC 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 (microphysics, described below). The pointwise conservation of energy for the microphysics step, Δ^{mp}E_{moist}=0 (where ${\mathrm{\Delta}}^{\mathrm{mp}}\u2022\approx {\left(\partial \u2022/\partial t\right)}^{\mathrm{mp}}\mathrm{\Delta}t$), is the basis for our microphysics parametrisation, which is now described. The purpose of the microphysics step is to determine values of ${\mathcal{S}}_{{\mathrm{b}}^{\prime}}$ and Ev−Co. There are no general closedform expressions for these quantities, so some equations need to be solved implicitly, which are established in this section.
The relationship between the microphysical 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{\u0303}}{\mathit{\rho}}{L}_{\mathrm{v}}{\mathrm{\Delta}}^{\mathrm{mp}}q$. The above constraint that Δ^{mp}E_{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 potentialtemperaturelike variable (hence only E_{b} appears in the microphysics 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
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 microphysics 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 pseudocode presented in Appendix B.
The task now is to determine (Ev−Co)Δt, which is the change in q (or minus the change in q_{c}) over the microphysics step. The following descriptions are made with reference to steps in the pseudocode in Appendix B (in the form of ApB, step x). Consider the following change in vapour over the microphysics step, which depends on whether q is super or subsaturated after the dynamical core is run for one time step:
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 supersaturation (Fig. 1a and ApB, 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 microphysics 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{\u0303}}{\mathit{\rho}}}^{\prime}$ and b^{′}. Changing b^{′} as prescribed by Eq. (8) will therefore change q_{s}. Equation (8) therefore represents a nonlinear 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{\u0303}}{\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{\u0303}}{\mathit{\rho}}}^{\prime}$ is constant over the microphysics step. When ${b}^{\prime}\ge \mathrm{0}$, this is done using a Newton–Raphson procedure while assuming ${\stackrel{\mathrm{\u0303}}{\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{\u0303}}{\mathit{\rho}}}^{\prime},{b}_{\mathrm{mp}}^{\prime}({q}_{\mathrm{mp}}\left)\right)\right)$. See ApB, 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{\u0303}}{\mathit{\rho}}{L}_{\mathrm{v}}\leftq{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 premicrophysics 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 nonlinear equation in (a) above is not used in this case. See ApB, 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.
 a.
 2.
The second case in Eq. (7) concerns evaporation of condensate into vapour (see ApB, 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 ApB, 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 microphysics step evaporates all available q_{c}, and see ApB, step 4b), and when the third argument is smallest, evaporation is limited by available buoyancy energy (panel d, where the microphysics step converts all available E_{b}, and see ApB, 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 nonlinear equations mentioned previously (but overall nonlinearity 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 supersaturation but when the threshold for condensation is not met or for subsaturation but when ${b}^{\prime}<\mathrm{0}$). In this case there are no microphysicsrelated changes to the fields.
The change in buoyancy due to microphysics, ${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:
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, ρ_{0}(z), is required. In the ABC model (Petrie et al., 2017), ρ_{0} 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}}\mathrm{exp}(z/H)$, where ρ_{00} is the surface air density (1.225 kg m^{−3}) and H is the scale height (9 km). For pressure, $p={p}_{\mathrm{0}}+{p}^{\prime}$, where p_{0} (the reference pressure) is in hydrostatic balance with ρ_{0}, $\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 ρ_{0}(z), is ${p}_{\mathrm{0}}\left(z\right)={p}_{\mathrm{00}}H{\mathit{\rho}}_{\mathrm{00}}g\left(\mathrm{1}\mathrm{exp}(z/H)\right)$, where p_{00} is the reference surface pressure. The relationship between surface reference pressure and surface reference density is p_{00}=Hρ_{00}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{\u0303}}{\mathit{\rho}}}^{\prime}$ to p^{′} and, using the above form of p_{o}(z), gives
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).^{1} 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, θ_{0}(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 θ_{00} is the surface potential temperature (300 K). The potential temperature is therefore
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
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{\u0303}}{\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 HydroABC, q and q_{c}, are each stored at the same location as $\stackrel{\mathrm{\u0303}}{\mathit{\rho}}$ and ${\stackrel{\mathrm{\u0303}}{\mathit{\rho}}}^{\prime}$. This choice makes sense from a numerical perspective as both q and q_{c} are multiplied by $\stackrel{\mathrm{\u0303}}{\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 HydroABC, 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{\u0303}}{\mathit{\rho}}\partial q/\partial t+q\partial \stackrel{\mathrm{\u0303}}{\mathit{\rho}}/\partial t+Bq\mathrm{\nabla}\cdot \left(\stackrel{\mathrm{\u0303}}{\mathit{\rho}}\mathit{u}\right)+B\stackrel{\mathrm{\u0303}}{\mathit{\rho}}\mathit{u}\cdot \mathrm{\nabla}q=\stackrel{\mathrm{\u0303}}{\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:
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 subtime 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 subtime 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 microphysics parametrisation of Sect. 2.3 using values of variables after the dynamical core is called.
In this section we show example behaviour of HydroABC. The initial conditions for the dry variables (u, v, w, ${\stackrel{\mathrm{\u0303}}{\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{\u0303}}{\mathit{\rho}}}^{\prime}$ is set using geostrophic balance, b^{′} is set using hydrostatic balance, and w is set by imposing a 3D divergencefree 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^{−1} 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^{−1}, 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.
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^{−1} (not shown), but again once B moderation is taken into account; this is equivalent to a horizontal mass transport of ±5 m s^{−1}. 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
HydroABC 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\phantom{\rule{0.125em}{0ex}}{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 fastmoving 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.
3.3 Excitation of waves by convection
The ABC model is capable of supporting Rossbylike, inertiogravity, 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 HydroABC, 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 HydroABC (right). The top panels are Hovmöller plots of w for the longitude–time domain. Both versions of the model show westward and eastwardpropagating 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 biFourier 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 lowerfrequency inertiogravity 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}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{\mathrm{1}}$. This is the wave frequency of the highestfrequency inertiogravity 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 higherfrequency 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 inertiogravity 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 highwavenumber acoustic waves have frequencies lower than the Lamb line (visible especially in the topright section of the bottomright panel). This is a numerical effect where the frequencies of highwavenumber waves are too low, leading to numerical dispersion (see e.g. Chap. 2 of Durran, 1999).
Comparing the dry ABC and HydroABC 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 HydroABC causes are (i) the excitation of stronger inertiogravity and acoustic modes across the spectrum and (ii) broadening of their frequencies. Effect (i) is visible at all wavenumbers and frequencies, especially for inertiogravity 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 topright panel of Fig. 4, there is a gradual onset of the strength of convection, which leads to a convolution of the spectra with this “lowfrequency” onset.
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 flowdependent (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 convectivescale 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 moistureaffected covariances and correlations between different variables and positions, as well as in different regimes in HydroABC, we generate an ensemble of HydroABC forecasts. Longitude–height slices from a 93member ensemble of convectivescale UM forecasts (Bannister et al., 2017) are each processed in the same way as described in Sect. 3 to give a 93member ensemble of HydroABC 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{\u0303}}{\mathit{\rho}}}^{\prime}$ (see Eq. 10, Sect. 2.4, where p and T depend on ${\stackrel{\mathrm{\u0303}}{\mathit{\rho}}}^{\prime}$ and b^{′} as shown in Eqs. 11 and 13). The example HydroABC 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 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 nonzero 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 nonzero 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 WRFbased study of Michel et al. (2011) and in the AROMEbased 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 counterintuitive. 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 nonzero, 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” setup in the initial conditions (covariances not shown, but correlations are shown below).
Figure 6 shows vertical correlations between a different selection of quantities to those in Fig. 5. The top three rows show autocorrelations of ${\stackrel{\mathrm{\u0303}}{\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{\u0303}}{\mathit{\rho}}}^{\prime}$ appears to comprise two layers of correlations (one below 6 km and the other above), where in each, ${\stackrel{\mathrm{\u0303}}{\mathit{\rho}}}^{\prime}$ is strongly correlated vertically but only weakly (and negatively) across each layer. In the convection zone there are some smallscale 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 sublayers, which are associated with the establishment of upper and lowerpressure 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 sublayering 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 largescale (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{\u0303}}{\mathit{\rho}}}^{\prime}$ at different heights (x axis). These are included in combination with the ${\stackrel{\mathrm{\u0303}}{\mathit{\rho}}}^{\prime}$ autocorrelations 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{\u0303}}{\mathit{\rho}}}^{\prime}/\partial z={b}^{\prime}$. This relates b^{′} to vertical gradients in ${\stackrel{\mathrm{\u0303}}{\mathit{\rho}}}^{\prime}$. Thus the property of HB is assumed when the vertical derivative of the ${\stackrel{\mathrm{\u0303}}{\mathit{\rho}}}^{\prime}\text{\u2013}{\stackrel{\mathrm{\u0303}}{\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{\u2013}{\stackrel{\mathrm{\u0303}}{\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{\u0303}}{\mathit{\rho}}}^{\prime}\text{\u2013}{\stackrel{\mathrm{\u0303}}{\mathit{\rho}}}^{\prime}\text{correlation}\right)/$ $\partial z\left({\stackrel{\mathrm{\u0303}}{\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{\u0303}}{\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 nearzero values of b^{′}–${\stackrel{\mathrm{\u0303}}{\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 nearzero values of b^{′}–${\stackrel{\mathrm{\u0303}}{\mathit{\rho}}}^{\prime}$ correlations (e and f), suggesting a clear breakdown of HB. In the first half of the forecast period, there are more nonhydrostatic 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 nonhydrostatic 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 nearzero values of b^{′}–${\stackrel{\mathrm{\u0303}}{\mathit{\rho}}}^{\prime}$ correlations (g and h), suggesting a propagated breakdown of HB. Bannister et al. (2011) and VetraCarvalho et al. (2012) similarly found evidence of nonhydrostatic patterns in the covariance patterns in rainy regions (compared to dry regions) of forecast errors in the convectivescale 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 nonconvecting regions, convectivescale 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 ensemblebased systems (e.g. ensemble Kalman filters or ensemblevariational methods; Bannister, 2017) and partially in hybrid systems (where the static background error covariances of traditional variational methods are averaged with ensemblederived covariances; Bannister, 2017). It is well known though that the use of ensemblebased covariances leads to rank deficiency and undersampling (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 nonconvecting 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 HydroABC 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}}^{\mathrm{max}}$, are important signatures of active convection but not necessarily of preconvective conditions. (vi) Hydrostatic imbalance (HI) at any horizontal location is defined as the largest absolute value of $\left(C\partial {\stackrel{\mathrm{\u0303}}{\mathit{\rho}}}^{\prime}/\partial z{b}^{\prime}\right)/$ $\left(\mathrm{rms}\left(C\partial {\stackrel{\mathrm{\u0303}}{\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 nonzero 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}}^{\mathrm{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 finescale but does not peak at the correct location of the future convection. The quantities w^{max}, ${q}_{\mathrm{c}}^{\mathrm{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 nonflowdependent, as in traditional variational schemes. Further work will be required to test the reliability of these indicators.
Described in this paper is how the simplified ABC model of convectivescale flow has been extended to include moist processes, namely evaporation and condensation, yielding a revised model called HydroABC. This is the next step in the aim to provide a lowcost framework to study convectivescale 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 nonlinear system of equations as set out in Sect. 2.3. The latent heat of vaporisation, L_{v} (set to 2500 J g^{−1}), is introduced and relates changes in vapour to changes in latent heat energy, Eq. (3).

Condensation happens when situations are supersaturated (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 nonlinear 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 subsaturated, 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{\u0303}}{\mathit{\rho}}}^{\prime}$ as set out in Sect. 2.4. The conditions used to relate b^{′} and ${\stackrel{\mathrm{\u0303}}{\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.
HydroABC 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 HydroABC; this is to keep the system as simple as possible while allowing highly nonlinear processes.
An example run of HydroABC shows a strong convective event driven by latent heat release. This forms a realistic anvil cloud (Sect. 3.1) and the excitation of inertiogravity 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 convectivescale 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 flowdependent statistics (found from a 93member ensemble) can depend on the presence of convection. It is found that there can be significant differences between the covariances in convecting and nonconvecting 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 waterrelated 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.
The formulation of HydroABC is based on a timestep split into a dynamical core and a microphysics parametrisation. HydroABC can also be represented as continuous equations, i.e. as Eqs. (1a–h) but with Eq. (1e) replaced with
which is an explicit form of the buoyancy equation. This is consistent with the basis of our microphysics 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
i.e. the righthand side of Eq. (A1). In order to show that this leads to an energyconserving system of equations, multiply Eq. (A1) by $\stackrel{\mathrm{\u0303}}{\mathit{\rho}}{b}^{\prime}/{A}^{\mathrm{2}}$:
Then use the following identity (Petrie et al., 2017, their Eq. 20):
with $\mathit{\gamma}={{b}^{\prime}}^{\mathrm{2}}/\left(\mathrm{2}{A}^{\mathrm{2}}\right)$. That is,
Making the substitution ${E}_{\mathrm{b}}=\stackrel{\mathrm{\u0303}}{\mathit{\rho}}{{b}^{\prime}}^{\mathrm{2}}/\left(\mathrm{2}{A}^{\mathrm{2}}\right)$ (see Sect. 2.2.2), gives
This is the continuous form of the equation for buoyancy energy. Next multiply Eq. (1g) by L_{v},
and make the substitution from Eq. (3),
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):
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}}$:
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:
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 righthand side. This issue is avoided by using the microphysics parametrisation step represented by Eq. (6).
Annotated pseudocode describing the dynamics/microphysics parametrisation is given in this part of the Appendix.
 1.
Propagate the HydroABC fields u, v, w, ${\stackrel{\mathrm{\u0303}}{\mathit{\rho}}}^{\prime}$, b^{′}, q, and q_{c} (described by Eq. 1) by one time step.
 2.
Calculate the premicrophysics 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 microphysics 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 nonlinear 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)}\phantom{\rule{0.25em}{0ex}}\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{\u0303}}{\mathit{\rho}}{L}_{\mathrm{v}}\leftq{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).
 a.
 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 semidefinite 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{\u0303}}{\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}$$
 a.
 5.
Overwrite the model's variables, b^{′}, q, and q_{c}, with their postmicrophysics values, ${b}_{\mathrm{mp}}^{\prime}$, q_{mp}, and ${{q}_{\mathrm{c}}}_{\mathrm{mc}}$ respectively.
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}}\phantom{\rule{0.125em}{0ex}}\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 HydroABC there is no explicit temperature variable, but we can derive a physically reasonable CAPElike 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
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
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):
Taking the difference between Eqs. (C2) and (C3) gives
where in the second line we have assumed mechanical equilibrium (that $\widehat{\stackrel{\mathrm{\u0303}}{{\mathit{\rho}}^{\prime}}}={\stackrel{\mathrm{\u0303}}{\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^{2} 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}:
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 supersaturated, the same nonlinear equations as used by HydroABC's microphysics 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:
which has units of joules per kilogram (J kg^{−1}). The related quantity convective inhibition (CIN) can also be found from the parcel's profile:
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 nonzero 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).
The HydroABC V2.0 code is written in Fortran 90 and is freely available via a GitHub repository (https://github.com/rossbannister/HydroABC_vn2.0, last access: 19 October 2023) and Zenodo (https://doi.org/10.5281/zenodo.7418510, Bannister, 2022). The HydroABC V2.0 user documentation is available in the same place. Example HydroABC initial condition data (as used in this paper) are available at https://doi.org/10.5281/zenodo.10025251 (Bannister and Zhu, 2023).
The contact author has declared that neither of the authors has any competing interests.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.
The authors would like to thank the National Natural Science Foundation of China for funding Jiangshan Zhu (grant ref 41875136) and the Natural Environment Research Council, which provides national capability funding for Ross Noel Bannister via the National Centre for Earth Observation (contract number PR140015). The authors would also like to thank Mike Cullen for correspondence during the development of HydroABC and the two anonymous reviewers for providing valuable and insightful feedback with their reviews.
This research has been supported by the National Natural Science Foundation of China (grant no. 41875136) and the National Centre for Earth Observation (contract number PR140015).
This paper was edited by Chiel van Heerwaarden and reviewed by two anonymous referees.
Bannister, R.: A review of operational methods of variational and ensemblevariational data assimilation, Q. J. Roy. Meteor. Soc., 143, 607–633, https://doi.org/10.1002/qj.2982, 2017. a, b
Bannister, R.: rossbannister/HydroABC_vn2.0: HydroABC_vn_2.0 (HydroABC_vn_2.0), Zenodo [code], https://doi.org/10.5281/zenodo.7418510, 2022. a
Bannister, R. and Zhu, J.: Ensemble of initial conditions for HydroABC, Zenodo [data set], https://doi.org/10.5281/zenodo.10025251, 2023. a
Bannister, R. N.: A review of forecast error covariance statistics in atmospheric variational data assimilation. II: Modelling the forecast error covariance statistics, Q. J. Roy. Meteor. Soc., 134, 1971–1996, 2008. a
Bannister, R. N.: The ABCDA system (v1.4): a variational data assimilation system for convectivescale assimilation research with a study of the impact of a balance constraint, Geosci. Model Dev., 13, 3789–3816, https://doi.org/10.5194/gmd1337892020, 2020. a
Bannister, R. N.: Balance conditions in variational data assimilation for a highresolution forecast model, Q. J. Roy. Meteor. Soc., 147, 2917–2934, https://doi.org/10.1002/qj.4106, 2021. a
Bannister, R. N., Migliorini, S., and Dixon, M.: Ensemble prediction for nowcasting with a convectionpermitting model–II: Forecast error statistics, Tellus A, 63, 497–512, 2011. a
Bannister, R. N., Migliorini, S., Rudd, A. C., and Baker, L. H.: Methods of investigating forecast error sensitivity to ensemble size in a limitedarea convectionpermitting ensemble, Geosci. Model Dev. Discuss. [preprint], https://doi.org/10.5194/gmd2017260, in review, 2017. a
Bannister, R. N., Chipilski, H., and MartinezAlvarado, O.: Techniques and challenges in the assimilation of atmospheric water observations for numerical weather prediction towards convective scales, Q. J. Roy. Meteor. Soc., 146, 1–48, https://doi.org/10.1002/qj.3652, 2020. a, b, c
Banos, I. H., Mayfield, W. D., Ge, G., Sapucci, L. F., Carley, J. R., and Nance, L.: Assessment of the data assimilation framework for the Rapid Refresh Forecast System v0.1 and impacts on forecasts of a convective storm case study, Geosci. Model Dev., 15, 6891–6917, https://doi.org/10.5194/gmd1568912022, 2022. a
Berre, L.: Estimation of synoptic and mesoscale forecast error covariances in a limitedarea model, Mon. Weather Rev., 128, 644–667, 2000. a
Clark, P., Roberts, N., Lean, H., Ballard, S. P., and CharltonPerez, C.: Convectionpermitting models: a stepchange in rainfall forecasting, Meteorol. Appl., 23, 165–181, 2016. a
Cullen, M. and Davies, T.: A conservative splitexplicit integration scheme with fourthorder horizontal advection, Q. J. Roy. Meteor. Soc., 117, 993–1002, 1991. a
Durran, D. R.: Numerical methods for wave equations in geophysical fluid dynamics, Springer, New York, ISBN 0387983767. 1999. a
Ehrendorfer, M.: A review of issues in ensemblebased Kalman filtering, Meteorol. Z., 16, 795–818, 2007. a
Errico, R. M., Bauer, P., and Mahfouf, J.F.: Issues regarding the assimilation of cloud and precipitation data, J. Atmos. Sci., 64, 3785–3798, 2007. a
Fabry, F. and Meunier, V.: Why are radar data so difficult to assimilate skillfully?, Mon. Weather Rev., 148, 2819–2836, 2020. a
Gill, A. E.: AtmosphereOcean Dynamics, Academic Press, ISBN 0122835220, 1982. a
Gustafsson, N., Janjić, T., Schraff, C., Leuenberger, D., Weissmann, M., Reich, H., Brousseau, P., Montmerle, T., Wattrelot, E., Bucanek, A., Mile, M., Hamdi, R., Lindskog, M., Barkmeijer, J., Dahlbom, M. Macpherson, B., Ballard, S., Inverarity, G., Carley, J., Alexander, C., Dowell, D., Liu, S., Ikuta, Y., and Fujita, T.: Survey of data assimilation methods for convectivescale numerical weather prediction at operational centres, Q. J. Roy. Meteor. Soc., 144, 1218–1256, 2018. a
Hamill, T. M., Whitaker, J. S., and Snyder, C.: Distancedependent filtering of background error covariance estimates in an ensemble Kalman filter, Mon. Weather Rev., 129, 2776–2790, 2001. a
Hohenegger, C. and Schär, C.: Atmospheric predictability at synoptic versus cloudresolving scales, B. Am. Meteorol. Soc., 88, 1783, https://doi.org/10.1175/BAMS88111783, 2007. a
Holton, J. and Hakim, G.: An Introduction to Dynamic Meteorology, 5th Edn., Academic Press, Waltham, MA, ISBN 9780123848666, 2013. a, b
Houtekamer, P. and Zhang, F.: Review of the Ensemble Kalman Filter for Atmospheric Data Assimilation, Mon. Weather Rev., 144, 4489–4532, 2016. a
Houtekamer, P. L. and Mitchell, H. L.: Data assimilation using an ensemble Kalman filter technique, Mon. Weather Rev., 126, 796–811, 1998. a
Houtekamer, P. L. and Mitchell, H. L.: Ensemble Kalman filtering, Q. J. Roy. Meteor. Soc., 131, 3269–3289, 2005. a
Kent, T., Bokhove, O., and Tobias, S.: Dynamics of an idealized fluid model for investigating convectivescale data assimilation, Tellus A, 69, 1369332, https://doi.org/10.1080/16000870.2017.1369332, 2017. a
Lee, J. C. K., Amezcua, J., and Bannister, R. N.: Hybrid ensemblevariational data assimilation in ABCDA within a tropical framework, Geosci. Model Dev., 15, 6197–6219, https://doi.org/10.5194/gmd1561972022, 2022. a
Leung, T. Y., Leutbecher, M., Reich, S., and Shepherd, T. G.: Atmospheric predictability: revisiting the inherent finitetime barrier, J. Atmos. Sci., 76, 3883–3892, https://doi.org/10.1175/JASD190057.1, 2019. a
Lorenc, A.: A study of ob monitoring statistics from radiosondes, composited for lowlevel cloud layers, Met Office NWP Forecasting Research Technical Report, 504, 1–32, 2007. a
Lorenz, E. N.: The predictability of a flow which possesses many scales of motion, Tellus, 21, 289–307, 1969. a
Ménétrier, B. and Montmerle, T.: Heterogeneous backgrounderror covariances for the analysis and forecast of fog events, Q. J. Roy. Meteor. Soc., 137, 2004–2013, 2011. a
Michel, Y., Auligné, T., and Montmerle, T.: Heterogeneous convectivescale background error covariances with the inclusion of hydrometeor variables, Mon. Weather Rev., 139, 2994–3015, 2011. a, b
Montmerle, T.: Optimization of the assimilation of radar data at the convective scale using specific background error covariances in precipitation, Mon. Weather Rev., 140, 3495–3506, 2012. a
Montmerle, T. and Berre, L.: Diagnosis and formulation of heterogeneous backgrounderror covariances at the mesoscale, Q. J. Roy. Meteor. Soc., 136, 1408–1420, 2010. a, b, c, d
Petrie, R. E., Bannister, R. N., and Cullen, M. J. P.: The “ABC model”: a nonhydrostatic toy model for use in convectivescale data assimilation investigations, Geosci. Model Dev., 10, 4419–4441, https://doi.org/10.5194/gmd1044192017, 2017. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r
Pielke, R.: Mesoscale Meteorological Modeling, Academic Press, San Diego, California, ISBN 9780123852380, 2002. a
Salby, M. L.: Fundamentals of atmospheric physics, vol. 61, Academic press, San Diego, California, ISBN 9786611049874, 1996. a, b, c
Skamarock, W. C., Klemp, J. B., Dudhia, J., Gill, D. O., Liu, Z., Berner, J., Wang, W., Powers, J. G., Duda, M. G., Barker, D. M., and Huang, X.Y.: A Description of the Advanced Research WRF Model Version 4.3, NCAR/TN556+STR, National Center for Atmospheric Research, https://doi.org/10.5065/1dfh6p97, 2021. a
Sun, J., Xue, M., Wilson, J. W., Zawadzki, I., Ballard, S. P., OnvleeHooimeyer, J., Joe, P., Barker, D. M., Li, P.W., Golding, B., Xu, M., and Pinto, J.: Use of NWP for nowcasting convective precipitation: Recent progress and challenges, B. Am. Meteorol. Soc., 95, 409–426, 2014. a
van Leeuwen, P. J.: Comment on “Data assimilation using an ensemble Kalman filter technique”, Mon. Weather Rev., 127, 1374–1377, 1999. a
VetraCarvalho, S., Dixon, M., Migliorini, S., Nichols, N. K., and Ballard, S. P.: Breakdown of hydrostatic balance at convective scales in the forecast errors in the Met Office Unified Model, Q. J. Roy. Meteor. Soc., 138, 1709–1720, 2012. a
Würsch, M. and Craig, G. C.: A simple dynamical model of cumulus convection for data assimilation research, Meteorol. Z., 23, 483–490, 2014. a, b
Yang, Y., Gao, S., Wang, Y., and Shi, H.: Impact of FeatureDependent Static Background Error Covariances for SatelliteDerived Humidity Assimilation on Analyses and Forecasts of Multiple Sea Fog Cases over the Yellow Sea, Remote Sens., 14, 4537, https://doi.org/10.3390/rs14184537, 2022. a
Yano, K.I., Ziemiański, M. Z., Cullen, M., Termonia, P., Onvlee, J., Bengtsson, L., Carrassi, A., Davy, R., Deluca, A., Gray, S. L., Homar, V., Köhler, M., Krichak, S., Michaelides, S., Phillips, V. T. J., Soares, P. M. M., and Wyszogrodzki, A. A.: Scientific challenges of convectivescale numerical weather prediction, B. Am. Meteorol. Soc., 99, 699–710, 2018. a
Zhang, M. and Zhang, F.: E4DVar: Coupling an ensemble Kalman filter with fourdimensional variational data assimilation in a limitedarea weather prediction model, Mon. Weather Rev., 140, 587–600, 2012. a
This relationship between b^{′} and θ^{′} does not imply that b^{′} has the full properties of potential temperature in the ABC model (e.g. that b^{′} is conserved during adiabatic processes). In Petrie et al. (2017) this relationship was used as a starting point in the derivation of the ABC equations from the Euler equations, and it is reused here as a pragmatic way of deriving temperature from the ABC variables.
 Abstract
 Introduction
 Dry ABC to HydroABC
 Example run of the model
 Convectiondependent covariances and identification of areas prone to convection
 Conclusions
 Appendix A: Continuous total energy equations of HydroABC
 Appendix B: Pseudocode for the microphysics parametrisation
 Appendix C: Computation of the convective available potential energy (CAPE)
 Code and data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
 Abstract
 Introduction
 Dry ABC to HydroABC
 Example run of the model
 Convectiondependent covariances and identification of areas prone to convection
 Conclusions
 Appendix A: Continuous total energy equations of HydroABC
 Appendix B: Pseudocode for the microphysics parametrisation
 Appendix C: Computation of the convective available potential energy (CAPE)
 Code and data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References