Articles | Volume 16, issue 21
Model description paper
31 Oct 2023
Model description paper |  | 31 Oct 2023

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

Jiangshan Zhu and Ross Noel Bannister

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är2007; Leung et al.2019; Lorenz1969), 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 Berre2010).

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 Meunier2020).

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 ϕ=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 ϕ=gH+ϕ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 (Bannister2020). 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 (Bannister2021). 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 ρ̃ (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 qc 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+qc, 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:


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, u=(u,v,w), f is the Coriolis parameter, ρ0 is a reference density, ρ̃ is the scaled density (ρ̃=(ρ0+ρ)/ρ0), p is the pressure perturbation (diagnosed purely from ρ̃), Sb 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 =/x/z is the two-dimensional vector derivative. The key differences between Hydro-ABC and the dry ABC system are Sb in Eq. (1e) and the additional equations, Eqs. (1g) and (1h). How Sb 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, dxdzρ, and total dry energy, dxdzEdry, where the integration is performed over the whole domain (Edry 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, qt, is the sum of the vapour and condensate mixing ratios, qt=q+qc. Adding Eqs. (1g) and (1h) gives the governing equation for qt: (ρ̃qt)/t+Bρ̃qtu=0. Multiplying this equation by the constant reference density ρ0, recognising that ρ=ρ0ρ̃, and integrating over the domain give the following equation for the total amount of water in the domain:

(2) t ρ q t d x d z + B ρ q t u d x d z = 0 .

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 Ek=ρ̃(u2+v2+w2)/2, the buoyancy energy is Eb=ρ̃b2/(2A2), and the elastic energy is Ee=Cρ̃2/(2B). These make up the total dry energy, Edry=Ek+Eb+Ee. In Hydro-ABC there is an additional latent heat energy which we define as

(3) E l = ρ ̃ L v q ,

where Lv 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, Emoist=Edry+El=Ek+Eb+Ee+El.

In its continuous form, Hydro-ABC is designed to conserve dxdzEmoist (see Appendix A). We, however, solve Hydro-ABC's Eq. (1) over a time step Δt in two parts: a “dynamical core” in which Sb, 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 qc). The evolution of the total energy is represented by the following sum:

(4) E moist t = E dry t dc + E l t dc E moist / t dc + E l t mp + E b t mp E moist / t mp ,

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 El) 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 Emoist/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 Sb, Ev, and Co

As explained above, Hydro-ABC solves Eq. (1) over a time step Δt in two stages, by first assuming Sb 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, ΔmpEmoist=0 (where Δmp/tmpΔ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 Sb 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: Δmpq=Ev-CoΔt. Further, the relationship between the change in El and the change in q is ΔmpEl=ρ̃LvΔmpq. The above constraint that ΔmpEmoist=0 may be put in the following way: ΔmpEmoist=ΔmpEk+ΔmpEb+ΔmpEe+ΔmpEl=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 Eb appears in the micro-physics terms of Eq. 4). The last equation then simplifies to ΔmpEb+ΔmpEl=0. Combining the above formulae with the form of Eb (which is given in the opening paragraph of Sect. 2.2.2) yields


Note that Δmpb=SbΔt. Equation (6) will now be used with further rules (below) to determine how our micro-physics scheme will modify q, qc, and b (which will reveal how Sb and Ev−Co are found). The description of the scheme may be read in combination with the pseudo-code presented in Appendix B.

Figure 1Graphical representation of Eq. (7), showing the change in q, qc, 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 qc) 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:

(7)Ev-CoΔt=-q-qscondensation(super-saturated):when qqs andb>0 or b<0 and ρ̃Lvq-qs/Eb>Γ,minqs-q1-e-Δt/τ,qc,b2/2A2Lvevaporation(sub-saturated):when q<qs and b>0,0otherwise,

where qs 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 qqs (red line) and qc will increase by the same amount (blue line). According to Eq. (6), associated with this step is an increase in b2 by the amount Δmpb2=2A2Lv(q-qs). This is due to the release of latent heat, which is converted to buoyancy energy. Writing Δmpb2 as bmp2-b2 (where bmp is the value of buoyancy after micro-physics and b is the value before), the equation for Δmpb2 translates to

    (8) b mp = + b 2 + 2 A 2 L v ( q - q s )

    (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 qs is shown in Sect. 2.4 to be a function of ρ̃ and b. Changing b as prescribed by Eq. (8) will therefore change qs. Equation (8) therefore represents a non-linear equation. This is solved by finding the root of the function f(qmp)=qmp-qs(ρ̃,bmp(qmp)) given that bmp is a function of qmp, i.e. bmp=+b2+2A2Lv(q-qmp), and assuming that ρ̃ is constant over the micro-physics step. When b0, this is done using a Newton–Raphson procedure while assuming ρ̃ is constant. The value of (Ev−Co)Δt for the first branch of Eq. (7) is then -q-qs(ρ̃,bmp(qmp)). See Ap-B, step 3a.

    • b.

      Step a, however, becomes potentially unphysical if b is large and negative and qqs 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 ρ̃Lvq-qs/Eb>Γ, where Γ≫1. This is the condition (which must be met only when b<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<0, the amount of latent heat released is based on the value of qs at the pre-micro-physics value of b and not on the value of qs 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 qs with timescale τ (Fig. 1b, red line) and qc relaxes downwards (blue line). According to Eq. (6), associated with this step is a change in b2 by the amount Δmpb2=-2A2LvEv-CoΔt (where Ev-Co>0). This translates to a new value bmp according to the following:

    (9) b mp = + b 2 - 2 A 2 L v Ev - Co Δ t ,

    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>0 because, in the case of b<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<0, Eq. (9) will yield an increased b (i.e. bmp-b>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 qc (panel c, where the micro-physics step evaporates all available qc, 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 Eb, 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 qs before evaporation lowers the value of b (and hence lowers qs). 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<0). In this case there are no micro-physics-related changes to the fields.

The change in buoyancy due to micro-physics, bmp-b, yields the value SbΔt in Eq. (1e), and the change in vapour, qmpq, gives the value (Ev−Co)Δt in Eqs. (1g) and (1h).

2.4 The saturated vapour mixing ratio

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

(10) q s = 1000 g kg - 1 380 Pa p exp 17.3 ( T - 273.2 K ) T - 35.9 K ,

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 qs) 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 qs we let this depend on z according to ρ0(z)=ρ00exp(-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=p0+p, where p0 (the reference pressure) is in hydrostatic balance with ρ0, dp0(z)/dz=-ρ0(z)g, where g is the acceleration due to gravity. The solution to this equation, given ρ0(z), is p0(z)=p00-Hρ00g1-exp(-z/H), where p00 is the reference surface pressure. The relationship between surface reference pressure and surface reference density is p00=Hρ00g, which is derived from the previous equation with the condition p0(z)0. ABC's equation of state, Eq. (1f), relates the ABC variable ρ̃ to p and, using the above form of po(z), gives

(11) p ( x , z , t ) = p 0 ( z ) + p ( x , z , t ) = p 00 - H ρ 00 g 1 - exp ( - z / H ) + C ρ 00 exp ( - z / H ) ρ ̃ ( x , z , t ) .

In Petrie et al. (2017), b is related to a potential temperature increment θ via θ=(θR/g)b, 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 θ(x,z,t)=θR+θ0(z)+θ(x,z,t). To derive a reference profile, θ0(z), we use the definition of the Brunt–Väisälä frequency N2=(g/θR)dθ0(z)/dz (see e.g. Holton and Hakim2013). The parameter A in Eq. (1e) takes the place of N, and so we solve the equation A2=(g/θR)dθ0(z)/dz. The solution to this equation is θ0(z)=θ00-θR+(A2θR/g)z, where θ00 is the surface potential temperature (300 K). The potential temperature is therefore

(12) θ ( x , z , t ) = θ R + θ 0 ( z ) + θ ( x , z , t ) = θ 00 + θ R g A 2 z + b ( x , z , t ) .

The standard relationship for T in terms of θ and p is T(x,z,t)=θ(x,z,t)p(x,z,t)/p00κ (Holton and Hakim2013), where κ is the constant R/cp (where R is the gas constant for dry air and cp is the specific heat capacity at constant pressure; the value of κ is 0.286). Using Eqs. (12) and (11) gives

(13) T ( x , z , t ) = θ 00 + θ R g A 2 z + b ( x , z , t ) p 00 - H ρ 00 g 1 - exp ( - z / H ) + C ρ 00 exp ( - z / H ) ρ ̃ ( x , z , t ) p 00 κ .

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 ρ̃.

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 qc, are each stored at the same location as ρ̃ and ρ̃. This choice makes sense from a numerical perspective as both q and qc are multiplied by ρ̃ 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 ρ̃q/t+qρ̃/t+Bqρ̃u+Bρ̃uq=ρ̃(Ev-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:

(14) q t + B u q = Ev - Co .

A similar equation is also found for qc.

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 qc (with Ev-Co=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, ρ̃, 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.5km 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 ρ̃ 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 qc=0.

3.1 Reference run

The rest of Fig. 2 (rows 2, 3, and 4) shows the evolution of w (left column) and q, qc (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=60min that localised upward motion starts to develop at ground level, most strongly at x=148.5km but also at 156km (also not shown). By 120min, 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 ∼200m 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 10m 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.

Figure 2Initial conditions of the example reference run of Hydro-ABC (A=0.01s-1, B=0.05, C=2×104m2s-2, f=10-4s-1, Δt=0.1 s, τ=1000 s, Lv=2500Jg-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 qc=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 360min. w is shown in the left column (m s−1), and q (shading) and qc (contours) are shown in the right column. The contours are the black (1 g kg−1), orange (6 g kg−1), and red (12 g kg−1) 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 240min, an anvil cloud of vapour and condensate has established itself between heights of ∼5.5 and 8.5km, which spreads out laterally in both directions. The zonal winds at these levels peak at around ±100m s−1 (not shown), but again once B moderation is taken into account; this is equivalent to a horizontal mass transport of ±5m 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 13km have developed on either side of the main plume. Between 240 and 360min, the smaller plumes are themselves impeded by downwelling air on each side of the main plume. By 360min 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, dxdzEmoist(t)=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=60min. 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 3Evolution 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 dxdzEmoist(t)/dxdzEmoist(0), where Emoist=Ek+Eb+Ee+El (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 (Salby1996, 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/(2π)=1.59×10-3s-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 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 Durran1999).

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 4Vertical wind in the longitude–time (a, b) and wavenumber–frequency (c, d) domains for level 30 in the model (7.7km height). The left panels are for dry ABC, and the right panels are for Hydro-ABC. The initial conditions of u, v, w, ρ̃, and b are identical for ABC and Hydro-ABC, and the additional initial conditions for q and qc are the same as those used in Fig. 2. In (c) and (d), the horizontal dashed lines are each at frequency A/(2π) (the pure gravity wave frequency), and the sloped dashed lines each have gradient df/dn=L-1BC (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×103 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. Lorenc2007; Montmerle and Berre2010; 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 360min). Each ensemble member is prepared with the same q and qc fields (as in Fig. 2), but the relative humidity (q/qs, the field that matters for convection) is different between members because qs is a function of b and ρ̃ (see Eq. 10, Sect. 2.4, where p and T depend on ρ̃ 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.5km 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 5Evolution of a selection of normalised vertical covariances in the convection zone (x=148.5km) at times between 40 and 360min. The rows are for ww, uw, qw, and qcw covariances. In each plot of p1p2 covariances, the y axis is for p1 and the x axis is for p2, 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.5km). Near the start of the integrations (t=40min), all covariances shown are too small to be detected on the scale used. At 80min, convection is just starting, which appears in a region of non-zero covariances within 4.5km of the surface. There are positive covariances between w and itself at different levels, and with q (and more weakly with qc), 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.5km layer. By 120min the non-zero covariance region has grown to within 910km 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, qc 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 ∼9km. This suggests that the region of positive qcw 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 qc 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 qc. 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 qc could have a positive correlation with the vertical wind. In fact statistics show that often q and qc are mutually positively correlated at points inside the cloud (not shown).

At 160min a shallow surface layer of negative qcw covariances has developed, which may be due to a circulation being induced by the convection, which is drawing in air with lower qc from the surroundings. There is an additional region of weak, but non-zero, ww and uw 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 240min 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.5km, including at locations within and beyond the moist “bubble” set-up in the initial conditions (covariances not shown, but correlations are shown below).

Figure 6Evolution of a selection of vertical correlations at times between 40 and 360min. The top set of three rows is for ρ̃ρ̃ correlations, and the bottom set is for bρ̃ correlations. For each set of three rows, the rows are in the convection zone (x=148.5km – the same position as used for Fig. 5), the moist zone (286.5km), and the dry zone (421.5km) respectively. In each plot of p1p2 correlations, the y axis is for p1, the x axis is for p2, 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 ρ̃ρ̃ 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ρ̃ correlation. Circles are associated with small bρ̃ 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 ρ̃ 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=40min the different zones show similar patterns. The model variable ρ̃ appears to comprise two layers of correlations (one below 6km and the other above), where in each, ρ̃ 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 80min 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 uw 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 160min in the convection zone the sub-layering has disappeared (probably because the shear zone has moved upwards), and the vertical length scales above 6km 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 200min 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 ρ̃ at different heights (x axis). These are included in combination with the ρ̃ 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ρ̃/z=b. This relates b to vertical gradients in ρ̃. Thus the property of HB is assumed when the vertical derivative of the ρ̃ρ̃ correlation function at a point in the top three rows of Fig. 6 is zero and is matched by the bρ̃ 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 ((ρ̃ρ̃ correlation)/ z(ρ̃)=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ρ̃ correlation magnitude (≤0.3), which is associated with HB, or a triangle otherwise. At 80min 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ρ̃ 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 160min in the convection zone, there are two stationary points highlighted (E and F). These do not correspond to near-zero values of bρ̃ 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 320min 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ρ̃ 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; Bannister2017) and partially in hybrid systems (where the static background error covariances of traditional variational methods are averaged with ensemble-derived covariances; Bannister2017). It is well known though that the use of ensemble-based covariances leads to rank deficiency and under-sampling (Houtekamer and Mitchell1998; van Leeuwen1999; Hamill et al.2001; Houtekamer and Mitchell2005; Ehrendorfer2007; Zhang and Zhang2012; Houtekamer and Zhang2016), 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/qs in the vertical column, where qs 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, wmax, and (v) condensate, qcmax, 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 Cρ̃/z-b/ rmsCρ̃/z+rmsb 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 u/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=10min, 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 ∼0km and drops at ∼285km, 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).

wmax (red, bottom), qcmax (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 wmax, qcmax, 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 7The geographical variation in various indicators/harbingers of convection at t=10min. 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 qc (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 Sb and is a source term for buoyancy (Eq. 1e) (latent heat is assumed to be exchanged directly only with buoyancy energy). Sb 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, Lv (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 super-saturated (q>qs), in which case q drops immediately (practically over a model time step) to a new qs 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<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, qs, sets a maximum q that can be supported. In real systems qs 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 qs in terms of ABC variables b and ρ̃ as set out in Sect. 2.4. The conditions used to relate b and ρ̃ to qs 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+qc, 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 (Berre2000; Bannister2008; 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

(A1) b t + B u b + A 2 w = - A 2 L v b ( Ev - Co ) ,

which is an explicit form of the buoyancy equation. This is consistent with the basis of our micro-physics parametrisation, Eq. (6), where Δmpb2=tb2Δt=2btbΔt=2bΔmpb. Plugging this into Eq. (6) gives


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 ρ̃b/A2:


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

(A2) ρ ̃ γ t + B u γ = ( ρ ̃ γ ) t + B ( ρ ̃ γ u ) ,

with γ=b2/(2A2). That is,


Making the substitution Eb=ρ̃b2/(2A2) (see Sect. 2.2.2), gives

(A3) E b t + B E b u + ρ ̃ b w = - ρ ̃ L v ( Ev - Co ) .

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

(A4) ρ ̃ L v q t + B ρ ̃ L v q u = ρ ̃ L v ( Ev - Co ) ,

and make the substitution from Eq. (3),

(A5) E l t + B E l u = ρ ̃ L v ( Ev - Co ) .

This is the continuous form of the equation for latent heat energy. The equations for kinetic and elastic energies (Ek and Ee 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, Emoist=Ek+Eb+Ee+El:


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:

(A8) t d x d z E moist = 0 .

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 1/b 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, ρ̃, b, q, and qc (described by Eq. 1) by one time step.

  • 2.

    Calculate the pre-micro-physics saturated specific humidity field, qs(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)>qs(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 b0, allow condensation to happen by solving the following non-linear simultaneous equations,


      and then setting qcmp=qc+q-qmp.

    • b.

      Else if ρ̃Lvq-qs/Eb>Γ, 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

    • c.

      Else do nothing (to avoid unrealistic situations).

  • 4.

    Else if b>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 (qs-q1-e-Δt/τ, qc, and b2/2A2Lv) and find which has the smallest value.

    • a.

      If (qs-q)1-e-Δt/τ 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

    • b.

      If qc 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

    • c.

      Else if b2/2A2Lv is the smallest, then evaporation is limited only by availability of buoyancy energy, Eb=b2/(2A2ρ̃) (there is enough condensate available), by transferring all buoyancy energy to latent heat, leaving

  • 5.

    Overwrite the model's variables, b, q, and qc, with their post-micro-physics values, bmp, qmp, and qcmc 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, T^, and the environmental temperature, Tenv (e.g. Sect. 7.4.1 of Salby1996; CAPE=gLFCLNBT^-Tenv/Tenvdz, 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 w^ is

(C1) w ^ t + B u ^ w ^ = D w ^ D t = D 2 z ^ D t 2 ,

where z^ is the vertical position of the parcel. Equation (C1) is the vertical acceleration of the parcel. We can substitute for w^/t+Bu^w^ using the vertical momentum Eq. (1c). The parcel equation is then

(C2) D 2 z ^ D t 2 = - C ρ ̃ ^ z + b ^ .

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):

(C3) - C ρ ̃ env z + b env = 0 .

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

(C4) D 2 z ^ D t 2 = - C ρ ̃ ^ - ρ ̃ env z + b ^ - b env b ^ - b env ,

where in the second line we have assumed mechanical equilibrium (that ρ̃^=ρ̃env). 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, b^(zi), 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, b^(z0)=b(z0) and q^(z0)=q(z0) respectively. Start with level i=1. By raising the parcel to the next level, the total buoyancy is assumed to be conserved: b^(zi)=b^(zi-1). The total buoyancy is the sum of the reference state and the perturbation, namely b^(zi)=b^0(zi)+b^(zi). By definition (Petrie et al.2017), A2 is the rate of increase in reference buoyancy with height, namely b^0(zi+1)-b^0(zi)=A2(zi+1-zi) in discrete form. Putting the last three equations together gives the parcel's perturbation value at zi:

(C5) b ^ ( z i ) = b ^ ( z i - 1 ) - A 2 ( z i - z i - 1 ) .

The parcel's mixing ratio is also assumed to be conserved: q^(zi)=q^(zi-1).

Should the parcel become super-saturated, the same non-linear equations as used by Hydro-ABC's micro-physics scheme to adjust b^(zi) and q^(zi) towards saturation are used here. This calculates revised values of b^(zi) and q^(zi) (call these b^mp(zi) and q^mp(zi) 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 b^mp(zi) and q^mp(zi) have been found, b^(zi) and q^(zi) are then overwritten respectively by these adjusted values. This process is repeated by replacing zizi+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 b^benv is called the level of free convection, zLFC, and the first level above that where b^benv is called the level of neutral buoyancy, zLNB. CAPE is then computed from the integral:

(C6) CAPE = z LFC z LNB b ^ ( z ) - b env ( z ) d z ,

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:

(C7) CIN = - 0 z LFC b ^ ( z ) - b env ( z ) d z .

A large CAPE is associated with large positive differences, b^(z)-benv(z), between zLFC and zLNB, 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 zLFC not be found (owing to b^<benv for all levels), CAPE would be zero and CIN would be found by integrating over the entire profile. Should a value of zLNB not be found (owing to b^>benv for all levels above zLFC), 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=10min. In the convection zone at 148.5km (left panel), b^>benv 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.5km (middle panel), is where CAPE remains large but where a value of zLFC is found. There is therefore a non-zero CIN associated with this position, but the moist instability still has influence. The other, at 286.5km (right panel), is just past the position where CAPE drops away (see Fig. 7). Here b^<benv 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 420km (not shown).

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


Code and data availability

The Hydro-ABC V2.0 code is written in Fortran 90 and is freely available via a GitHub repository (, last access: 19 October 2023) and Zenodo (, Bannister2022). The Hydro-ABC V2.0 user documentation is available in the same place. Example Hydro-ABC initial condition data (as used in this paper) are available at (Bannister and Zhu2023).

Author contributions

JZ and RNB jointly designed the moisture framework. JZ modified the dry ABC code to include the moist processes described in this paper, ran the model, and produced Figs. 27 and C1. RNB wrote the paper and produced Fig. 1.

Competing interests

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 Hydro-ABC and the two anonymous reviewers for providing valuable and insightful feedback with their reviews.

Financial support

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).

Review statement

This paper was edited by Chiel van Heerwaarden and reviewed by two anonymous referees.


Bannister, R.: A review of operational methods of variational and ensemble-variational data assimilation, Q. J. Roy. Meteor. Soc., 143, 607–633,, 2017. a, b

Bannister, R.: rossbannister/Hydro-ABC_vn2.0: Hydro-ABC_vn_2.0 (Hydro-ABC_vn_2.0), Zenodo [code],, 2022. a

Bannister, R. and Zhu, J.: Ensemble of initial conditions for Hydro-ABC, Zenodo [data set],, 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 ABC-DA system (v1.4): a variational data assimilation system for convective-scale assimilation research with a study of the impact of a balance constraint, Geosci. Model Dev., 13, 3789–3816,, 2020. a

Bannister, R. N.: Balance conditions in variational data assimilation for a high-resolution forecast model, Q. J. Roy. Meteor. Soc., 147, 2917–2934,, 2021.  a

Bannister, R. N., Migliorini, S., and Dixon, M.: Ensemble prediction for nowcasting with a convection-permitting 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 limited-area convection-permitting ensemble, Geosci. Model Dev. Discuss. [preprint],, in review, 2017. a

Bannister, R. N., Chipilski, H., and Martinez-Alvarado, 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,, 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,, 2022. a

Berre, L.: Estimation of synoptic and mesoscale forecast error covariances in a limited-area model, Mon. Weather Rev., 128, 644–667, 2000. a

Clark, P., Roberts, N., Lean, H., Ballard, S. P., and Charlton-Perez, C.: Convection-permitting models: a step-change in rainfall forecasting, Meteorol. Appl., 23, 165–181, 2016. a

Cullen, M. and Davies, T.: A conservative split-explicit integration scheme with fourth-order 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 ensemble-based 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.: Atmosphere-Ocean 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 convective-scale 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.: Distance-dependent 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 cloud-resolving scales, B. Am. Meteorol. Soc., 88, 1783,, 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 convective-scale data assimilation, Tellus A, 69, 1369332,, 2017. a

Lee, J. C. K., Amezcua, J., and Bannister, R. N.: Hybrid ensemble-variational data assimilation in ABC-DA within a tropical framework, Geosci. Model Dev., 15, 6197–6219,, 2022. a

Leung, T. Y., Leutbecher, M., Reich, S., and Shepherd, T. G.: Atmospheric predictability: revisiting the inherent finite-time barrier, J. Atmos. Sci., 76, 3883–3892,, 2019. a

Lorenc, A.: A study of ob monitoring statistics from radiosondes, composited for low-level 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 background-error 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 convective-scale 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 background-error 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 non-hydrostatic toy model for use in convective-scale data assimilation investigations, Geosci. Model Dev., 10, 4419–4441,, 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/TN-556+STR, National Center for Atmospheric Research,, 2021. a

Sun, J., Xue, M., Wilson, J. W., Zawadzki, I., Ballard, S. P., Onvlee-Hooimeyer, 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

Vetra-Carvalho, 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 Feature-Dependent Static Background Error Covariances for Satellite-Derived Humidity Assimilation on Analyses and Forecasts of Multiple Sea Fog Cases over the Yellow Sea, Remote Sens., 14, 4537,, 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 convective-scale numerical weather prediction, B. Am. Meteorol. Soc., 99, 699–710, 2018. a

Zhang, M. and Zhang, F.: E4DVar: Coupling an ensemble Kalman filter with four-dimensional variational data assimilation in a limited-area 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 re-used here as a pragmatic way of deriving temperature from the ABC variables.

Short summary
We describe how condensation and evaporation are included in the existing (otherwise dry) simplified ABC model. The new model (Hydro-ABC) includes transport of vapour and condensate within a dynamical core, and it transitions between these two phases via a micro-physics scheme. The model shows the development of an anvil cloud and excitation of atmospheric waves over many frequencies. The covariances that develop between variables are also studied together with indicators of convective motion.