Dynamics of an idealized fluid model for investigating convective-scale data assimilation

Abstract

An idealized fluid model of convective-scale numerical weather prediction, intended for use in inexpensive data assimilation experiments, is described here and its distinctive dynamics are investigated. The model modifies the rotating shallow water equations to include some simplified dynamics of cumulus convection and associated precipitation, extending and improving the model of Würsch and Craig. Changes to this original model are the removal of ad hoc diffusive terms and the addition of Coriolis rotation terms, leading to a so-called 1.5-dimensional model. Despite the non-trivial modifications to the parent equations, it is shown that this shallow water type model remains hyperbolic in character and can be integrated accordingly using a discontinuous Galerkin finite element method for nonconservative hyperbolic systems of partial differential equations. Combined with methods to ensure well-balancedness and non-negativity, the resulting numerical solver is novel, efficient and robust. Classical numerical experiments in the shallow water theory, such as the Rossby geostrophic adjustment and flow over topography, are reproduced for the standard shallow water model and used to highlight the modified dynamics of the new model. In particular, it exhibits important aspects of convective-scale dynamics relating to the disruption of large-scale balance and is able to simulate other features related to convecting and precipitating weather systems. Our analysis here and preliminary results suggest that the model is well suited for efficiently and robustly investigating data assimilation schemes in an idealized ‘convective-scale’ forecast assimilation framework.

Keywords:

• shallow water models
• data assimilation
• numerical weather prediction
• convection dynamics
• discontinuous Galerkin finite element method

1. Introduction

Numerical weather prediction (NWP) models solve non-linear partial differential equations (PDEs) that describe atmospheric motions on many scales, whilst parameterizing unresolved processes at the smaller scales as a function of the resolved state. In the context of NWP, data assimilation (DA) involves incorporating meteorological observations in the forecast model in a dynamically consistent manner to provide the ‘optimal’ initial condition for a forecast of the future atmospheric state, taking into account errors in both observations and previous forecasts (Kalnay, 2003). Optimality of the initial state is crucial in such a highly non-linear system with limited predictability. Indeed, significant gains in the accuracy of NWP can be attributed to improvements in assimilation algorithms and observing systems.

Until recently, operational NWP models were running with a horizontal resolution larger than the size of most convective disturbances, such as cumulus cloud formation, which were accordingly parameterized. Despite the coarse resolution leaving many ‘subgrid’-scale dynamical processes unresolved, there has been a great deal of success in weather forecasting owing mainly to the dominance of large-scale dynamics in meteorology (Cullen, 2006). Variational DA algorithms have successfully exploited this notion that atmospheric dynamics are close to a balanced state (e.g. hydrostatic and semi-/quasi-geostrophic balance), resulting in analysed states and forecasts that remain likewise close to this balance (Bannister, 2010).

Increasing computational capability has led in recent years to the development of high-resolution models at national meteorological centres in which some of the convective-scale dynamics are explicitly (or at least partially) resolved (e.g. Done et al., 2004; Baldauf et al., 2011; Tang et al., 2013). This so-called ‘grey-zone’, the range of horizontal scales in which convection and cloud processes are being partly resolved dynamically and partly by subgrid parameterizations, presents a considerable challenge to the NWP and DA community (Hong and Dudhia, 2012). Current regional NWP models are running at a spatialgridsize on the order of 1 km with future refinement inevitable, and smaller scale processes are known to interfere with DA algorithms based on the aforesaid balance principles (Vetra-Carvalho et al., 2012). As such, high-resolution NWP benefits hugely from having its own DA system, rather than using a downscaled large-scale analysis (Dow and Macpherson, 2013).

To aid understanding of and facilitate research into such large and complex operational forecast assimilation systems, simplified models can be utilized that represent some essential features of these systems yet are computationally inexpensive and easy to implement. This allows one to investigate and optimize current and alternative assimilation algorithms in a cleaner environment before making insights or considering implementation in a full NWP model with real observing systems (Ehrendorfer, 2007). Systems of ordinary differential equations (ODEs), such as the L63 model (Lorenz, 1963) and its successors Lorenz, 1986; Lorenz, 1996; (Lorenz and Emanuel, 1998; Lorenz, 2005), continue to be the basis for numerous DA studies (e.g. Neef et al. (2006); nee2009; sub2012; Bowler et al. (2013); Fairbairn et al. (2014)). They provide chaotic dynamics on a range of scales yet their low dimensionality means that they are computationally cheap and easy to implement in an idealized forecast assimilation system. The gap in the complexity of such ODE models and the primitive equation models of operational forecasting is, however, too large. Shallow water type models attempt to bridge this gap. They capture interactions between waves and vortical motions in rotating stratified fluids and have received attention in DA research for the ocean and atmosphere (e.g. Zhu et al., 1994; Žagar et al., 2004; Salman et al., 2006; Stewart et al., 2013). Here, we extend and analyse a modified shallow water model originally proposed by Würsch and Craig (2014) for DA research (herein WC14).

Convective (cumulus) clouds are characterized by highly buoyant, unstable air that accelerates upwards in a localized region to significant heights (see, e.g. Houze, 1993a). If the air reaches a sufficient height, precipitation forms and subsequently falls through the convective column, reducing the buoyancy and turning the updraft into a downdraft (along with associated effects from latent heat release). The model of WC14 captures some aspects of this life cycle of single-cell convection, while following the classical shallow water dynamics in non-convecting and non-precipitating regions. The binary ‘on–off’ nature of convection and precipitation is inherently difficult to resolve in NWP models, requiring highly non-linear functions that pose further issues for convective-scale DA algorithms. Thus, the inclusion of switches, in the form of threshold heights, provides a relevant analogy to operational NWP and is an important aspect of the modified model.

The difference between the model proposed here and that of WC14 is twofold. First, incorporating a meridional velocity component and Coriolis terms means that dynamics associated with rotating fluids, such as geostrophy, are present in the model. Second, and more importantly, the diffusion terms used to stabilize the model of WC14 are removed, resulting in a hyperbolic system of PDEs. Accordingly, the model can be integrated robustly using a discontinuous Galerkin finite element method (DGFEM) for hyperbolic systems, cf. Rhebergen et al. (2008), coupled with the method of Audusse et al. (2004) to ensure well-balancedness. This novel framework ensures non-negativity of the layer depth and the ‘rain mass fraction’ variable; it is also more versatile for analysis than using a leading order finite volume approach. While the DG formulation includes higher order discretization in space, the model and methodology is demonstrated here at leading order in a series of test simulations chosen to illustrate the model’s distinctive dynamics.

The purpose of this paper is to introduce the model and the above-mentioned numerical solver, and consequently to investigate its distinctive dynamics. We consider this of scientific interest in itself, but also as a prerequisite for its use in DA experiments. In the next section, the physical motivation and mathematical description of the model are given. A key aspect of the model is that, despite the modifications, it remains hyperbolic, thus permitting the use of a powerful class of numerical methods for such PDE systems. Sections 3 and 4 introduce a new scheme for the numerical integration and illustrate the modified dynamics of the model with respect to the classical shallow water theory. We conclude with a summary of the key features of the dynamics of the model and some comments on its intended use in an idealized forecast assimilation framework.

2. Model description

2.1. Classical shallow water

Shallow water (SW) flows are ubiquitous in nature and their governing equations have wide applications in the dynamics of rotating, stratified fluids. The shallow water equations (SWEs) are considered a useful tool for modelling dynamical processes of the Earth’s atmosphere and oceans. They approximately describe inviscid, incompressible free surface fluid flows under the assumption that the depth of the fluid is much smaller than the wavelength of any disturbances to the free surface, i.e. a fluid in which the vertical length scale is much smaller than the horizontal length scale.

Interesting dynamical features of the SWEs are gravity waves, vortical motions and shocks. Models based on the SWEs capture the interaction between fast gravity waves and the slowly varying geostrophic vortical mode. Gravity waves are known to play an important role in the initiation of atmospheric convection, particularly in the presence of orography, suggesting a model based on the SWEs is appropriate for investigating convective-scale DA. By definition, shock waves occur wherever the solution is discontinuous. Such discontinuities in the model variables (or their spatial derivatives) are mathematical idealizations of severe gradients, akin to fronts in an atmosphere. As such, propagation of shock waves in the model can be thought of as the propagation of atmospheric fronts (Parrett and Cullen, 1984; Frierson et al., 2004; Bouchut et al., 2009).

The standard shallow water model on a rotating Cartesian f-plane (2dRSW) in which dynamical variables do not depend on one of the spatial coordinates (here the y-coordinate, so that ) can be written as (see, e.g. Zeitlin (2007)): (1a) (1b) (1c) where is the space- and time-dependent fluid depth, is the prescribed underlying topography (so that is the free surface height), u(x, t) and v(x, t) are velocity components in the zonal x- and meridional y-direction, f is the Coriolis parameter (typically in the mid-latitudes) and g is the gravitational acceleration. This system of equations, together with specified initial and, where appropriate, boundary conditions, determine how the flow evolves in time. The effective pressure p(h), following the terminology of isentropic gas dynamics, has the standard form: . It is useful to introduce the equations in this form to illustrate the modifications described in the next section.

Physically, this model extends the one-dimensional SWEs by adding transverse flow v and Coriolis effects. The existence of transverse flow with no variation in the y-direction means the model should not be considered one- or two-dimensional, but rather one-and-a-half-dimensional (e.g. Bouchut et al., 2009). This set-up offers more complex dynamics associated with rotating fluids (e.g. geostrophy) than a purely 1D model whilst remaining computationally inexpensive, a crucial factor for a ‘toy’ model.

2.2. Modified shallow water

The model introduced by WC14 extends the one-dimensional SWEs to mimic conditional instability and include idealized moisture transport via a ‘rain mass fraction’ r or, alternatively, ‘precipitated water fraction’. We use similar physical concepts and argumentation here but employ a mathematically cleaner approach without diffusive terms which results in a hyperbolic system of PDEs. Other ‘moist’ SW models have been developed for atmospheric dynamics on the synoptic scale, perhaps most famously by Gill (1982) and more recently by, e.g. Bouchut et al. (2009), Zerroukat and Allen (2015). Our interest in the WC14 model stems from its simplicity in incorporating convective motions, namely rapidly evolving updrafts, downdrafts and idealized precipitation effects, without the need for explicitly considering temperature and other thermodynamic properties.

Heuristically, atmospheric moist convection can be thought of as a two-fluid problem, in which one fluid can transform itself into another simply through vertical displacement (Stevens, 2005). It is this concept that is attractive in the WC14 model and that we seek to capture in our idealized ‘convective–scale’ model: the single-layer SWEs (1) are modified when the height of the fluid crosses certain thresholds. When the fluid exceeds these threshold heights, which can be seen as switches for the onset of convection and precipitation, different mechanisms kick in and alter the classical shallow water dynamics. In these modified regions, the behaviour of the flow is transformed from the standard shallow water dynamics to a simplified representation of cumulus convection. The changes to the governing equations are as follows. The mass (1a(1a) ) and hv-momentum (1c(1c) ) equations are unchanged. The hu-momentum Equation (1b(1b) ) is altered by the effective pressure and the inclusion of a ‘rain water mass potential’, . To close the system, an evolution equation for the ‘rain mass fraction’ r is required, including source and sink terms (2d(2d) below). The modified rotating shallow water (modRSW) model is described by the following equations: (2a) (2b) (2c) (2d) where P and Q are defined via the effective pressure by: (3a) (3b) with denoting the derivative of p with respect to its argument h, and:(4)

The constants () and (dimensionless) control the removal and production of rain, respectively, () converts the dimensionless r into a potential in the momentum equation and controls the strength of the feedback, and (m) are critical heights pertaining to the onset of convection and precipitation. For and r initially zero, it is clear that the model reduces exactly to the classical shallow water model (1).

Figure 1. Schematic of the pressure term P(h; b) in (3): the modified pressure above the threshold is lower than the standard pressure , thus forcing the fluid to rise where .

The modification to the standard SWEs first occurs when the free surface height exceeds the threshold in (3). The fundamental dynamics of cumulus convection are the dynamics of buoyant air: air motions in all convective clouds emerge in the form of vertical accelerations that occur when moist air becomes locally unstable (i.e. less dense) than its environment (see, e.g. Markowski and Richardson (2011b)). Initiation of deep convection requires that air parcels reach their level of free convection (LFC), the height at which the air parcel achieves positive buoyancy due to latent heat release from condensation, thus forcing it further upwards through the atmosphere. Associated with the rapid ascent (and subsequent descent) of air in a localized region is the adjustment of the mass field in and around the cloud due to perturbations of a characteristic pressure field (Houze, 1993a). Thus, it can be expected intuitively that buoyancy cannot be instigated without a simultaneous disturbance to the pressure field (Houze, 1993b). This mechanism is exemplified by the threshold height which can be thought of as the LFC: exceedance of forces fluid in that region to rise by modifying the pressure terms (3). The pressure at a given height above , namely , is lower than the standard pressure p(h) at that same height (see the schematic in Fig. 1). Owing to this relative reduction in pressure, the fluid experiences a reduced restoring force due to gravity and should therefore rise. Thus, the changes to the pressure terms (3) instigate positive buoyancy and a representation of conditional instability.

Model ‘rain’ is produced (i.e. the ‘rain mass fraction’ r increases) when the fluid exceeds a second threshold height (higher to ensure that precipitation forms at some time after the onset of convection), in addition to positive wind convergence (). This convergence condition is synonymous with the upward displacement of an air parcel from the surface and subsequent convective updraft. In three-dimensional models, horizontal moisture convergence, , for some moisture field q and horizontal velocity , is often used to parameterize bulk convection and is also a forecasting diagnostic for the initiation of deep moist convection (Markowski and Richardson, 2011a). It is well known that moisture convergence is correlated with horizontal wind convergence ; thus, the condition is conceptually credible and ensures that air is still rising for precipitation to form. The term in (2d(2d) ) and (4(4) ) controls how much ‘rain’ is produced and is a tunable parameter. Once there is model ‘rain’ in the system, it feeds back to the hu-momentum Equation (2b(2b) ) via the term, and precipitates via a linear removal term involving the tunable parameter . In nature, as precipitation forms and subsequently falls through a cloud, it reduces and eventually overcomes the positive buoyancy, thus turning an updraft into a downdraft. This gradient term imitates this effect and can be considered as a contribution to the pressure, increasing it locally where there is model ‘rain’ present. As a consequence, we expect the presence of ‘rain’ to enhance the restoring force and therefore suppress the updraft, eventually leading to the collapse of the convective column and limiting the growth of convection in the model. The strength of this feedback is controlled by the tunable parameter , with higher values leading to enhanced suppression.

We stress that what we refer to as model ‘rain’, namely the model variable r, is an artefact of our model and clearly not the same as nature’s rain. In essence, it is the mass fraction of precipitated water in the system and engenders some highly non-linear effects of precipitation in a simplified modelling environment. Accordingly, one can think of as the mass fraction of precipitable water (i.e. the whole column is a source for precipitation), with the total mass being conserved. Thus, removal of ‘rain’ by the sink term in (2d(2d) ) does not remove mass from the system, rather it transfers it between precipitated and precipitable water so that it may precipitate again at a later time. This is not a realistic feature, however, it means that there is a continual source of model ‘rain’ and so the model does not work for a limited time only, a crucial point when considering its use in idealized DA experiments. The nature of these parameterizations, viz. the tunable parameters , , and , are by their construction ad hoc, but as demonstrated by our numerical simulations in Section 4 and Würsch and Craig (2014), they provide a plausible way to parameterize the idealized transport of moisture in the modified system.

The essential thermodynamic properties central to moist convection (namely latent heat release due to condensation) are in some sense hidden in our model (2). It should be noted though that the effects of these properties are not absent, rather they are modelled indirectly by the modification of the pressure terms. This achieves some simplified dynamics of convection associated with buoyancy, as demonstrated in the numerical experiments in Section 4, without the explicit inclusion of temperature and condensation. Together with the idealized precipitation process, we argue that our simplified approach provides a ‘toy’ model with interesting dynamics nonetheless.

2.3. Hyperbolicity

Hyperbolic systems of PDEs arise from physical phenomena that exhibit wave motion or advective transport. Such systems have a rich mathematical structure and have been extensively researched from both an analytical (e.g. Whitham, 1974) and numerical perspective (e.g. LeVeque, 2002). The classical SWEs are a well-known example of a system of hyperbolic PDEs, being a special case of isentropic gas dynamics. Here, we show that the modRSW model remains hyperbolic despite the non-trivial modifications and non-conservative products (NCPs).

A system of n PDEs is hyperbolic if all the eigenvalues , , of its Jacobian matrix are real and the Jacobian is diagonalizable (i.e. its eigenvectors form a basis in ). To show hyperbolicity (and facilitate numerical implementation in the next section), the modRSW model (2) is expressed in non-conservative vector form:(5)

where:(6)

and P, Q and given by (3) and (4(4) ), respectively. It is non-conservative in the sense that the system cannot be written in divergence form, i.e. the NCP cannot be expressed in terms of a flux function (there is no function such that ). The Jacobian matrix of the system (5(5) ) is given by:(7)

and its four eigenvalues are:(8)

Clearly, are real. Since is non-negative and P(h, b) is non-decreasing (hence ), the term under the square root is non-negative. Hence, are real and, since there are repeated eigenvalues, we conclude that the modRSW model is (weakly) hyperbolic.

Hyperbolic systems are often studied analytically via the method of characteristics. This leads to a transformation of variables into a new set of Riemann variables that propagate along characteristic curves in (x, t)-space (Whitham, 1974). Although this is in principle possible for the modRSW model, the complexity of the system results in abstruse expressions for Riemann variables, offering little insight analytically. But as the prime purpose here is to provide a physically plausible numerical forecast model for conducting idealized DA experiments, further Riemann analysis is neglected. However, one aspect relating to the wave speeds (determined by the eigenvalues) deserves a further comment. It is well known that waves travelling through saturated regions of convection slow down (e.g. Harlim and Majda (2013)), and simplified models of a moist atmosphere should naturally reflect this. For example, the SW model of Bouchut et al. (2009) for a large-scale moist atmosphere has lower wave speeds in ‘moist’ regions compared to dry regions. For comparison, the eigenvalues of the classical shallow water system (1) are:(9)

For the modRSW model (2), is smaller when , since then , and smaller for when is sufficiently small (specifically, less than gh), both relative to the standard shallow water case with . Hence, modulation of the wave speeds by convection is captured in the model too.

Properties such as moist enthalpy and potential vorticity conservation, which are present in moist SWEs derived from the vertically averaged primitive equations in pressure coordinates (see Bouchut et al. (2009)), are absent in our model. However, since our goal is to provide a ‘toy’ model that exhibits some basic features of convecting and precipitating weather systems for use in inexpensive and idealized DA experiments, these properties are of secondary importance.

3. Numerical formulation

3.1. Methodology

There exists a powerful class of numerical methods for solving hyperbolic problems, motivated by the need to capture shock formation in the solutions, a consequence of non-linearities in the governing equations. Efficient and accurate finite volume schemes for systems of conservation laws are very well developed (e.g. LeVeque, 2002; Toro, 2009). For shallow water models, there are well-balanced schemes that deal accurately with topography and Coriolis effects, maintaining steady states at rest and non-negative fluid depth h(x, t) (Audusse et al., 2004); bou2007. However, the nature of the modRSW model (namely, the presence of NCPs including step functions) requires careful treatment beyond the typical methods for conservation laws. The DGFEM developed by Rhebergen et al. (2008) offers a robust method for solving systems of non-conservative hyperbolic PDEs of the form (5(5) ) and is developed here for the modRSW model (2).

The first step of any finite element method is to convert the PDE of interest into its equivalent weak formulation using the standard test function and integration approach. However, the presence of NCPs in the governing equations complicates this somewhat because the weak solution in the classical sense of distributions does not exist when the solution becomes discontinuous (Rhebergen et al., 2008). To overcome the absence of a weak solution for systems of equations of the form (5(5) ), Rhebergen et al. (2008) employ DLM theory (after Dal Maso, LeFloch, and Murat; Dal Maso (1995)) for NCPs which defines an NCP as a bounded measure in such a way to enable the weak solution to be defined. This is achieved by considering a single NCP , where g is a smooth function but u may admit discontinuities, and defining a smooth regularization of the discontinuous u:(10)

where is the Dirac measure at the discontinuity and is a Lipschitiz continuous path connecting the model states across the discontinuity, an artefact of the regularization. In DGFEM, the computational states are generally continuous on each element but discontinuous across an element boundary. It is in this context that the framework afforded by the DLM theory (and culminating in (10(10) )) appears naturally in the weak formulation and subsequent discretization.

Here, we provide a summary of the scheme developed for the modRSW model; further technical material is appended and is referenced accordingly. For full details of the methodology for general systems, including the key theorems employed from the DLM theory, the reader is directed to Rhebergen et al. (2008).

3.2. Discretization

The one-dimensional flow domain is divided into N open elements for with nodes/edges . Element lengths may vary. Formally, one can define a tessellation of the N elements :(11)

where overbar denotes closure of an element with its boundary , i.e. . This simply means that the elements cover the whole domain and do not overlap. The space-DGFEM weak formulation is obtained by (i) multiplying each equation of the system (5(5) ) by an arbitrary test function , generally continuous on each element but discontinuous across an element boundary; and (ii) integrating (by parts) over each element and summing over all elements. The space discretization is achieved by replacing the exact model states and test functions w by approximations in terms of polynomial basis function expansions, with the order of the polynomials determining the order of the scheme.

In the following, repeated i, j-subscript indices are used for the summation convention with denoting components of vectors, k-subscript denotes values in element and L, R-superscript denotes limiting values to the left/right of an element edge (e.g. , and ). In one space dimension and considering cell only at a given t, the weak form reads (from Equation A11 in Rhebergen et al., 2008):(12)

where and are the numerical flux terms given by: (13a) and the NCP flux through an element edge is:(13b)

Here, is the standard HLL numerical flux,(14) is the ij-th element of the matrix , and are the fastest left- and right-moving signal speeds (cf. (8(8) )) in the solution of the Riemann problem, determined by the eigenvalues of the Jacobian of the system: (15) (16a) Since the goal here is a toy model for DA research, it is preferable to keep the scheme as computationally efficient as possible and acknowledge higher order accuracy as of secondary importance. The lowest order space discretization uses piecewise constant basis functions to approximate the model states (so-called ‘zero-th’ order (DG0); equivalent to the finite volume method in one dimension). In principle, the topography b can be treated as a model variable ( with ) such that the topographic source term in (6(6) ) is then treated as an NCP. However, hitherto less well-known issues with well-balancedness for DG0 discretizations with varying topography mean this approach is unsatisfactory; instead, we discretize the topographic source term directly using the established method of Audusse et al. (2004), resulting in a well-balanced scheme at lowest order that efficiently preserves non-negativity of fluid depth h and rain hr. It is first necessary to isolate the topographic source term in (6(6) ) from the other terms pertaining to rotation and removal of rain: . Then, is discretized as a standard linear extraneous forcing term and via the method of Audusse et al. (2004).

Using piecewise constant basis functions and alternately in each element (since the test function is arbitrary), the semi-discrete space-DGFEM scheme for element reads:(16b)

where are reconstructed states to the left and right of node , and is the discretized topographic source term. See Appendix 2 for further details pertaining to these reconstructions, , and the scheme of Audusse et al. (2004).

Figure 2. Time evolution of the height profile for the standard shallow water case I (left), case II with convection and no rain with (middle) and case III with convection and rain for finite (right). Non-dimensional simulation details: .

The contribution from DLM theory (10(10) ) is apparent throughout the flux terms, as is its dependence on the regularization path connecting the left state to the right state. Here we employ a linear path . It is clear from (14(14) ) that in the absence of NCPs ( for all i, j) the numerical flux reduces exactly to the standard flux. However, for , the NCP contributions of the form in (10(10) ) must be calculated. The NCP flux (14(14) ) for the modRSW model is:(17)

where contains the contribution from the NCP integral expressions:(18)

and , are expressions containing Heaviside functions associated with the instantaneous thresholds and . The average of a quantity is denoted by and the jump of a quantity across a node is denoted by . The full derivation of the NCP flux is given in Appendix 3.

4. Numerical experiments and dynamics

This section presents some numerical experiments which have been chosen to highlight the dynamics of the modified rotating shallow water model (2) compared to those of the classical model (1). The experiments are based on: (i) a Rossby adjustment scenario and (ii) non-rotating flow over topography, both of which have a rich history in shallow water theory including known exact steady-state solutions. To illustrate the effect that exceeding the threshold heights has on the dynamics, a hierarchy of model ‘cases’ is employed:

• Case I: always (effectively setting ). The model (2) reduces to standard RSWEs (1) if initially.

• Case II: always, but may exceed . This is considered a ‘stepping stone’ to the full model to isolate the effect of the first threshold exceedance. Thus, given exceedance and the consequent modification to the gradient of the pressure (3a(3a) ), we expect the fluid to be forced upwards (a ‘convective updraft’).

• Case III: may exceed both (and ). This is the full modRSW model with convecting and rain processes to be used for idealized convective-scale DA research.

For the modRSW model to have credibility as a shallow water-type model, it is crucial that it reproduces, in case I, known results of the standard SWEs. The existence of exact steady-state solutions thus provides a benchmark to test this and the solutions can be used as reference states to compare the subsequent modifications introduced by cases II and III. We expect simulations in cases II and III to display markedly different behaviour compared to the ‘dry’ system, and will elucidate these distinctive dynamics with reference to the physical basis described in Section 2.2.

The non-dimensionalized equations (Appendix 1) are implemented on a domain of unit length using the mixed NCP-Audusse numerical scheme derived in the previous section and the forward Euler time discretization. Neuman outflow (cf. LeVeque, 2002) boundary conditions are applied in all simulations. This means that the fluid is allowed to leave the flow domain in a physically consistent manner, essentially setting the domain to be infinitely large. In this case, the required information is typically extrapolated from the interior solution. This is achieved by extending the computational mesh to include so-called ‘ghost’ elements and . Values in these elements are set at the beginning of each time step in a way that takes into consideration the boundary conditions, and the updating algorithm is then exactly the same in every element. Care needs to be taken when implementing outflow conditions to ensure that the specified boundary information does not contaminate the interior solution. Outgoing waves should propagate out of the domain without generating spurious reflections from the artificial boundary. The most simple, yet sufficient, approach uses a zero-order extrapolation, that sets the same value to the model variables in the first and last two elements of the domain at the start of each time step. This does not prevent reflections completely but any artificial effects are negligible compared to the solution. Further simulation details for each experiment are given in figure captions and the main text.

4.1. Rossby adjustment (unbalanced velocity)

The following experiment, motivated by Bouchut et al. (2004), explores Rossby adjustment dynamics in which the evolution of the free surface height is disturbed from its rest state by a transverse jet, i.e. fluid with an initial constant height profile is subject to a localized v-velocity distribution. In order to adjust to this initial momentum imbalance, the height field evolves rapidly, emitting inertia gravity waves and shocks that propagate out from the jet and eventually reach a state of geostrophic balance. The shape of the initial velocity profile of the jet v(x) is that employed by Bouchut et al. (2004):(19)

and the initial conditions are , and . The bottom topography b is zero throughout the domain.

Figure 3. Hovmöller plots for the Rossby adjustment process with initial transverse jet: case I (left), II (middle) and III (right). From top to bottom: h(x, t), u(x, t), v(x, t) and r(x, t). Non-dimensional simulation details: same as Fig. 2.

Figure 4. Evolution of h and r for the Rossby adjustment process with initial transverse jet: case I (left), II (middle) and III (right). Top row: Hovmöller plots for h. Subsequent rows: profiles of h (black line; left axis) and r (blue line; right axis) at different times denoted by the dashed lines in the top row. Non-dimensional simulation details: same as Fig. 2.

Figure 5. Hovmöller plots for the Rossby adjustment process with initial transverse jet, highlighting the conditions for the production of rain: case III. From left to right: , , and r(x, t). Non-dimensional simulation details: same as Fig. 2.

Figure 6. Top row: Hovmöller diagram plotting the evolution of the departure from geostrophic balance : light (deep) shading denotes regions close to (far from) geostrophic balance. Subsequent rows: profiles of fv (red) and (black) at different times denoted by the dashed lines in the top figure. For case I (left), II (middle), and III (right). Non-dimensional simulation details: same as Fig. 2.

Figure 7. Flow over topography (, and ): profiles of , b (black; left y-axis), exact steady-state solution for the SWEs (red dashed; as derived in Appendix 4) and rain r (blue; right y-axis) at different times: case I (left), II (middle) and III (right). The dotted lines denote the threshold heights . Non-dimensional simulation details: .

Figure 8. Hovmöller plots for flow over topography (), highlighting the conditions for the production and subsequent evolution of rain: case III. From left to right: , and r. Non-dimensional simulation details: same as Fig. 7.

Figure 9. Same as Fig. 7 but with two orographic ridges: , , and . Non-dimensional simulation details: same as Fig. 7.

Figure 10. Same as Fig. 8 but with two orographic ridges. Non-dimensional simulation details: same as Fig. 7.

Snapshots of the time evolution of the height field are shown in Fig. 2. In case I, we observe two low-amplitude gravity waves propagating to the left and right of the jet core, in agreement with the results of Bouchut et al. (2004) for the standard shallow water theory. Doubling the number of elements reduces the error by a factor 2 (not shown), as expected for a DG0 scheme, verifying numerical convergence. Thus, the model reduces analytically and numerically to the classical rotating shallow water model when the fluid does not exceed the threshold heights and .

For case II, exceedance of modifies the pressure terms, triggering positive buoyancy and leading to a convective updraft. However, no ‘rain’ is produced as is not exceeded. It may be the case that, as , the solution diverges in case II (especially, as since there is no restoring force provided by the downdraft. However, numerical diffusion at the element nodes plays a key role at lowest order where the gradients are steep (i.e. at shocks or significant updrafts), and prevents continuous growth of the convective columns, even in case II. We expect that, in case III, given exceedance and convergence (), ‘rain’ is produced ( contribution) and then slowly removed (i.e, transformed back to precipitable water due to ), providing a downdraft to suppress convection. The strength of the downdraft and consequent suppression of the height field is controlled directly by the parameter. This enhanced suppression is apparent in Figs. 2 and 4, comparing cases II and III: as rain is produced the vertical extent of the updraft in case III is diminished, yet it remains a coherent convective column. Physically, this is due to the feedback of r in (2b(2b) ) and provides justification of the conceptual arguments put forward in Section 2.

The evolution of all four model variables for each case is illustrated in Fig. 3 and detailed further for the fluid depth and rain in Fig. 4. The inertia–gravity waves, indicated by a sharp contour gradient, in the h and u fields in Fig. 3 are clearly apparent as they propagate from the jet core. In cases II and III, the wave fronts exceed the threshold and become convection-coupled waves. The shallow left wave propagates slightly slower than in the standard shallow water case from to 0.5 before leaving the domain. The right wave becomes strongly coupled to the deep convection and slows down. This confirms the wave speed analysis in the ‘Hyperbolicity’ part of Section 2: convection-coupled waves are slower than their ‘dry’ counterparts, and in particular slower in case II than case III. This is also corroborates the numerical simulations of Bouchut et al. (2009).

Multicellular convection (probably the most common form of convection in mid-latitudes) is characterized by repeated development of new cells along the gust front and enables the survival of a larger scale convective system (Markowski and Richardson, 2011c). A basic representation of this is achieved here: the initial convective column subsides around and a new updraft develops in its place with the associated production of rain. The downdraft from the subsiding column instigates a gravity wave that propagates leftward and initiates a region of light convection and rain away from the initial disturbance, another key aspect of atmospheric convection. This is apparent in the top left corner of the Hovmöller plots for h and u in Fig. 3 for cases II and III and the h and r profiles at in Fig. 4.

Figure 5 shows fluid height and positive wind convergence alongside the evolution of r. The production of rain requires both exceedance and convergence, hence we see rain forming in regions where these two processes coincide. It should be noted here that the amount of rain produced and the speed at which it subsequently precipitates is controlled by the parameters and , respectively. Different values would lead to different solutions, not just for hr but all variables, due to the coupling in (2b(2b) ). Moreover, the rate of rain production is directly proportional to the strength of convergence; this explains why there is more rain produced in the deep convection-coupled wave than in the smaller updraft associated with the left-propagating gravity wave.

The Rossby adjustment scenario (Blumen, 1972; Arakawa, 1997) describes how an initial momentum imbalance adjusts to a state of geostrophic balance between the pressure gradient and rotation. Shallow water flow in perfect geostrophic balance satisfies (to leading order with quadratic terms neglected):(20)

In the standard shallow water theory, the geostrophic mean state (i.e. ) is rapidly achieved via the emission of gravity waves (in some cases forming shocks) from the jet core (Bouchut et al., 2004). The shift from large- to convective-scale NWP is in some sense a shift from balanced to unbalanced dynamics. Traditional DA systems developed for large-scale NWP exploit the fact the mid-latitude dynamics at the synoptic scale are close to geostrophic and hydrostatic balance. However, this balance is no longer manifest at smaller scales where rotation no longer dominates and vertical accelerations modulate the flow. Hence, an interesting point here, in the context of convective-scale dynamics and DA, is the disruption of these large-scale balances in the model. By construction of the effective pressure (3a(3a) ), and hence its gradient, a breakdown of the balance (21(21) ) is to be expected in cases II and III, and the numerical results verify this. The top row of Fig. 6 plots the difference (21(21) ) as a function of space and time for the three cases, illustrating where a state close to geostrophic balance is achieved (light shading) and where this balance is broken (deep shading); subsequent rows show profiles of fv and at different times.

In case I, the height field adjusts by emitting shocks from the jet core and quickly approaches the expected balanced state with the Coriolis acceleration fv. Bouchut et al. (2004) note that oscillations may persist for some time in the jet core. Exceedance of the first threshold causes the fluid in that region to rise and instigates deep convection. The gradient of the height field is severely altered and so we see the breakdown of geostrophic balance in the jet (case II: Fig. 6, middle column). The same is true for case III – the height field is qualitatively similar to case II and thus geostrophic balance is not achieved. The leftward propagation of the gravity wave is also manifest here from as a region far from geostrophic balance.

The modRSW model thus exhibits a range of dynamics in which flow is far from geostrophic in the presence of convection whilst remaining ‘classical’ in the shallow water sense in non-convecting and non-precipitating regions. The breakdown of such balance principles is a fundamental feature of convective-scale dynamics and is therefore a desirable feature of the model and its subsequent use in idealized DA experiments.

4.2. Flow over topography

We consider non-rotating (infinite Rossby number) flow over an isolated parabolic ridge defined by:(21)

where is the height of the hill crest, a is the hill width parameter and its location in the domain. Such flow over topography has been extensively researched (see, e.g. Baines, 1998) and is often used as a test case in numerical studies owing to the range of dynamics (dependent on Froude number ), including shocks, and the existence of analytical non-trivial steady-state solutions. Here, we consider supercritical flow with . In this regime, the fluid depth increases over the ridge (as opposed to subcritical flow () in which the depth decreases over the ridge) and a shock wave propagates at a height above the rest depth to the right of the ridge. Such a set-up caters for the present purpose of illustrating the modifications via the hierarchy of model cases as the fluid rises naturally and exceeds the chosen thresholds above the rest height.

The initial conditions are: , , . Since there is no rotation, the transverse velocity v is zero always and the dynamics are purely one-dimensional in space. For standard shallow water flow (case I), the exact steady-state solution is found by solving a third-order equation in h (Houghton and Kasahara, 1968):(22)

Note that although b is a function of x, it is considered a parameter when solving for h. This is obtained by considering the steady-state system (i.e. (1) with and ) and then solving for h conditional on . For modRSW flow, such an analytical equation for the steady-state solution does not exist when (cases II and III). However, it is possible to derive a system of ordinary differential equations (ODEs) in h and r and solve for their steady states for all three cases, which can then be used as a benchmark for the numerical PDE solution for large t for all three cases. The ODE solution for case I matches the analytical solution (23(23) ) (not shown); see Appendix 4 for full details.

Figure 7 shows the evolution of the free surface height and rain r for the three cases. In case I, flow over the ridge reaches the known exact steady-state solution (red-dashed line), thus confirming that correct solutions of the classical shallow water model have not been violated. The ‘convection’ threshold (and later ) is exceeded in two regions: (i) directly over the ridge and (ii) downstream from the ridge where the wave propagates to the right (cases II and III, respectively; Fig. 7), and the long-time numerical PDE steady-state solution (black solid line) for these cases converges to the ODE solution (red-dashed line). As with the previous experiment, the extent of the updraft in case III is slightly reduced owing to the contribution to the hu-momentum equation when r is positive. The extent of this suppression is less than the Rossby adjustment scenario, reflecting the value of in this simulation. We emphasize here that a different choice of (and indeed and ) leads to different dynamics relating to the convection and precipitation. Values chosen here are for illustrative purposes, highlighting the modified the dynamics. When using the model for idealized DA experiments, these parameters can be tuned to yield different configurations as desired.

It is apparent from Fig. 7 that the wave that triggers the downstream updraft becomes a convection-coupled wave and subsequently propagates slower than for the standard SW flow, as was observed in the Rossby adjustment experiment and anticipated by the wave-speed analysis. Rain is produced in and advected with the convective column as it propagates downstream from the ridge and slowly precipitates. Such lee-side enhancement and propagation of deep convection downstream from a ridge is a characteristic phenomenon of orographically induced clouds (Houze, 1993c). Figure 8 plots exceedance and wind convergence alongside r and, as with the Rossby adjustment scenario, illustrates the conditions required for the production of rain. Generating rain both requires and is proportional to positive wind convergence, so we see more rain where this is greater. This corroborates the physical argument put forward in Section 2.2 that rain is produced only when the fluid is rising and the amount of rain is controlled by the strength of the updraft.

Figures 9 and 10 show corresponding results with two orographic ridges. Again, the steady-state solution is achieved in all three cases, whilst the inclusion of a second obstacle for the fluid introduces more complex and nonlinear dynamics with multiple rapidly evolving regions of convection and precipitation.

5. Conclusion

We have presented an idealized fluid model, based on the rotating SWEs and the model of WC14 for cumulus cloud dynamics, intended for use in inexpensive DA experiments at convective scales. Changes to the dynamics are brought about by the exceedance of two threshold heights, akin to (i) the LFC () and (ii) the onset of precipitation (). When the fluid exceeds these heights, the classical shallow water dynamics are altered to include a representation of conditional instability (leading to a convective updraft) and idealized moisture transport with associated downdraft and precipitation effects.

The mathematical modifications to the parent equations described herein, and the physical arguments behind the changes, are strongly motivated by the model of WC14 but improve upon it in two ways. First, the inclusion of a meridional velocity component and Coriolis effects means that dynamics associated with rotating fluids are present in the model. Second and, more importantly, the diffusion terms in WC14 have been removed. The dynamics of WC14 are highly sensitive to these diffusion terms, which are tuned to stabilize the model for a specific set-up and are the dominant controlling factor of the system’s dynamics. As such, the original numerical implementation is not robust to alterations to, e.g. the bottom topography, the gridsize and model parameters, each change requiring ad hoc tuning of the diffusion coefficients and integration time step.

Despite these modifications, the resulting model is shown to be (weakly) hyperbolic, thus constituting a hyperbolic system of non-conservative PDEs. The numerical methodology of Rhebergen et al. (2008) has been developed here for our model and is shown to deal robustly with the threshold nature of the NCPs, whilst well-balancedness is ensured by discretizing the non-zero topography via the method of Audusse et al. (2004).

Classical numerical experiments in shallow water theory, based on the Rossby geostrophic adjustment problem and non-rotating flow over topography, have been reproduced here and used to illustrate the modified dynamics of the model. Crucially, the model reduces exactly to the standard SWEs in non-convecting, non-precipitating regions. This is clear from the model formulation in Equations (2)–(4(4) ), and further confirmed by the numerical model which reproduces known shallow water results in case I. The model also exhibits important aspects of convective-scale dynamics relating to the disruption of large-scale balance principles which are of particular interest from a DA perspective (Bannister, 2010). The Rossby adjustment scenario clearly illustrates the breakdown of geostrophic balance in the presence of convection and precipitation, while the breakdown of hydrostatic balance is implicitly enforced by the modified pressure when the LFC is exceeded. Furthermore, the experiments simulated here have illustrated other features related to convecting and precipitating weather systems, such as the initiation of daughter cells away from the parent cell by gravity wave propagation, and convection downstream from an orographic ridge.

Although based on the model of WC14, the absence of artificial diffusion terms from the governing equations results in a mathematically cleaner formulation with conservation of total mass (‘dry’ plus ‘rain’), and a markedly different dynamical behaviour emerges. With the addition of rotation (and consequent Rossby adjustment dynamics) and analysis of steady-state solutions for flow over topography, we have developed and tested a robust numerical solver and investigated the model’s distinctive dynamics in advance of its use in idealized DA experiments.

DA research using idealized models is primarily carried out in a so-called ‘twin’ experiment setting, whereby the same numerical model is used to generate a ‘nature’ run (which acts as a surrogate truth and is used to generate pseudo-observations) and the forecast. Preliminary results from an idealized forecast assimilation system demonstrate the model’s suitability for conducting inexpensive experiments to evaluate DA schemes in the presence of convection and precipitation (Kent, 2016), and further investigations will presented elsewhere. A basic forecast assimilation framework (briefly comprising scripts for the numerical solver, idealized forecast assimilation routines, plotting and data analysis) is deposited online and is free to access (Kent, 2017). In particular, it is hoped that the numerical solver arising from this study provides a useful tool to the community and facilitates other studies in the field of convective-scale DA research.

Acknowledgements

This work has been undertaken within a CASE award funded by the Engineering and Physical Sciences Research Council and Met Office. We thank Gordon Inverarity for numerous helpful discussions and comments on pre-submission drafts, Michael Würsch for his support at the outset of this work and one anonymous reviewer whose comments helped improve the readability of this manuscript, particularly relating to the construction and explanation of the model.

Additional information

Funding

This work was supported by Engineering and Physical Sciences Research Council [grant number 1305398].

Notes

No potential conflict of interest was reported by the authors.

Appendix 1

Non-dimensionalized modRSW equations

It is useful to work with the non-dimensionalized equations. This parameterizes the problem, yielding non-dimensional parameters that characterize the modelled system and embody its dynamics. The dimensionless coordinates and variables are related to their dimensional counterparts by characteristic scales , and :(23)

Then the dimensionless time coordinate and relevant derivatives are:(A1)

Substituting these into the model equations and defining the non-dimensional effective pressure yields the following dimensionless system (with the hats dropped): (A2) (A3a) (A3b) (A3c) where:(A3d) (A4) (A5) , and the following parameters are introduced:(A6)

The Rossby number, Ro and Froude number, Fr, control the strength of rotation and stratification, respectively, compared to the inertial term .

Appendix 2

Numerics: discretizing the topographic source term

In Audusse et al. (2004), a well-balanced scheme is derived for solving the SWEs with non-flat topography. The two main developments to achieve this are: (i) using reconstructed computational states and to the left and right of an element edge in the numerical flux instead of cell-centred values and ; and (ii) discretizing the topographic source term by considering the leading order balancing requirement for nearly hydrostatic flows. A summary is given here; for full details the reader is referred to Section 2 of Audusse et al. (2004).

In the asymptotic limit for nearly hydrostatic flows the leading order fluid depth h adjusts so as to satisfy the balance of momentum flux and momentum source terms:(A7)

This also ensures that the ‘lake at rest’ property (i.e. the trivial steady-state solution and ) is satisfied. Integrating over element yields an approximation to the topographic source term in the form of a flux:(B1)

The reconstructions for the leading order fluid depth: (B2) (B3a) are truncated to ensure non-negativity of the depth: . Note that a modifed CFL condition imposes a time step restriction also required to ensure non-negativity. Thus, the reconstructed computational states to the left and right of node are:(B3b)

The fluxes in (13) and (14(14) ) are evaluated using these reconstructions and the discretized topographic source term in (17(17) ) is:(B4)

for P defined in (3).

Appendix 3

Numerics: NCP flux derivation

Here, we derive fully the NCP flux (18(18) ) and (19(19) ) for the modRSW model (2). A linear Lipschitz continuous path is defined which connects the left state to the right state: . The model states are approximated elementwise by piecewise constant (DG0) basis functions: . The average of a quantity is denoted by and the jump of a quantity across a node is denoted by . For hyperbolic systems with NCPs of the form (5(5) ), the NCP numerical flux is (Equation (44) in Rhebergen et al., 2008):(B5)

where and are the fastest left- and right-moving signal velocities in the solution of the Riemann problem, and the Riemann ‘star-state’ is given by:(C1)

The integrands involve calculations from the rows of the ‘non-conservative’ matrix in Equation (6(6) ). We define and in terms of the Heaviside function (): (C2) (C3a) and use the following properties of in the derivation:(C3b)

For , the NCPs are zero since the first and third rows of the matrix have zero entries only. In this case the integrals in the flux (C1(C1) ) are zero. Thus, when and the flux is and , respectively. The middle state, , requires further manipulation after substituting (C2(C2) ) into (C1(C1) ):(C4)

The NCP flux reduces to the well-known HLL flux for conservative systems, as alluded to in Section 3.

For , the integrand to be calculated is:(C5)

where we recall that denotes the jump of a quantity across a node, . Integrating over yields the hu component of (19(19) ):(C6)

The expression to be inserted in the flux function (C1(C1) ) then becomes:(C7)

Thus, for , the numerical flux is:(C8)

while for :(C9)

Noting from (C5(C5) ) that the NCP flux for reduces to the HLL flux when the integrand is zero, an expression is obtained from (C1(C1) ) and (C2(C2) ) for the cases when there are non-zero NCPs:(C10)

For , the integrand includes the term, the switch dependent on model variables h, u and topography b. We have that:(C11)

and so the integral to be computed in the flux is:(C12)

To proceed, we set and consider defined by Equation (4(4) ) but with and , so that is a function of :(C13)

It is apparent that depends on the end points of only, and is therefore independent of . If then , and if then . Thus, is equivalent to . It should be noted that this argument is valid for piecewise constant numerical profiles only, i.e. cell averages. A scheme that approximates continuous profiles using means and slopes would require greater consideration.

First, we compute the integral of over [0, 1]:(C14)

where and . When , this integral is trivial:(C15)

For , a change of variable and integration yields:(C16)

Hence,(C17)

Intuitively, this makes sense: when and (i.e. and ), the rain threshold has not been exceeded, meaning no rain is produced, and the above integral is zero.

Proceeding in the same manner, we compute the integral of the product over [0, 1]:(C18)

Again, when , this integral is trivial:(C19)

For , a change of variable and integration yields:(C20)

Hence,(C21)

Equation (C13(C13) ) now reads:(C22)

Thus, for , the numerical flux is:(C23)

while for :(C24)

and finally for :(C25)

This completes the calculations; the NCP flux in vector form is summarized as follows:(C26)

where is the HLL numerical flux:(C27)

and arises due to the NCPs:(C28)

Appendix 4

Numerics: steady-state modRSW solutions

Here, we derive a system of ODEs in h and r to be solved for non-trivial steady states for flow over topography. To facilitate this, we consider a system of equations for h, u, and r: (C29) (D1a) (D1b) where:(D1c)

Steady-state solutions are found by considering time-independent flow (): (D2) The first of these steady-state equations gives immediately a solution of u in terms of h:(D3a)

which is then substituted into the remaining equations, yielding a system of 2 ODEs to solve for h and r. Using (D4(D4) ) and noting that:(D4)

the system in terms of h and r reads: (D5) (D6a) We seek a system of the form , where , prime denotes derivative with respect to x, and , are given from the equations set. If is non-singular (and hence invertible), then we can solve numerically for using, e.g. a simple finite difference scheme.

The system (D6) is expanded as follows: (D6b) (D7a) where if and zero otherwise and the terms in square brackets are components of and :(D7b)

The term (given in (4(4) )) requires further manipulation; re-writing in terms of the Heaviside function we have:(D8)

Thus, the system reads where and is solved using an explicit forward Euler finite difference scheme: . The value at is required to compute the Heaviside of the height gradient in (D9); all other components in are evaluated using values at level j. To start marching through space, note that , so that . Then proceed as usual for .