The “ABC model”: a non-hydrostatic toy model for use in convective-scale data assimilation investigations

. In developing methods for convective-scale data assimilation (DA), it is necessary to consider the full range of motions governed by the compressible Navier–Stokes equations (including non-hydrostatic and ageostrophic ﬂow). These equations describe motion on a wide range of timescales with non-linear coupling. For the purpose of developing new DA techniques that suit the convective-scale problem, it is helpful to use so-called “toy models” that are easy to run and contain the same types of motion as the full equation set. Such a model needs to permit hydrostatic and geostrophic balance at large scales but allow imbalance at small scales, and in particular, it needs to exhibit intermittent convection-like behaviour. Existing “toy models” are not always sufﬁcient for investigating these issues. A simpliﬁed system of intermediate complexity derived from the Euler equations is presented, which supports dispersive gravity and acoustic modes. In this system, the separation of timescales can be greatly reduced by changing the physical parameters


Introduction
Advances in computer power have enabled numerical weather prediction (NWP) models to operate at higher resolutions than has previously been possible. In 2009, the Meteorological Office (Met Office) upgraded the resolution of its Unified Model (UM; Davies et al., 2005) for the UK domain from 12 to 1.5 km ). Resolutions of this degree are expected to resolve the large-and synopticscale features well. Bryan et al. (2003) found that models with resolutions of 100 m are necessary to provide meaningful simulations of convection. Resolutions of O (100 m) are not yet affordable over the UK domain with current computer resources, although research experiments with the UM over smaller domains with 200 m resolution have shown marked benefit (Lean et al., 2008). Models of O ( 1 km) resolution are known as convective-scale models because they are capable of resolving some convection explicitly; thus, they do not require a full convection scheme. For instance, it is possible to explicitly represent features such as thunderstorms O (10 km) and mesoscale convective systems (MCSs) O (10-100 km), though not necessarily resolve their internal structure (e.g. Bryan et al., 2003;Clark et al., 2005;Lean et al., 2008). Convective-scale forecasting can facilitate more accurate and earlier indications of extreme or haz-Published by Copernicus Publications on behalf of the European Geosciences Union. ardous weather, e.g. severe convection (Lean et al., 2008), which is of clear benefit.
As NWP moves towards the convective scale, it is appropriate to examine the data assimilation (DA) scheme underpinning the forecast. The DA process combines meteorological data from a variety of sources including satellites, radar, weather stations, and radiosondes with a previous forecast (a background state) to produce an analysis. The NWP model is then integrated forward from the analysed state. The DA scheme that combines the observed and background data should provide an analysis that is consistent with the observations and the model broadly within the specified "bounds" of the observation errors, the background (prior) state errors, and the model errors (if a model error scheme is incorporated into the DA system). The development of our toy model is a step towards a detailed and technical investigation of the convective-scale DA problem though its utility is not limited to this application.
Convective-scale DA introduces new issues. The errors in the larger-scale flow are still present, but in addition there will be errors on the small scales resolved by the convectivescale model which will have a different correlation structure. A pragmatic solution is to rely on a larger-scale DA system to correct the large-scale errors and thus allow convective-scale DA to focus on small scales. The model introduced in this paper is intended to allow development of methods of assimilating information over all scales. Detailed reviews of the issues are given by Park and Županski (2003), Dance (2004), Sun (2005), and Lorenc and Payne (2007).
The current Met Office operational large-scale DA scheme enforces hydrostatic balance as a strong constraint and exploits geostrophy as a weak constraint in the background error covariance model (Lorenc et al., 2000;Bannister, 2008). The use of the hydrostatic balance relationship is valid for flows where the aspect ratio is much less than 1, e.g. Holton (2004) and Vallis (2006). In regions of convection, the aspect ratio increases, and thus hydrostatic balance may no longer be a good approximation. Vetra-Carvalho et al. (2012) demonstrated that hydrostatic balance breaks down in the UM when it is run at 1.5 km horizontal resolution in regions of convection. At midlatitudes and high latitudes, the geostrophic assumption is accurate for large-scale flows where the Rossby number is small (e.g. Holton, 2004). At the convective scale, the Rossby number is not small, and therefore the use of geostrophic balance is no longer appropriate. It is therefore important that these balances are relaxed in convective-scale DA. Some variational DA methods, such as those termed "EnsVar" (Lorenc, 2013;Liu and Xue, 2016;, use information from an ensemble to represent background error covariance information without, in principle, the need to impose balances via a prescribed background error covariance matrix. These methods though suffer from noise in the sampled error covariance matrix and thus rely on methods to mitigate its effect (namely by localisation, which is known to destroy balances when they are relevant; Kepert, 2009;Bannister, 2015). The sampled (and localised) error covariance matrix in these methods is often hybridised with a prescribed background error covariance matrix (Clayton et al., 2013), which does impose balances. This brings attention back to the validity of such balances when such methods are applied at convective scales and hence to simplified systems where this issue can be studied closely.
Operational systems have to resolve features at both the synoptic and the convective scales, requiring a large number of grid points. Such systems are very expensive to run and are not ideal tools for research purposes. The wide range of timescales means that semi-implicit integration schemes are required for efficiency, e.g. Davies et al. (2005), and the nonlinear coupling between acoustic and gravity waves through the equation of state makes analysing the small-scale behaviour difficult (Thuburn et al., 2002). Thus, it would be useful to have a simplified model which describes a variety of regimes but without the extreme separation of timescales and the full non-linear coupling between acoustic and gravity waves present in the real system. A simplified system that has these properties allows problems such as the convectivescale DA problem to be explored in a practical but physically realistic way.
Perhaps the simplest non-linear model of convection is the well-known Lorenz 63 system (Lorenz, 1963), which describes convection with only three variables. These are (respectively) proportional to the strength of the convective motion, the size of the temperature differences between the upand downwelling air, and the degree of deviation from linearity of the temperature profile. The resulting three ordinary differential equations are easily integrated numerically, but they miss the representation of the complex spatial aspects of the problem required to mirror real forecasting problems. Würsch and Craig (2014) discuss the lack of availability of suitable simplified models of convection for DA research, and they note that people have tended to run full NWP models for this purpose but in idealised settings (see references in Würsch and Craig, 2014, for examples). These models, however, remain complicated and expensive to run. Würsch and Craig (2014) developed a simplified model for purposes of convective-scale DA research. Their model is based on the one-dimensional shallow water model, modified to account for the phase transitions of cloud formation and precipitation -essential processes in the formation of cumulus convection. Although their model has shown to be very useful for this purpose, its one-dimensionality makes it impossible to tackle questions relating to the breakdown of hydrostatic balance and to simulate our inability in practical situations to resolve vertical structures from observations.
The simplified system derived in this paper is intended to be run in vertical slice geometry (longitude/height), so that many fewer degrees of freedom are needed than in an operational three-dimensional system. The equations are modified so that the speed of the acoustic and gravity waves can be controlled, and thus the normally large separations in timescales can be reduced. The equation of state is also modified so that the degree of coupling between the acoustic and gravity waves is reduced. The modifications are designed so that energy is conserved in the equations, which is necessary for realistic behaviour. In order to study the dynamically related breakdown of balance, no moisture is included, but intermittent convection-like behaviour is still seen (e.g. via gravity wave breaking). These simplifications permit large-scale balanced flows and sporadic small-scale nonhydrostatic flows (including convection) to coexist within the framework of a simplified and practical model.
The structure of this paper is as follows: Sect. 2 provides a derivation of the toy model equations which are analysed in terms of a scale analysis and energy conservation properties. Section 3 describes the numerical implementation of the model. Section 4 provides a linear analysis of the equations. Section 5 shows the results of an idealised integration which illustrates how the model can be used to simulate different flow regimes. Section 6 provides a summary and some concluding remarks. Future work will exploit this model in testing different approaches to convective-scale DA, as piloted in . For reference, Table 1 summarises the symbols used throughout this paper.

Derivation of the model equations
The model is derived from the compressible 3-D Euler equation (Eq. 1); see, e.g. Holton (2004), Pielke (2001), and Vallis (2006): Equation (1a) is the momentum equations, where t is time; u = (u, v, w) comprises zonal (u), meridional (v), and vertical (w) components; p is pressure; g is the acceleration due to gravity; and ρ is density. The f -plane assumption is made, and k is the vertical unit vector. Equation (1b) is the compressible mass continuity equation. Equation (1c) is the adiabatic thermodynamic equation, where θ is potential temperature. Equation (1d) is the equation of state, where p 00 = 1000 hPa, κ = R/c p is a constant, with c p the specific heat capacity at constant pressure, and R the gas constant for dry air.
From this set of equations, we wish to construct a toy model that has large-scale geostrophically and hydrostatically balanced flow, permits intermittent convective-like behaviour, and is of practical use for investigating issues that arise in the convective-scale DA problem (e.g. it is cheap to integrate).

Modifications to the 3-D Euler equations
In order to derive a model with the properties outlined above, Eq. (1) is modified in two stages. Firstly, a set of physically based approximations is made, and secondly a set of "toy model" simplifications is made. The latter set does not attempt to replicate the real system; rather, it is intended to retain desired physical characteristics of the real system but simplify the computational implementation. In order to simplify the system, it will be assumed that the model is periodic in the zonal direction and homogeneous in the meridional direction (i.e. the variables are functions of longitude, height, and time only).

Physically based modifications
The variables are decomposed such that they have a basic state and perturbation component, e.g. Pielke (2001): (2) Here, applies to any model variable except θ (for θ , see below). The basic state (subscript 0) is a function of height only, and the perturbation (primed) is a function of longitude (x), height (z), and time (t). Potential temperature contains also a constant reference value (subscript R): The wind components u, v, and w have zero reference state values; therefore, the prime notation is dropped for the winds. For convenience, explicit reference to the arguments x and z will be dropped in much of the following derivation. The basic state is assumed to satisfy hydrostatic balance: and the equation of state is The Brunt-Väisälä frequency, N, is defined as The pressure gradient terms in Eq. (1a) are represented with Eq.
(2), products of perturbations are neglected, and it is assumed that ρ 0 ρ in the momentum equations. A buoyancy variable, b = b 0 (z) + b (x, z), is introduced for convenience; it is related to θ by Scaled density (Eq. 14) b Buoyancy (Eq. 7) q Passive tracer concentration Long. pos. variables scaled by horiz. length scales Horiz., vert. wavenumbers n x , n z Horiz., vert. wavenumber indices σ , σ R , σ g , σ g , σ a , σ a Wave frequency, specifically Rossby, gravity, acoustic Horiz., vert. domain sizes c g , c a Gravity and acoustic group speeds bu eff Effective buoyancy Combining these physically based approximations gives the following equations:

The "ABC model" modifications
It is desirable to reduce the stiffness of the system of Eq. (8) so that it can be integrated explicitly with a time step that is not too small. The following toy model modifications are made so that the toy equations retain the basic properties desired; i.e. they are geostrophically and hydrostatically balanced on the large scale but permit intermittent convectionlike behaviour on the small scale that is unbalanced. The modifications are as follows.
1. We control the gravity waves by replacing N by the tunable parameter, A (which has units of s −1 ). This is the pure gravity wave frequency (Sect. 4.4).
2. We control the acoustic waves by multiplying the divergent term of the compressible continuity equation by the dimensionless parameter B (where 0 < B ≤ 1). To ensure energy conservation (Sect. 2.3), it is required that B also multiplies the advective components of the momentum and thermodynamic equations. Acoustic waves can have frequencies that are normally much higher than gravity waves, but choosing a small B can help to reduce the acoustic wave frequencies.
The effect of these parameters on the wave speeds will be demonstrated by numerical linear analysis in Sect. 4.6. The acoustic and gravity waves in the real atmosphere are coupled through the equation of state (Thuburn et al., 2002). This coupling can be reduced by using a linearised and simplified equation of state. Linearising Eq. (8f) about the basic state gives where we have used θ R + θ 0 ≈ θ R (see Appendix A). This is used in two ways to give modifications 3 and 4 below.

4.
Secondly and separately, for the purposes of simplifying the equation of state, neglecting buoyancy perturbations in the linearised equation of state (Eq. 9) gives This is a means of decoupling gravity and acoustic waves. Further, setting gives the simplified equation of state: where C is taken to be a global constant and has units of N m kg −1 = m 2 s −2 . The quantity √ BC is the pure sound wave speed in this system (Sect. 4.5).
5. Reference density ρ 0 is taken to be a constant and not a function of height.
Combining modifications 1 to 6 gives the final form of the toy model equations: Note that Eq. (15d) conserves mass following the flow modulated by B, i.e. Bu, but total mass remains conserved (Sect. 2.3.1). We also include the following tracer transport equation for diagnostic purposes: where q is the tracer concentration. Note that the advection term is not multiplied by B in Eq. (16) (as B will be generally chosen as B ≤ 1, we allow advection of the tracer to have its full effect so that tracer transport can be seen over an integration of a few hours). We refer to these simplified equations as the "ABC model", reflecting the three tunable parameters.

Scale analysis of the ABC model
A scale analysis of Eq. (15) is performed by nondimensionalising the equations using characteristic values. Our scale analysis deviates from standard analyses in two ways: (i) we allow different characteristic length scales for each variable (in the horizontal and vertical), and (ii) we do not assume incompressibility (see below for further explanation). For the characteristic values, we set u = Uu * , v = Vv * , w = Ww * , ρ = P ρ * , ρ ∼ 1, and b = B b * . For the characteristic horizontal length scales, we set (respectively, for each variable) Often in scale analyses the characteristic vertical speed, W, is written in terms of other characteristic variables by using the incompressible continuity equation in a 2-D (longitude-height) system ∂u/∂x + ∂w/∂z = 0. Scaling this gives W ∼ UL V w /L H u . We do not use this relation, as some of the flows considered are highly compressible.
Using these definitions in the system of Eq. (15) and introducing the Rossby number, Ro = U/f L H u , the aspect ratio, the vertical-to-zonal wind ratio, W U = W/U, and the meridional-to-zonal wind ratio, V U = V/U, gives the following non-dimensionalised equations: When the first three terms of Eq. (17a) and (17b) are small (often achieved with small Ro), the geostrophic relationships emerge. Expressed back in terms of the dimensional variables, they are Under similar circumstances, Eq. (17c) defines the hydrostatic relationship. Expressed back in terms of the dimensional variables, it is

Conservation of mass and energy
As the toy model of Eq. (15) is no longer based on standard thermodynamics, we must demonstrate that it forms a physically reasonable set. To this end, we now show that it conserves mass and energy.

Conservation of mass
Noting the definition in Eq. (14), multiplying the continuity equation, Eq. (15d), by the constant ρ 0 gives the equation for the evolution of density perturbations. Adding the zerovalued term ∂ρ 0 /∂t then produces the equation for the evolution of density: ∂ρ/∂t + B∇ · (ρu) = 0. Since the model uses periodic boundary conditions zonally and zero vertical wind conditions at the top and bottom boundaries (Sect. 3.2), the divergence theorem shows that the equations conserve mass:

A useful "identity" used to demonstrate conservation of energy
Dividing the continuity equation shown in Sect. 2.3.1 by ρ 0 gives the equation for ρ evolution: ∂ ρ/∂t + B∇ · ( ρu) = 0. Using this equation and expanding ∂( ργ )/∂t + B∇ · ( ργ u), for an arbitrary time-and space-varying scalar field γ , we find Equation (20) is treated as an identity in the forthcoming energy analysis.

Buoyant energy
Multiplying the thermodynamic Eq. (15e) by ρb /A 2 and using Eq. (20) where the perturbation buoyant energy, E b , is

Elastic energy
Multiplying the continuity Eq. (15d) by C ρ /B, we find where the perturbation elastic energy, E e , is

Total combined energy and its conservation
Adding Eqs. (23), (24), and (26) shows that the combined energy, Integrating Eq. (28) over the whole domain for the total combined energy gives This toy model is set up to have periodic boundary conditions in the x direction, no variation in the y direction, and zero vertical wind at the top and bottom boundaries (see Sect. 3.2). The divergence theorem then leads to conservation of total combined energy (see Appendix B): 3 Numerical implementation of the ABC model Now that a physically reasonable set of toy model equations has been formed, we provide the details of how they are treated numerically.

Model discretisation
The toy model uses a similar grid to that of the southern UK (SUK) version of the UK Met Office's Unified Model (UM) but with some differences given below. In the horizontal, the Table 2. Upper and lower boundary conditions of each prognostic model variable; z = 0 is the lower boundary position and z = L z is the upper boundary position.
Lower Upper SUK model covers a domain of 540 km in longitude and 432 km in latitude with a resolution of 1.5 km on an Arakawa C grid. In the vertical, it has 70 vertical levels up to approximately 40 km on an irregularly spaced Charney-Phillips grid (Lean et al., 2008). The toy model grid is shown in Fig. 1. The differences from the SUK are that the toy model is periodic in the zonal direction, is homogeneous in the meridional direction, and uses regularly spaced vertical levels up to a lid of about 15 km. The toy model uses only 60 levels (level spacing δz ≈ 250 m) and has flat orography.
This grid is a natural discretisation of the equations which does not require a significant number of interpolations. There are approximately 10 5 variables in the state space of the toy system.

Boundary conditions
The vertical boundary conditions that we use are summarised in Table 2. At the lower boundary, the horizontal winds are zero (no-slip conditions) and the vertical wind is zero to conserve total mass and energy. At the upper and lower boundaries, the vertical derivative of density is zero. For the equations to have the capability to support hydrostatic balance, Eq. (19) implies that b should be zero at the vertical boundaries. At the upper boundary, the horizontal winds are chosen to maintain consistency with the boundary conditions of p and b through thermal wind balance, and the vertical wind is again zero to conserve total mass and energy.

Time integration scheme
The time integration is evaluated using a split-explicit, forward-backward scheme (Cullen and Davies, 1991), and here we give a description of this scheme applied to the system of Eq. (15). The forward-backward scheme operates over a time step t and comprises two stages: an adjustment stage and an advection stage. The adjustment stage operates over a sub-time step δt, where δt = t/n and n is typically a small positive integer (in this implementation, n = 2). The adjustment stage contains two parts: the forward part and the backward part. Let t be the time at the start of the t time step and let t i be shorthand for t + iδt. The following is a description of the ith sub-time step.
In the forward part of the forward-backward scheme, the momentum and thermodynamic equations are evaluated omitting the advective terms. The u and v equations are considered simultaneously to find the adjustment due to the Coriolis and pressure gradient terms. Then, the w momentum and b equations are dealt with simultaneously to find the adjustment due to buoyancy, pressure gradient, and vertical wind. The forward part of the adjustment stage gives an implicit approximation to u, v, w, and b at the next sub-time step.
The equations for u (Eq. 15a) and v (Eq. 15b), omitting the advective terms, are discretised as Solving the system of Eq. (31) for u(t i+1 ) and v(t i+1 ) gives where α f and β f are defined by The equations for w Eq. (15c) and b Eq. (15e) omitting advective terms are discretised as Solving Eq. (34) for w(t i+1 ) and b (t i+1 ) gives where α A and β A are defined by Equations (32a), (32b), (35a), and (35b) are the discretised forms of the split-explicit equations that are evaluated in the forward part of the forward-backward scheme in the adjustment stage. The spatial derivatives are left in continuous form but are discretised in the numerics using standard centred finite differences.
In the backward part of the forward-backward scheme, the continuity Eq. (15d) is evaluated using the wind and buoyancy data calculated in the forward part, i.e.
The term in brackets on the right-hand side is equal to ∇ · ( ρ(t i )u(t i+1 )) but has been expanded in Eq. (37) to allow the second term to use the upstream gradient of ρ(t i ). After integration of n steps over the full t, the value of p is known and the values of the variables u, v, w, and b are known but without the effect of advection.

Advection stage
The advection stage advects the fields u, v, w, and b calculated in the adjustment stage using the sub-time-stepaveraged windsū andw, which are taken to be valid over the full t. Let φ be any of u, v, w, or b ; then, the advection step is given by whereū = (ū,w) T . As with Eq. (37), the upwind gradient of φ is computed in Eq. (38). Note that for the tracer advection, φ = q, Eq. (38) is used with B = 1.

Overall properties
The spatial derivatives evaluated by the forward-upstream scheme are first-order accurate (Press et al., 2007) and the time integration which utilises a split-explicit and forwardbackward scheme is also first-order accurate (Ames, 1958;Gadd, 1978). The stability of the forward-backward scheme increases the time steps which are permitted by the Courant-Friedrichs-Lewy (CFL) criterion (Ames, 1958). The splitexplicit scheme has been used in early implementations of the UK Met Office's NWP model due to its ability to conserve mass (Gadd, 1978;Cullen and Davies, 1991).

Linear/normal mode analysis of the ABC model
In this section, a normal mode analysis of the toy model equations is performed. This follows a similar procedure used for the shallow water equations in Sect. 6.4 of Daley (1991) and in Sect. 2.4 of Cullen (2006). In this procedure, the model equations are first linearised to permit a mathematical analysis of small amplitude perturbations. The normal modes are the characteristic solutions of the linear equations. Each mode has a characteristic frequency and wavelength (hence, the analysis is done in terms of spectral modes) and is a specific combination of different variables. There are three different types of normal mode solution, namely Rossby-like modes (whose normal mode patterns obey geostrophic and hydrostatic balances), gravity modes (buoyancy-driven modes which have characteristic horizontal divergence), and acoustic modes (the most rapidly oscillating modes which are made up of compression waves). The linear/normal mode analysis allows us to probe the dispersion relations (mode frequencies as a function of wavenumber) and the balanced/unbalanced character of the linear modes. For simplicity, this analysis is performed on a continuous periodic domain of sizes L x and L z .

Linearisation
The non-linear model of Eq. (15) is linearised about a state of rest (see Appendix A). It is convenient to write the model equations in terms of velocity potential, χ (describing the divergent part of the horizontal flow), and stream function, ψ (describing the rotational or solenoidal part). The Helmholtz theorem (see, e.g. Salby, 1996, which is generically given by "horizontal wind = ∇ h χ + k × ∇ h ψ", where k is the vertical unit vector and ∇ h = (∂/∂x, ∂/∂y, 0)) allows the horizontal wind to be written as Noting the lack of y dependence in the ABC model, this gives u = ∂χ /∂x and v = ∂ψ/∂x. The linearised model equations are then ∂ ∂t ∂ ∂t

Normal mode analysis
Now take the following functional dependence for a particular frequency σ , horizontal wavenumber k, and vertical wavenumber m: Substituting Eq. (41) into Eq. (40) and expressing the resulting set of equations in matrix form gives This is an eigenvalue equation where L is a real and symmetric matrix (due to the choice of factors in Eq. 41) and thus will have real eigenvalues. For each distinct choice of horizontal and vertical wavenumber (k, m), L has five eigenvalues, denoted σ R , σ g , σ g , σ a , and σ a , where σ R = 0, σ g = −σ g , and σ a = −σ a .
The three distinct modes are the Rossby-like mode (subscript R), two inertia gravity modes (g and g ), and two acoustic modes (a and a ). The algebraic form of the R mode is simple and is discussed in Sect. 4.3 below, but the forms of the remaining modes are very complicated and thus are considered only firstly in "pure" forms (Sect. 4.4 and 4.5) and then numerically in the wave speed analysis (Sect. 4.6).

The Rossby-like mode
The normalised R mode has σ R = 0 and is www.geosci-model-dev.net/10/4419/2017/ Geosci. Model Dev., 10, 4419-4441, 2017 where K = (AL z [kBC + L x f 2 ] + L 2 x f 2 m 2 C)/(L 2 x f 2 m 2 BC). This mode, as we shall show, supports geostrophic balance defined by Eq. (18a) and (18b). Firstly, relation (Eq. 18a) in terms of the variables defined in Eq. (41) and for wavenumber k is which is consistent with Eq. (45). Secondly, and trivially, relation (Eq. 18b) is equivalent to ∂χ /∂x = 0, which is also consistent with Eq. (45). There is no vertical wind associated with the R mode. There remains a buoyancy component for this mode to support hydrostatic balance defined by Eq. (19). Relation (Eq. 19) in terms of the variables defined in Eq. (41) and for wavenumbers k, m is which is consistent with Eq. (45).

The pure gravity waves
Following Kalnay (2002), pure gravity waves can be investigated by neglecting rotation and pressure perturbations (by Eq. 13 density perturbations are therefore neglected too). We anticipate that the gravity waves will be sensitive to A given that A is related to the static stability parameter N (the Brunt-Väisälä frequency). Under these conditions, Eq. (43) has two eigenvalues, σ g = A and σ g = −A, representing the pure gravity wave frequencies. In the limit of A = 0, no gravity waves are supported.

The pure acoustic waves
Following Kalnay (2002), pure acoustic waves can be investigated by neglecting rotation, gravitation, and stratification (i.e. set f = 0, g = 0, A 2 = 0, and b = 0). Under these conditions, Eq. (43) has three eigenvalues, σ = 0 (which is an incompressible mode that does not interest us here), σ a = √ BC (k/L x ) 2 + (m/L z ) 2 , and σ a = −σ a , the latter two representing the pure acoustic wave frequencies. The pure acoustic wave speed in the horizontal (for example) is ∂σ a /∂k, which becomes √ BC in the small-scale limit. In the limit that B = 0 or C = 0, the system becomes incompressible and no acoustic waves are supported.

Wave speed analysis experiments
In Sect. 4.4 and 4.5, we demonstrated how the pure gravity and acoustic waves depend upon the parameters A, B, and C. The analysis there was simplified (by explicitly neglecting processes that are not directly associated with gravity and acoustic waves, respectively) in order to derive analytical expressions. Here, we look at the gravity and acoustic wave speeds in a more detailed way without making the approximations made before. These reveal the normal modes of the linearised system of Eq. (40) (see, e.g. Thuburn et al., 2002), which now includes rotation, gravitation, and stratification. We show how the wave speeds behave in the linearised system, and as a function of wavenumber, and of parameter values. To reduce the stiffness of the system, we would like the speeds of the gravity and acoustic modes to have value ∼ U or ∼ V, the characteristic speeds of the horizontal wind components, and thus the results of this subsection are important for choosing parameter values for suitable model runs.
The standard values of the parameters that we use for this section are A = 0.02 s −1 (estimated from a typical value of the Brunt-Väisälä frequency), B = 1.0, C = 10 5 m 2 s −2 (estimated from initialising data), and f = 10 −4 s −1 , and for simplicity, periodicity is assumed in the x and z directions. Figure 2 shows the horizontal group speeds for the gravity (c g = ∂σ g /∂k; Fig. 2a) and acoustic (c a = ∂σ a /∂k; Fig. 2b) waves as a function of the integer index, n x (characterising the horizontal wavenumber k = 2n x π/L x ) for a range of parameter values (the integer index, n z , characterising the vertical wavenumber m = 2n z π/L z , is fixed at n z = 3). Note that these k and m values are slightly different from those used in Sect. 4.1 to 4.5.
Gravity waves in the approximated system are found to be stationary (Sect. 4.4), but gravity waves in the full system are not; see Fig. 2a. There is a strong sensitivity of c g to A (larger A, faster gravity waves), and the fastest gravity waves have large horizontal and large vertical scales (small n x and n z ). There is also a sensitivity of c g to BC, which is especially evident at large vertical scales and over large and intermediate horizontal scales (not shown, but see Sect. 5.3.3).
Acoustic waves in this system have different characteristics to the gravity waves in many respects. Acoustic waves can be much faster, but their speed may be controlled via the strong sensitivity of c a to BC, and the fastest acoustic waves have small horizontal and small vertical scales (large n x and n z ); see Fig 2b. It is at these small scales that the acoustic waves saturate to the value √ BC as found in Sect. 4.5. The sensitivity of c a to A is weak for the smaller values of A tested but moderate for the largest value of A tested (not shown).
The ability of the parameters A, B, and C to change the speed of both the horizontal gravity and acoustic waves has been demonstrated in Fig. 2. The buoyancy frequency parameter A primarily controls the gravity waves, and the product BC primarily controls the acoustic waves. The vertical gravity and acoustic wave speeds respond in a similar way to the parameters as the horizontal waves (not shown). In addition to modifying the acoustic wave speeds, B and C have other, separate effects: B slows advection round the domain, and C influences the hydrostatic and geostrophic balance relationships (see Sect. 2.2). It is permissible to alter the value of only one or any combination of these parameters depending on the required result. These results are used in the next section to help choose suitable parameter values.

Reference parameters
In this section, some desired dynamical characteristics that are seen in the real atmosphere are demonstrated in this simplified setup. It is required that the model mimics the multiscale behaviour of the real atmosphere, i.e. displays hydrostatic and geostrophic balance on the large scale while permitting imbalance and intermittent convective-like behaviour on the small scale, while allowing an explicit solver. The results from the linear analysis of Sect. 4.6 give a taste of how the wave speeds depend on A, B, and C, but the values that we settle on as reference values are A = 2 × 10 −2 s −1 , B = 10 −2 , and C = 10 4 m 2 s −2 . Figure 3a shows the frequencies and the magnitudes of the horizontal and vertical wave speeds for the gravity and acoustic waves for these reference parameters. The acoustic wave frequencies in Fig. 3a are always higher than those of the gravity waves (the latter, with an upper bound of A). The frequencies of the gravity and acoustic waves for n z = 3 (left) are of the same order, but the acoustic wave frequencies for the extreme case of n z = 59 (right) have much higher values by more than an order of magnitude. These are classic dispersion curves for these modes in the atmosphere (e.g. Fig. 14.9 of Salby, 1996), and they allow us to estimate that the highest frequency that the model will encounter with these parameters is ∼ 0.25 s −1 (4 s period, from the right panel of Fig. 3a). The parameter values though will (in Sect. 5.3) be changed by an order of magnitude each. Over these extended parameter values, the maximum wave frequency is found to be ∼ 1.6 s −1 (0.625 s period). This allows us to set the time steps of our model (Sect. 3.3.1), which we choose as t = 0.1 s and δt = 0.05 s. We use these values for all experiments.
The ability to control speeds in order to make the gravity and acoustic wave speeds comparable is more effective than the ability to control frequencies. Comparing, for instance, Figs. 3b and 2 shows how the gravity and acoustic wave speeds have been reduced to comparable values (a maximum of 10 m s −1 with the reference parameters, compared with a maximum of 1000 m s −1 for the parameters tested for Fig. 2). The speed of 10 m s −1 applies in the horizontal (Fig. 3b) and in the vertical (Fig. 3c). As well as using a t that is smaller than the shortest wave period (see above), we can use the maximum wave speed to also check that t is consistent with the grid spacings ( x and z) via the CFL condition (Durran, 1999). Typically, this states that the Courant number, Co, satisfies 0 ≤ Co ≤ 1 for stability of the numerical solution. Given that x = 1500 m, z = 250 m, and t = 0.1 s, the Courant number is Co = t × (u max / x + u max / z) = 0.1 × (10/1500 + 10/250) ≈ 0.005, which satisfies the CFL condition for the reference parameters. The maximum wave speed for the other parameter sets studied in Sect. 5.3 is ∼ 30 m s −1 , which still results in a small Courant number.

Idealised initial conditions
The model was first initialised with idealised initial conditions to ensure that the model behaves reasonably with the reference parameter values. In this run, the initial conditions are zero for all variables apart from ρ , which takes the form of the Gaussian described in the caption (Fig. 4a). The ρ and (u, v) fields for up to 6 h are shown in Fig. 4. In the real atmosphere, such a positive density perturbation induces anticyclonic motion as geostrophic adjustment develops, and a similar response is seen in the toy model (the "vertical" components of the arrows represent meridional wind, which is out of the page on the right and into the page on the left). After 3 h (Fig. 4b), the horizontal wind is significantly divergent, indicating that the ρ perturbation is being dissipated by gravity waves which smooth out the initial perturbation, whose maximum value has reduced to about a third of its original value. After 6 h (Fig. 4c), the flow is mainly rota- Figure 3. Gravity and acoustic wave properties for the reference parameters A = 2 × 10 −2 s −1 , B = 10 −2 , and C = 10 4 m 2 s −2 . The panels are frequencies (a) and the magnitudes of the horizontal (b) and vertical (c) wave speeds. In panels (a) and (b), values are a function of horizontal wavenumber, n x ; the left column is for n z = 3, and the right column is for n z = 59. In panel (c), values are a function of vertical wavenumber, n z ; the left column is for n x = 10, and the right column is for n x = 350. tional (there is a weak convergent flow near the centre of the domain) and the ρ perturbation has moved to the boundaries. Figure 4 can be used to verify the wave speeds determined by linear analysis. Consider Fig. 4b, where the edge of the feature has propagated approximately 80 km over the 3 h. This gives an approximate horizontal gravity wave speed of ∼ 7 m s −1 , which is around the maximum horizontal gravity wave speed found from the linear analysis in Fig. 3b.

Intermittent convection-like behaviour
Convective motion in the atmosphere is difficult to model and to assimilate, as it is often intermittent and associated with small-scale divergence. In the real atmosphere, it is usually driven by latent heating, but our simple model is dry, and thus we rely on other processes such as wave breaking to drive such motion. Intermittent convection-like motion is a desirable property of our model in order be a useful system to study convective-scale data assimilation. An indication of the presence of convection is vertical motion, and vertical motion necessarily indicates an imbalance; see Eq. (45). We look at the w field from an integration of the model firstly with the reference parameters. The initial conditions of the model were created from the following procedure.
-Take values of u, v from a latitude/height slice of an output the Met Office's convective-scale (1.5 km grid) UM (this is the same model used by the Met Office during the 2012 Olympics, as indicated by Golding et al., 2014, and has the same horizontal resolution and grid staggering as our model). These fields are adjusted to eliminate the discontinuity imposed by the periodic boundary conditions 1 .
-Calculate ρ by integrating the geostrophic balance equation (Eq. 18a) on each level.
x is the horizontal distance from the western boundary, = 150 km is the relaxation distance, is the size of the discontinuity in u or v (i.e. the magnitude of the difference in u or v between the western and eastern boundaries in the raw UM data), and δ is the magnitude of a typical increment of u or v between neighbouring grid boxes. This procedure is performed separately for each vertical level.
-Calculate w from the continuity equation for zero threedimensional divergence.
Each variable is then incremented (independently for each level) so that its horizontal mean is zero, and finally ρ is set as ρ = 1+ ρ . The model's initial conditions are then nearly balanced (the incrementing will disrupt the hydrostatic balance slightly), and unbalanced motion (including convection-like behaviour) then develops. Figure 5 shows w over a 6 h integration of the model using the reference parameters. The initial conditions in Fig. 5a show vertical winds that are of relatively small scale in the horizontal, with elongated structures over the lowest 5 km or so in the middle of the domain. These are of course not generated by this model but are derived from the UM data. An indication of the kind of behaviour that this model is capable of generating is shown in Fig. 5b and c, for 3 and 6 h into the forecast, respectively. The most striking aspect of the w field at 3 and 6 h is that the scales of the features are even smaller than those at t = 0. Additionally, the magnitudes of w are smaller with very small regions of moderate values especially in the eastern part of the domain. The similarity in these qualitative aspects of Fig. 5b and c shows that this kind of behaviour is not merely transient. We regard these plots as indicators of intermittent convection-like behaviour, which is studied further below.

Systematic exploration of model behaviour over parameter space
In order to test novel approaches to data assimilation, it is desirable to run the model in different flow regimes, which we investigate by varying the parameters A, B, and C systematically, each over a 3 h model run. We are particularly interested in understanding how the parameter values affect the degree of convection and of imbalance, and we start this investigation by introducing the diagnostics for the reference parameters.

Reference parameters
We settle on four kinds of diagnostic for each parameter set, which are shown in Fig. 6 for the reference parameters. Figure 6a is the vertical wind speed, and Fig. 6b is the effective buoyancy, bu eff . The latter is defined as the (non-constant) stability bu eff = ∂ b 0 (z) + b (x, z) /∂z, where ∂b 0 (z)/∂z = A 2 . A positive (negative) effective buoyancy indicates statically (un)stable air, so negative and small positive values suggest convective activity (contours for negative values of w and bu eff are dashed). Figure 6c is the distribution of tracers after 3 h. The tracers were initialised at t = 0 on a grid of 20 points distributed throughout the domain (small rectangular regions in Fig. 6c) and the distribution after 3 h provides an indication of the history of the wind behaviour. These fields are labelled with the minimum, maximum, and root mean squared values. Figure 6d indicates the degree of relative geostrophic imbalance (blue lines and left scale) and hydrostatic imbalance (red lines and right scale) averaged over the domain, at half-hour intervals over the integration. These quantities are found, respectively, using Eqs. (18a) and (19) to give where rms indicates the root mean squared value of the quantity in brackets over the domain. The fields ρ , v, and b are filtered before computing these diagnostics by removing scales (i) below 100 km (to give the solid lines in Fig. 6d), (ii) below 10 km (to give the dashed lines), and (iii) below 1 km (i.e. unfiltered, to give the dotted lines). This gives us an indication of how the degree of imbalance is affected by scale.
There are variations of upward and downward vertical motion over the domain (Fig. 6a), but there are no regions that are specifically more convectively active than others. The bu eff diagnostic is fairly uniformly small over most of the domain (Fig. 6b) but does have more variability in the uppermost 5 km of the domain where it is weakly negative in a thin layer at 14 km (sandwiched between two strongly stable layers) during this snapshot. There is a small amount of disturbance of the tracer field after 3 h (Fig. 6c).
The Rossby number is estimated to be small (Ro ∼ 0.06), and the geostrophic imbalance is found to be moderate for the reference run (Fig. 6d), which stabilises to around 0.45 when only large scales are present, but higher, to around 0.7 to 0.8, when smaller scales are included (see the blue lines and the left-hand scale in Fig. 6d). The hydrostatic imbalance also increases as the scales shorten (see the red lines and the righthand scale in Fig. 6d) but is much lower than the geostrophic imbalance (0.025 for the smallest scales). By estimating the magnitudes and length scales of the fields in this run, the scale analysis in Sect. 2.2 does show that the last two terms in Eq. (17a) and (17c) to be much larger than the other terms by about 3 and 6 orders of magnitude, respectively.
In order to check that the linear analysis of Sect. 4 is relevant to the non-linear model integration, Table 3 compares the horizontal gravity wave speeds found from the linear analysis applied to the reference values, to the propagation speed of anomalies in the horizontal divergence field found from time sequences of model output (not shown; we refer to this as "feature tracking"). Horizontal divergence, ∇ h (u, v) = ∂u/∂x + ∂v/∂y = ∂u/∂x = ∇ 2 h χ (where χ is the velocity potential used in Sect. 4.1), is often associated with gravity wave activity. As is shown in Sect. 4, the wave speeds are dispersive in this system as so are a function of the horizontal and vertical wavenumbers. The relevant wavenumbers are found by looking at the characteristic length scales of the fields (second and third columns in Table 3), which (d) Geo. and hydro. imbalance for all experiments conducted in this paper correspond, respectively, to horizontal wavenumber index n x = 3 and vertical wavenumber index n z = 2. The corresponding horizontal gravity wave speeds from the linear analysis are given in the fourth column, and the measured speed from the feature tracking is given in the fifth column. The values are comparable (in this and the other experiments), suggesting that the linear analysis is indeed a useful guide to the behaviour of the non-linear model.
Energy in the continuous system of equations was proven to be exactly conserved in Sect. 2.3, but the discretisation and numerical integration scheme introduce errors which will lead to non-conservation. Figure 7a (black solid line) shows that these errors do lead to a small loss of energy over the 3 h (less than half of a percent of the initial energy), which we assume is acceptable.

Changes to the parameter A
Recall that the parameter A controls the gravity wave frequency and speed. In this section, two 3 h integrations are done: one with A decreased by an order of magnitude (A−; Fig. 8, left panels) and one with A increased by an order of magnitude (A+; Fig. 8, right panels).
A− appears to result in more active w values than in the reference run, and A+ appears to result in slightly less active w values (Fig. 8a). The effective buoyancy (Fig. 8b) has more structure than in the reference run, with bands of lowered bu eff appearing in A− (with patches of slightly negative bu eff in the lower part of the domain which are too small to show as contours in the left plot of Fig. 8b), while A+ has no negative values at all. The increased vertical motion in A− is seen in the tracer fields (Fig. 8c), which have been transported more vertically in A− and slightly less vertically in A+ than the reference run. These results make physical sense given that A controls the static stability of the fluid.
There are differences and similarities in the geostrophic and hydrostatic imbalances between A−, A+, and the reference run (Fig. 8d). In terms of geostrophic balance, A− is seen to be more balanced, and A+ is slightly less balanced than the reference run, and they all have similar Rossby numbers ∼ 0.06. This observation seems counter-intuitive since we would expect the higher gravity wave speed of A+ to lead to a state of more, not less, balance (although the adjustment timescale -indicated by the continuous blue curves in Figs. 6d and 8d -is shorter for A+). The effect of A on the gravity wave speeds is discussed in Sect. 4 and is confirmed further in Table 3, which compares the linear analysis and feature tracking analysis approaches. A− does indeed show a significant decrease in the gravity wave speed, and A+ shows an increase, and the degree of change is roughly consistent between the two approaches. In terms of hydrostatic balance, A− is slightly less balanced and A+ is much less balanced than the reference run. A− is consistent with our expectation, but A+ is not. Table 3. Summary of the characteristic horizontal and vertical length scales found for each of the experiments (estimated after 3 h model integrations), the horizontal gravity wave speeds corresponding to these scales as found from the linear analysis, and an estimate of the gravity wave speed found by tracking features in the horizontal divergence field (∂u/∂x).  Even though we have ensured that the Courant number is small for all experiments (Sect. 4.7), we ask the question: is it possible that these unexpected balance results are artefacts of numerical imprecision? Is it possible that, as A is changed by an order of magnitude for A− and A+, the changes in speed of the waves (especially gravity) or changes in the vertical transport may result in a loss of precision and hence lead to unreliable balance diagnostics?
The first issue (imprecision due to the change in the wave speeds and thus the Courant number) is first studied by investigating the dependence of the balance diagnostics to changes of the integration time step t. It is found that t can be increased up to ∼ 50 times (i.e. to ∼ 5 s, thus increasing the Courant number) before there are any noticeable changes to any of the balance diagnostics (not shown). This suggests that any loss of precision due to a too-long t used in the main experiments can be ruled out. The Courant number can also be increased by reducing the grid size. It was found that the grid size can be halved (in the horizontal and vertical directions simultaneously) without qualitatively changing the balance results. In these higher-resolution experiments, it is found (not shown) that the measured degree of imbalance reduces in all runs, but the degree of imbalance relative to each experiment is unchanged (i.e. A− still has more geostrophic balance than the reference run, and A+ has slightly less, and this qualitative agreement is true also for the experiments to be discussed in later sections when run at higher resolution). This suggests that, although the higher-resolution runs have a different precision to the standard resolution runs, the fact that the results are qualitatively the same indicates that our results are probably robust.
The second issue (imprecision due to change in the vertical motion) is studied by noticing that there is an association between greater vertical motion, and increased energy loss (we assume that energy loss is in turn associated with imprecision). The less statically stable A− run results in more vertical motion and more energy loss than the reference run (4 % loss over the 3 h, solid red line in Fig. 7a) and the more stable A+ run results in less vertical motion and less energy loss (0.2 % loss, dotted red line). Can this mean that variations in imprecision lead to unreliable diagnostics and hence are responsible for the counter-intuitive balance results summarised above? The higher-resolution runs mentioned above are more successful at conserving energy (the higher-resolution counterparts of Fig. 7a are shown on the same scale in Fig. 7b). A−, e.g. loses ∼ 2.5 times less energy at the higher-resolution than A− at standard resolution (solid red lines). This suggests that the numerical integrations (and hence the balance diagnostics) for the higher-resolution runs are more reliable. Since the high-resolution balance diagnostics are qualitatively the same as those of the standard resolution results shown, this suggests that our results are not anomalous. This is not a definitive conclusion, but it does demonstrate some robustness of our results. The counterintuitive balance results though remain unexplained at this stage, although it should be noted that the balance diagnostics themselves may only be meaningful to compare within the same system of parameter values rather than between different systems.

Changes to the parameter B
Recall that the parameter B (with C) controls mainly the acoustic wave speed. Two further 3 h integrations are done: one with B decreased by an order of magnitude (B−; Fig. 9, left panels), and one with B increased by an order of magnitude (B+; Fig. 9, right panels).
B− results in similar magnitude w values as the reference run, and B+ results in a slight increase, but there is little change to the structures of the w field (Fig. 9a). The effective buoyancy is largely unaffected by the changes in B (Fig. 9b). A similar story applies to the tracers for B−, but the tracers for B+ do show increased vertical transport (Fig. 9c), which is consistent with the larger root mean squared w for B+. B− (B+) has slightly more (less) geostrophic and hydrostatic balance, and the Rossby number remains small, ∼ 0.07 (∼ 0.11). As discussed in Sect. 5.3.2, one would normally assume that a faster wave speed, as in the B+ run, would result in more balanced fields, but this is not the case here. The scale analysis in Sect. 2.2 reveals though that it is the product BRo, rather than just Ro, that is the quantity that scales terms that knock the system out of balance, and BRo is smaller (larger) in B− (B+) than in the reference run.
Changing B affects mainly the acoustic wave speed, but it can effect the gravity wave speed too. This is shown by the linear analysis and feature tracking analysis in Table 3, which shows a consistent decrease in the horizontal gravity wave speed for B− and an increase for B+. Although this effect of B on the gravity wave speeds is large at these scales, its effect on the acoustic waves is much greater (for the scales in Table 3, the acoustic waves vary between 0.1 m s −1 for B− and 30 m s −1 for B+). Note that for smaller vertical scales the influence of B on the gravity wave speed is much smaller.
Changing the B parameter has a dramatic effect on the energy conservation (Fig. 7a), where B− yields the least energy loss of all experiments (on the scale used, the numerical loss of energy is indistinguishable from a perfectly conservative scheme -solid blue line in Fig. 7a), but B+ results in one of the most erroneous runs (an 8 % loss in energy over 3 h, dotted blue line). Reducing the time step of the integration to 1/10 of the value used in the main runs improves the energy conservation only marginally (not shown), but halving grid lengths (in the horizontal and vertical) improves the conservation (e.g. just over 2 % loss over 3 h for B+, dotted blue line in Fig. 7b).

Changes to the parameter C
Recall that the parameter C (with B) controls mainly the acoustic wave speed. Two further 3 h integrations are done: one with C decreased by an order of magnitude (C−) and one with C increased by an order of magnitude (C−). The initial conditions for C− and C+ each differ from those used before, as the procedure used to generate balanced initial conditions described in Sect. 5.2 from UM data depends on parameter C.
The C− and C+ results are not shown because they are virtually indistinguishable from the B− and B+ runs, respectively (including the relative balance results). There are two differences though. The first is that ρ is scaled by C −1 (when C is decreased (increased) by an order of magnitude, ρ (not shown) is increased (decreased) by an order of magnitude compared to the B− (B+) runs, with the field structures remaining the same). This is seen in the scale analysis Eq. (17), where C and P (the characteristic value of ρ ) always appear together as a product. This is how the C− and C+ runs maintain the same level of geostrophic and hydrostatic balances as the B− and B+ runs, respectively. The second difference is seen in the numerical scheme's energy loss. The B+ run (at the original model resolution) loses energy significantly (8 %), but the C+ run loses less (1 %, dotted green line of Fig. 7a), and the C− run loses even less (around 0.1 % , solid green line). This is another beneficial effect of introducing the B parameter (i.e. the same value of a particular desired acoustic wave speed √ BC can be achieved by decreasing B and increasing C).

Conclusions
A set of simplified mass-and energy-conserving equations has been derived which allows the gravity and acoustic wave characteristics to be controlled with three parameters: A (the pure gravity wave frequency), B (the modulation of the divergent term in the continuity equation, and of the advection terms in other equations), and C (the proportionality constant for the toy model's equation of state). The term √ BC is the pure small-scale acoustic wave speed. The introduction of B allows the acoustic wave speed to be reduced so that it is comparable to the gravity wave speed, hence allowing explicit integration schemes to be used to approximate the so- lution of the equation set (such as the split-explicit, forwardbackward scheme used here). The linearised equations support a zero-frequency Rossbylike mode and dispersive gravity and acoustic modes. The system is shown to behave in a way that reflects aspects of the atmosphere, namely geostrophic adjustment, convective behaviour influenced by buoyancy, and scale-dependent geostrophic and hydrostatic imbalances. Some of the results concerning the effect of changing A on the degree of balance are counter-intuitive. Although we believe that the numerical results are robust, further work could be done to understand these results. It may be simply that the diagnostics (Eq. 48) cannot be compared between sets of parameters. The model has no water vapour, which simplifies the scheme considerably (although water vapour and moist processes could be added if required). The energy is not perfectly conserved in practice, which is due to the finite discretisation of the model (and to the choice of integration scheme), although numerical energy loss is assumed to be acceptable in most runs.
The purpose of developing this model is to facilitate research into ways of modelling the background error covariance matrix (B) used in convective-scale data assimilation. The B matrix is normally modelled with guidance from large-scale dynamics, namely that geostrophic balance is dominant, and hydrostatic balance is exact. These assumptions are probably not applicable at convective scales (as shown by Berre, 2000, Bannister et al., 2011, Vetra-Carvalho et al., 2012, Bannister, 2015, and as we have seen here, where more imbalance is present at the smaller scales). A key idea which will be explored in a forthcoming paper is to use the normal mode structure of the linearised equations to define the B matrix rather than relying on imposed balances. It is hoped that this will have physically appropriate structures and the correct degree of balance at different scales (preliminary work has been done by . Code and data availability. The model is written in Fortran 90, and the plotting code is written in Python. This software is open source and freely available on a GitHub repository (Petrie et al., 2017). The initial conditions used to start the model runs studied in this paper are available from the same repository. and the "surface" S therefore represents the boundary of the model (in a plane, the divergence theorem reduces to an area integral on the left-hand side and a closed line integral on the right-hand side, and is equivalent to Green's theorem). The meaning of A depends upon the application.
To prove conservation of mass, let A = ρu, in which case Eq. (B1) gives where L x and L z are the length and height of the model, respectively, and the four terms on the right-hand side represent contributions from the four sides of the model's domain (starting from the lowermost/west-most point and integrating clockwise around the domain). The model has zero vertical wind values at the top and bottom boundaries (Table 2), which removes the second and fourth terms on the right-hand side. Furthermore, swapping the integration limits on the third integral (which introduces a minus sign to that term) and noting that the fields are periodic in the horizontal (ρ(L x , z)u(L x , z) = ρ(0, z)u(0, z)) leads to the first and third terms cancelling. The right-hand side is therefore zero. The equation describing the evolution of the total mass is where the last line follows from the second line using the mass continuity equation given in Sect. 2.3.1, and the final zero results from Eq. (B2) being zero as shown above. This proves conservation of total mass in the ABC model. To prove conservation of energy, the second and third terms in Eq. (29) must be shown to be zero. For the second term, use Gauss' divergence theorem with A = B(E k +E b )u, and for the third term use A = C ρ ρu . Since each of these is proportional to u, the same arguments used for mass conservation can be used to show that the second and third terms of Eq. (29) do not contribute. This proves that the rate of change of total energy integrated throughout the domain, as in Eq. (29), is zero, and therefore energy is conserved.