A two‐fluid single‐column model of the dry, shear‐free, convective boundary layer

A single‐column model of the dry, shear‐free, convective boundary layer is presented in which non‐local transports by coherent structures such as thermals are represented by the partitioning of the fluid into two components, updraught and environment, each with a full set of prognostic dynamical equations. Local eddy diffusive transport and entrainment and detrainment are represented by parametrizations similar to those used in eddy diffusivity mass flux schemes. The inclusion of vertical diffusion of the vertical velocity is shown to be important for suppressing an instability inherent in the governing equations. A semi‐implicit semi‐Lagrangian numerical solution method is presented and shown to be stable for large acoustic and diffusive Courant numbers, though it becomes unstable for large advective Courant numbers. The solutions are able to capture key physical features of the dry convective boundary layer. Some of the numerical challenges posed by sharp features in the solution are discussed, and areas where the model could be improved are highlighted.


INTRODUCTION
Motivated by current challenges in representing cumulus convection in weather and climate models (e.g. Holloway et al., 2014;Gross et al., 2018),  recently proposed an approach in which a full set of prognostic equations for density, momentum, potential temperature, and moisture is used for the resolved-scale average dynamics of convective updaughts, as well as for their environment. (Downdraughts and multiple classes of convective updraught can also be included with their own prognostic equations.) The resulting governing equations resemble those used in modelling multi-phase flows (e.g. Drew, 1983;Städtke, 2007), so we call this approach the "multi-fluid" approach. An attractive and novel feature of the multi-fluid approach is that, in a numerical model, it allows the resolved-scale average dynamics and transport of the convective updraughts to be handled by the dynamical core. Processes such as entrainment, however, must still be parametrized. In order for the multi-fluid approach to be practically useful, suitable numerical solution methods must be developed. Also, suitable schemes must be formulated for the subfilter-scale fluxes and for the sources and sinks of each type of fluid that represent processes such as entrainment and detrainment. This paper presents a first step in this direction in the form of a two-fluid single-column model of the dry convective boundary layer.
The multi-fluid governing equations can be derived systematically by conditional filtering. The idea is analogous to the filtering procedure used to derive the equations of large-eddy simulation (LES) but, in addition, it makes use of a set of quasi-Lagrangian fluid labels to pick out different regions of the fluid, for example convective updraughts and their environment. At any point in the fluid exactly one of the labels takes the value 1 while the rest are 0. Multiplying one of the dynamical equations (e.g. the mass or momentum equation) by one of these labels before applying a spatial filter gives the corresponding equation for the filter-scale average of the fluid component picked out by that label. The resulting equations include terms representing the effects of subfilter-scale fluxes as well as terms representing the relabelling of fluid parcels.  provide a detailed derivation of the conditionally filtered equations and discuss the relationship between relabelling and entrainment/detrainment.  have documented the conservation properties of the multi-fluid equations and shown that they possess physically reasonable normal modes, thus giving some confidence that the multi-fluid equations do provide a physically sound basis on which to model convective flows.
Under suitable approximations the multi-fluid governing equations reduce to those for a typical mass flux convection scheme . They also include subfilter-scale terms that could account for local turbulent transports such as those commonly modelled by a downgradient eddy diffusion in the boundary layer. Thus, the multi-fluid equations have clear conceptual links to existing widely used approaches. At the same time, one of the attractions of the multi-fluid approach is the possibility that, in a numerical model, the dynamics of convective updraughts and downdraughts could be handled by the dynamical core, leaving only smaller-scale processes such as entrainment and detrainment, local turbulent fluxes, and microphysics to be parametrized. Such an approach would be a significant departure from current practice, and would shift the traditional, but artificial, distinction in weather and climate modelling between "physics" and "dynamics." By allowing the dynamical core to handle the resolved-scale dynamics of convection, the multi-fluid approach has the potential to overcome several limitations of conventional convection parametrization schemes. For example, a dynamical memory of the state of convection (e.g. Plant and Craig, 2008), including cold pools (e.g. Grandpeix and Lafore, 2010), would be included in a natural way, convective systems could propagate to neighbouring grid columns without having to switch off and reform (e.g. Grandpeix and Lafore, 2010), and compensating subsidence would not be parametrized in the convecting grid column (e.g. Krueger, 2001;Kuell et al., 2007) but could occur where required by the dynamical equations.
In order to realize this potential of the multi-fluid approach, it is necessary to develop suitable forms of the parametrized terms, along with suitable numerical methods for solving the governing equations, and to show that they capture the relevant physical processes. The aim of this paper is to document an initial step in that development process for a simple, but non-trivial, single-column problem and so demonstrate a proof of concept for the multi-fluid approach.
For our purpose it is appropriate to start with the simplest relevant problem: a single-column two-fluid model of the dry convective boundary layer. In the convective boundary layer, there is a significant contribution to the vertical potential temperature transport from large coherent eddies or thermals, which is up the mean gradient in the upper half of the boundary layer. This transport and the resulting mean potential temperature structure cannot be accurately modelled only by a downgradient eddy diffusion (e.g. Holtslag and Boville, 1993); the ability of the two-fluid model to capture this non-local and upgradient transport will be a valuable test of the approach. Also, beginning with the dry convective boundary-layer case allows us to take advantage of the close conceptual similarity between the two-fluid model and the eddy diffusivity mass flux (EDMF) model of Soares et al. (2004) and Siebesma et al. (2007). In both models local turbulent eddy fluxes are represented in terms of an eddy diffusivity. In the EDMF model, non-local fluxes are represented by a steady entraining plume mass flux scheme, while in the two-fluid model they are represented by the second, convecting, fluid. Thus, we are able to adopt or adapt the parametrizations of a number of processes from the EDMF approach for use in the two-fluid approach. The conceptual similarity of the multi-fluid and EDMF approaches is further emphasized by the recent extension of the EDMF approach by Tan et al. (2018) to include prognostic equations for updraught properties; their model is then essentially a two-fluid anelastic model.
For the moist convective boundary layer, EDMF schemes have been extended to allow multiple categories of updraught Sušelj et al., 2012). In the dry case, however, a single updraught appears to be sufficient. For this reason we limit ourselves here to two fluid partitions in the multi-fluid scheme. By restricting attention to a dry single-column model, the dynamical and numerical issues are simpler than in three dimensions and we avoid the complications associated with moist processes. Nevertheless, there are still some non-trivial issues to address and valuable lessons to be learned, as discussed below. The governing equations for the dry single-column two-fluid model are presented and discussed in Section 2.
Note that we should not expect the two-fluid model to outperform the EDMF model on this test problem. The eddy turnover time (a few minutes) is much shorter than the time-scale of evolution of the boundary-layer mean fields (hours), so the dynamical memory of the convective updraughts in the two-fluid model will provide no advantage. Also, in this single-column case, compensating subsidence must necessarily occur in the same grid column, and there can be no propagation of convection to neighbouring grid columns. Thus, if the two-fluid model can produce similar results to EDMF, then such an outcome may be regarded as an initial proof of concept and demonstration of feasibility of the multi-fluid approach.
The EDMF scheme of Siebesma et al. (2007) is formulated in such a way that the volume fraction or area fraction of updraughts is never used explicitly. In the multi-fluid approach, the volume fraction of updraughts is part of the solution, and suitable values of updraught volume fraction and vertical velocity are crucial for obtaining adequate mass fluxes and convective transports. Section 3 discusses the determination of a suitable profile of updraught volume fraction, as well as the calibration of entrained and detrained values of vertical velocity and potential temperature.
Since the eventual intended application of the multi-fluid approach is in weather and climate models, we use numerical methods for the resolved-scale dynamics based on the ENDGame dynamical core (Wood et al., 2014) used operationally in the Met Office Unified Model. It is a semi-implicit semi-Lagrangian scheme, which should provide a stable treatment of fast waves and advection with long time steps. Since vertical eddy diffusion and entrainment and detrainment are all fast processes on the scales of interest, the time integration scheme is extended to give an implicit treatment of these terms fully coupled to the dynamics. The details are provided in Section 4.
Some sample results are presented in Section 5, showing that the model captures key physical features of the dry convective boundary layer. A valuable test of the well-posedness of a model or parametrization scheme and its numerical solution is that its solution should converge with increasing resolution. Section 5 also discusses the convergence of the two-fluid boundary-layer model. Section 6 summarizes the conclusions and discusses some areas where further development could improve the model.

GOVERNING EQUATIONS
The governing equations are based on the conditionally filtered compressible Euler equations, as given by . They are restricted to a single column, so that the horizontal velocity vanishes and the solution is independent of the horizonal coordinates. There are n = 2 fluid components, i = 2 representing the convecting or updraught fluid component and i = 1 representing the rest of the fluid. The governing equations then reduce to the following.
Here, i is the filter-scale volume fraction of fluid component i, and Equation 1 represents the fact that the volume fractions of the different fluid components must sum to unity. Equation 2 expresses conservation of mass, with i the filter-scale density of fluid component i, and w i the filter-scale vertical velocity of fluid component i. The prognostic thermodynamic variable is chosen to be the potential temperature rather than the entropy to retain a similarity to ENDGame. Neglecting diabatic heating, the filter-scale potential temperature i satisfies Equation 3, where D i ∕Dt = ∕ t + w i ∕ z means the 'material' rate of change following fluid i, and F i SF is the subfilter-scale flux of within fluid component i. The vertical momentum equation is Equation 4. In the pressure gradient term Π = (p∕p 0 ) is the Exner pressure, where p 0 = 10 5 Pa is a constant reference pressure and = R∕c p , with R the gas constant for dry air and c p the specific heat capacity at constant pressure. The gravitational acceleration is Φ∕ z, where Φ is the geopotential, and F w i SF is the subfilter-scale flux of w within fluid component i. The terms d i allow for the fact that the net pressure gradient force on fluid i is not exactly c p i Π∕ z but also includes mutual pressure forces (typically a "drag") between fluids 1 and 2 as well as other subfilter-scale pressure variations. Conservation of momentum requires ∑ i d i = 0. In order to write Equation 4, a perfect gas equation of state Equation 5 has been assumed, and, as is usually done, some subfilter-scale contributions to Equation 5 have been neglected. The thermodynamic equation and momentum equation are written in advective form in anticipation of a semi-Lagrangian discretization, again as in ENDGame.
In order to represent entrainment and detrainment, fluid parcels are relabelled. This relabelling appears in the conditionally filtered equations as the term  ij , the local rate at which fluid component j is relabelled as fluid component i; thus  21 is the entrainment rate and  12 is the detrainment rate, both expressed as mass per unit volume per unit time. The fluid that is relabelled carries its potential temperaturêi j and its vertical velocityŵ ij with it, thus affecting the budgets of potential temperature and vertical momentum for both fluids.  give a full derivation and further discussion; also see Section 2.3.
Note that the two fluid components share the same Exner pressure field Π. This formulation has several important consequences : the number of equations matches the number of unknowns; the pressure may be interpreted as a Lagrange multiplier associated with the constraint Equation 1; and inter-fluid acoustic modes are filtered so that the linearized equations possess physically reasonable normal modes. The Exner pressure Π may be diagnosed from the prognostic variables i i and i by taking i times Equation 5 and summing over i: after which i is obtained by back-substitution in Equation 5 and We have chosen to work with a fully compressible governing equation set, rather than make a Boussinesq or anelastic approximation, because modern high-resolution weather prediction models are increasingly based on fully compressible governing equations. Consequently, the numerical methods used must stably handle acoustic waves even at very high acoustic wave Courant numbers (Section 4).
The terms on the right-hand sides of Equations 2-4 must be parametrized. The parametrizations used here, and discussed in the following subsections, are mostly based on those in the literature on EDMF schemes. However, it is highly desirable that a model should be based on a well-posed set of differential equations (Cullen et al., 2001;Arakawa and Wu, 2013), distinct from the numerical methods used to solve those equations, so that it makes sense to talk about "the exact solution" of the differential equations, and convergence of a numerical solution to the exact solution (Williamson, 2008). Guided by this principle, some of the parametrized terms have been modified; see also the discussion in Section 6.

Boundary-layer depth and vertical velocity scale
Following Soares et al. (2004), Siebesma et al. (2007), Neggers (2009) and Neggers et al. (2009), several of the parametrized terms are specified in terms of the boundary-layer depth z * and an associated convective vertical velocity scale w * . A method is therefore required to diagnose the boundary-layer depth. Following these authors, we take z * to be the height at which the updraught velocity goes to zero. The convective vertical velocity scale is given by w * = (gz * Q * ∕ 00 ) 1∕3 , where Q * is the surface potential temperature flux and 00 = 300 K is a characteristic surface potential temperature. In the present work Q * is a specified constant.

Vertical diffusion coefficient
As in the EDMF approach, the local subfilter-scale flux of potential temperature F i SF is represented by a downgradient eddy diffusive flux. Here the eddy diffusive fluxes are applied separately in both the environment and updraught fluid components, rather than to the mean fluid potential temperature. For simplicity, and to aid comparison with previous work, the same diffusion coefficent is used in the updraught and the environment. Soares et al. (2004), Siebesma et al. (2007), Neggers (2009) andNeggers et al. (2009) use a specified vertical profile of eddy diffusivity, based on Holtslag (1998), where u * is the friction velocity, k = 0.4 is von Kármán's constant, andz = min(z∕z * , 1) is the height normalized by the boundary-layer depth. The upper bound of 1 onz is a convenient way to impose zero diffusion above the boundary layer. Following Siebesma et al. (2007), we consider the shear-free case and so take the background horizontal velocity to be zero and set u * = 0.
To ensure that the diffusivity remains non-zero and the potential temperature gradient remains finite at the surface, we follow the standard approach of including a roughness length z 0 , leading to a modified definition ofz: The numerical experiments discussed below use z 0 = 0.1 m. The specified potential temperature flux Q * is imposed as the bottom boundary condition for the potential temperature diffusion in both fluids.
Our initial implementation included vertical eddy diffusion only on the potential temperature fields 1 and 2 . However, it was found that the numerical solutions often became unstable, even when additional numerical damping mechanisms were included. Subsequent analysis (Appendix A) revealed that the two-fluid equations with unequal basic state flow in the two fluids are subject to a Kelvin-Helmholtz-like instability; this is an inherent instability of the continuous governing equations, not a numerical instability. The stability analysis suggests that the inclusion of a sufficiently strong vertical diffusion on w, parametrizing the effects of the F w i SF terms, would stabilize the flow, at least for small enough vertical wavelengths. Therefore, we include vertical diffusion of w 1 and w 2 as well as 1 and 2 . We use the same profile of the vertical diffusion coefficient (Equation 7) for w as for . The vertical diffusion of w uses the bottom boundary condition w i = 0 in both fluids.

Entrainment and detrainment
Mass flux schemes are known to be particlularly sensitive to the parametrization of entrainment and detrainment (e.g. Romps, 2016, and references therein). EDMF schemes for the convective boundary layer have used a variety of forms for the fractional entrainment per unit depth , including a constant value (Angevine, 2005), a specified function of normalized depth (Soares et al., 2004;Siebesma et al., 2007), a rate inversely proportional to updraught speed , and a rate inversely proportional to a diagnosed turbulence length-scale (Witek et al., 2011b). Many EDMF schemes prescribe the vertical profile of mass flux or of updraught volume fraction within the boundary layer or subcloud layer. When this is done, the vertical profile of detrainment is implied by the updraught mass budget and so does not need to be explicitly parametrized. The scheme of Angevine (2005) and Angevine et al. (2010) is one example that does explicitly parametrize detrainment; the fractional detrainment rate per unit depth is a simple profile of z which, in the absence of condensation, peaks at the boundary-layer top.
In the two-fluid model we do not prescribe either the mass flux profile or the updraught volume fraction profile; they both emerge as part of the solution. Therefore we are obliged to parametrize both entrainment and detrainment. The governing Equations 2-4 are expressed in terms of the absolute entrainment rate  21 and detrainment rate  12 per unit time rather than the fractional entrainment and detrainment rate per unit depth. The  ij are related to and by where is the updraught mass flux. 1 Entrainment and detrainment rates are notoriously difficult to infer directly from high-resolution LES, partly because they depend critically on how air is labelled as belonging to updraught or environment (e.g. Romps, 2010;Dawe and Austin, 2013;Yeo and Romps, 2013). Entrainment rates may also be estimated indirectly from LES using a bulk plume model where is some materially conserved scalar. However, a number of approximations are involved in Equation 12, including statistical steadiness of the flow, small updraught volume fraction, and the assumption that the value of in the entrained air is equal to the horizontal mean . Romps (2010) argues that this method systematically underestimates the actual entrainment.
In the present work, the entrainment and detrainment rates are expressed in terms of corresponding time-scales ij , wherêi j is the density of the relabelled air. This form is convenient for the approximate linearization of the equations needed for their semi-implicit time integration (Section 4). For simplicity we takêi j = j . Combining Equation 13 with 9 and 11 and noting that 1 ≈ 2 and 1 ≈ 1 gives allowing previously published formulations for to be adapted for the two-fluid model. Here we use a form based on that used by Cheinet (2003): where c = 0.4 is a dimensionless constant, so that In practice we have found it helpful for capturing an initial spin-up from rest to define the inverse time-scale 1∕ * by w 2 max ∕z * , where w 2 max is the maximum speed in the updraught profile, rather than the more obvious w * ∕z * ; this 1 Note that here M is defined in terms of the absolute updraught velocity w 2 rather than a perturbation from the mean ascent, as is often done.
has the effect of reducing the entrainment rate during the first few minutes, but makes little difference thereafter. There are fewer studies or data on which to base a parametrization of detrainment for the dry convective boundary layer. A crucial aspect of any such parametrization is that it must provide enough detrainment near the top of the boundary layer to detrain all of the air arriving in updraughts and so prevent the build-up of a large volume fraction of fluid 2 at the boundary-layer top. Siebesma et al. (2007) use a prescribed mass flux profile and suggesting Based on LES results, Rio et al. (2010) propose a parametrization in which air is detrained when it becomes negatively buoyant, with the detrainment rate increasing strongly as the updraught slows: where b is the updraught buoyancy. In this formulation the detrainment rate peaks very strongly at the top of the boundary layer.
In practice we have found that a combination of Equations 19 and 20, slightly modified and with tuned coefficients, works well. The contribution given by Equation 19 helps to maintain a suitable mass flux profile in the midboundary layer, while the contribution given by Equation 20 prevents a build-up of fluid 2 at the boundary-layer top. Expressing the resulting detrainment profile in terms of a detrainment time-scale gives wherew = max(0.01, w 2 ∕w * ) and c = 0.7 is a dimensionless constant, with buoyancy defined as b = g( 2 − 1 )∕ 00 . The factor 1.001 in the denominator of the first term and the bounding ofw away from zero are needed to avoid a singularity in the detrainment rate at z = z * . It is common in the literature to use different entrainment rates for different variables. For example, if is taken to be the entrainment rate for potential temperature, then Angevine (2005) uses 2 as the entrainment rate for vertical velocity while Siebesma et al. (2007) and Neggers et al. (2009) use 0.5 as the entrainment rate for vertical velocity. This practice is defensible, and even natural, if the entrainment rate is defined by an equation like 12. However, the multi-fluid equations are based on the intuitive idea that the entrainment rate is a rate of relabelling of fluid mass. In this context it is not permitted, and does not make sense, to have different entrainment rates for different variables. Instead, the multi-fluid framework permits the properties of the entrained air̂2 1 and w 21 to differ from those of the environmental air 1 and w 1 ≈ 0. Similar comments apply to detrainment. Here we take the properties of the entrained air to be linear combinations of the mean properties of the updraught and environment air: where i = 1, 2 and j = 3 − i. For entrainment of w, we take b w 2 = 0.5 and for detrainment b w 1 = 1. These choices are justified by studies of the vertical momentum budgets of updraughts (e.g. Sherwood et al., 2013;Romps and Charn, 2015), which suggest that the entrainment of vertical momentum is smaller than would be predicted based on entrainment of stationary environmental air. Use of Lagrangian particles to estimate entrainment and detrainment in our own LES results (Section 3) also supports the choice b w 2 ≈ 0.5 and b w 1 ≈ 1. Finally, with this choice, entrainment has approximately the same effect on the vertical momentum budget as the entrainment of stationary environmental air at a rate 0.5 , as used by Siebesma et al. (2007) and based on their LES results and a bulk plume model. For we set b 2 = 1 and b 1 = 1. Our LES results (Section 3) suggest values of b i closer to 0.5. However, this choice led to excessively large updraught potential temperatures and updraught advective potential temperature fluxes; the model responded by moving the potential temperature minimum downward to an unrealistic height z∕z * ≈ 0.2 in order for the diffusive potential temperature flux to compensate for the excessive advective flux.

Updraught-base properties
In the mass flux approach, some scheme is needed to specify the properties at the base of the updraught. In the literature on EDMF schemes, the updraught base value u of some materially conserved variable is typically specified as where an overbar indicates the grid cell horizontal mean and w ′ ′ | s is the surface flux of . There is considerable variation in the values used for the turbulent velocity scale W, which may be set to w * (Angevine, 2005;Neggers et al., 2009;Sušelj et al., 2012), to an estimate of the subgrid standard deviation of w at the lowest model level (Siebesma et al., 2007), or to the square root of the subgrid turbulent kinetic energy (TKE; Soares et al., 2004;Witek et al., 2011b). The values of the dimensionless scaling factor range from 0 to 10. There is also considerable variation in the value used for the updraught base vertical velocity, which may be set to some constant multiple of w * (Angevine, 2005;Sušelj et al., 2012), to some constant multiple of the subgrid standard deviation of w at the lowest model level ), or to zero (Witek et al., 2011b;Tan et al., 2018).
In the two-fluid model, we must specify the volume fraction or mass fraction of fluid 2 at the lower boundary as well as the updraught properties. In practice the mass fraction must be specified in the lowest model layer, where it is set to This choice is based on values used in EDMF schemes and estimated from LES (Section 3), followed by tuning to improve the agreement between the two-fluid model and LES results. Following Equation 24 with parameters similar to those used by Siebesma et al. (2007), we impose with̃= 1.5 (̃≈ ∕ 1 ) and W given by s w , the subgrid standard deviation of w, (Holtslag and Moeng, 1991) evaluated at the lowest model layer. Once again the inclusion of the roughness length inz avoids a singularity at z = 0 when u * = 0, as assumed here. These updraught-base properties are imposed by modifying the entrainment and detrainment terms for 1 1 , 2 2 , 1 , and 2 in the lowest model layer so as to force these quantities to satisfy Equations 25 and 26 while conserving the mean mass = 1 1 + 2 2 and potential temperature content * = 1 1 1 + 2 2 2 , where the notation * indicates a mass-weighted mean. As noted by , conserving the total potential temperature content has the advantage of also conserving the total internal energy and preserving the pressure (Equation 5). No adjustment is made to w 1 or w 2 , since both are zero at the bottom boundary. Imposing Equations 25 and 26 at the updraught base can be interpreted as a sorting of the subgrid PDF to assign the most buoyant fraction of the fluid to fluid 2, anticipating that the dynamics will then cause it to ascend in an updraught.

Pressure drag
Taking a horizontal mean (again indicated by an overbar) of the flux form of the vertical momentum equation, before conditional filtering, and neglecting transience of the mean state, leads to (Schumann and Moeng, 1991), where here p ′ is the departure of the pressure from a profile in hydrostatic balance with . Schumann and Moeng (1991) conclude that the non-hydrostatic pressure gradient accelerates rising thermals in the lower half of the convective boundary layer and decelerates them in the upper half. The LES results of Siebesma et al. (2007) show that the square of the updraught vertical velocity is roughly proportional to the vertical velocity variance, and they use this, in conjunction with Equation 28, to parametrize the effect of the non-hydrostatic pressure on the updraught vertical momentum.
Here we do not include such a contribution to the updraught momentum equation for the following reason.  in flux form, summing over i using Equation 32 below, and neglecting transience gives where Π ′ is the departure of Π from a profile in hydrostatic balance with * . Equation 29 is the two-fluid model analogue of Equation 28, and shows that the effects of the horizontal mean non-hydrostatic pressure are explicitly resolved in the two-fluid model and so do not need to be parametrized. We have confirmed that the balance expressed by Equation 29 holds in our numerical solutions, with the F w i SF terms providing a relatively small contribution.
LES experiments (Romps and Charn, 2015) have shown that a pressure perturbation dipole typically exists across individual buoyant thermals that opposes their upward acceleration. Such a drag is not captured explicitly by Equation 29 and so must be parametrized via the d i terms in Equation 4. Romps and Charn show that their data are consistent with a drag given approximately by where the thermals are approximated as spherical with radius R. If the thermal radius is assumed proportional to z * within the dry convective boundary layer, then these results suggest that the drag can be parametrized by taking and, for momentum conservation, The experiments discussed below use = 1.

UPDRAUGHT, ENTRAINMENT AND DETRAINMENT PROPERTIES
An important factor in the performance of the two-fluid model is the simulated profile of 2 , the volume fraction of convecting fluid. It has a strong influence on the non-local contribution to the transport and hence on the resulting boundary-layer potential temperature profile as well as the boundary-layer-top entrainment flux.
Observations of the convective boundary layer (e.g. Young, 1988;Miao et al., 2006, and references therein) suggest that the fractional coverage of large coherent updraughts or thermals is of the order 10-20% or more, though the diagnosis is very sensitive to how updraughts are defined. Siebesma et al. (2007) used updraught fractions of 0.01, 0.03 and 0.05 to diagnose updraught properties from their LESs to calibrate their EDMF scheme. However, the updraught fraction implied by the prescribed mass flux profile used in their EDMF scheme is significantly larger. Several EDMF schemes (Soares et al., 2004;Witek et al., 2011b;2011a) use a constant updraught volume fraction in the boundary layer of order 0.1.
Updraught volume fraction and other properties can be diagnosed from LES of the convective boundary layer, but again the calculation is sensitive to how updraughts are defined. Couvreux et al. (2010) propose a sampling criterion that captures well the updraughts from the surface to the top of the convective boundary layer. A passive but decaying tracer is released at the surface; throughout the boundary layer, air with positive vertical velocity and with tracer mixing ratio greater than the mean plus one standard deviation at that height is labelled as updraught air. They found an updraught fractional coverage decreasing slightly with height from 0.2 to 0.13.
Given this uncertainty in the fractional coverage of updraughts, we also estimated values from our own incompressible Boussinesq LES. The case studied is case 1 of Siebesma et al. (2007). It has an initial potential temperature profile (z) = (297.2 − 1.95 × 10 −3 z) K (with z in metres), and a constant surface potential temperature flux Q * = 0.06 Km/s is applied so that a convective boundary layer grows gradually into the stable background. The domain size is 4.8 km × 4.8 km horizontally and 3 km deep, and the horizontal and vertical grid spacings are 25 and 10 m, respectively. Updraught fractions (and other properties) were estimated following the method of Couvreux et al. (2010). The profile of updraught fraction estimated in this way, averaged over the sixth hour of simulation, is plotted in Figure 1a, and is in the range suggested by earlier studies. The updraught vertical velocity and mass flux are also shown for comparison with the results of the two-fluid model below. Figure 2 shows updraught and downdraught values of and w, along with values for entrained and detrained air. The entrained and detrained values are estimated as the average values on Lagrangian particles that change their label from environment to updraught or vice versa. The estimated entrained valueŵ 21 is roughly mid-way between w 1 and w 2 , while the detrained valuesŵ 12 are close to w 2 . This result gives some support for the values of b w i chosen in Section 2.3. For potential temperature both the entrained valueŝ2 1 and detrained valueŝ1 2 are roughly mid-way between 1 and 2 . As noted in Section 2.3, values of b i based on these estimates do not give good results in the two-fluid model. A likely explanation for this discrepancy is as follows. The Lagrangian particle method is known to suffer from excessive switching of particle labels between updraught and environment (e.g. Yeo and Romps, 2013) leading to larger estimates for entrainment and detrainment rates than other methods. The entrainment and detrainment rates diagnosed in this way in our LES are indeed significantly larger, especially in the middle of the boundary layer, than those parametrized in the two-fluid model. These large label-switching rates will also bias the values of̂diagnosed from the LES towards less extreme values.
Further details and results from this work will be discussed at length elsewhere.

NUMERICAL METHODS
The governing equations are solved using a semi-implicit, semi-Lagrangian solution method. This choice was intended to provide a stable treatment of fast waves such as acoustic waves, and to permit large vertical advective Courant numbers so that the time step would not be restricted by the vertical velocity in strong updraughts. Also, the scheme used resembles a one-dimensional version of the ENDGame scheme (Wood et al., 2014) used operationally at the Met Office, which should facilitate the adoption of the multi-fluid approach if it proves successful. The staggering of variables on the grid is a generalization of the Charney and Phillips (1953) grid, with i and w i staggered relative to i i and Π.
The semi-implicit, semi-Lagrangian discretization of Equations 2-4 may be written concisely in the form Here, R i comprises the terms in the mass continuity equation for fluid i that depend only on the (known) variables at time step n, while L i comprises the terms that depend only on the (unknown) variables at time step n + 1. Analogous definitions apply for the and w equations. The operator  represents transport by standard semi-Lagrangian advection. In our initial experiments we found that semi-Lagrangian interpolation of i i led to very large conservation errors ( i i is much less smooth than i itself); therefore instead we use a one-dimensional version of the conservative large-time-step advection scheme SLICE (Zerroukat et al., 2009), indicated by the operator  . Like a standard semi-Lagrangian scheme, SLICE works by computing departure points for trajectories that arrive at grid points, but it achieves conservation by remapping cell integrals of the advected quantity rather than interpolating point values.
In more detail, the terms in Equations 33-35 are given by These equations apply for i = 1, 2, with j = 3−i. Superscripts n and n + 1 indicate time steps, and the and = 1 − coefficients are off-centring parameters for the semi-implicit time scheme. The numerical experiments discussed below use = 0.65. All vertical derivatives are evaluated using centred finite differences.
Note that vertical diffusion is a fast process, so for stability we treat it with a backward time step. To improve the balance between diffusion and other processes, the diffusion contribution is fully coupled to the semi-implicit semi-Lagrangian stepping of other terms, not treated in a time-split way. Entrainment and detrainment are also fast processes. They are treated in a semi-implicit semi-Lagrangian way with the same off-centring weights as the dynamical terms in the equations; in this way entrainment and detrainment are sampled at the beginning and end of each air parcel trajectory.
In order to solve Equations 33-35, a quasi-Newton method is used. This is essentially a Newton method (e.g. Press et al., 1988) in which the Jacobian is approximated, and in which the linear system that results at each iteration is not solved exactly. After some number l of Newton iterations, Equations 33-35 will not be satisfied exactly, but there will be some residuals given by where the superscripts on the L terms and  and  operators indicate that they are evaluated using the iteration l estimates for the fields at step n + 1. The Newton method seeks increments (indicated by a prime) to the estimated step n + 1 values, designed to reduce the residuals, by solving the linear system Here, where the  i are tridiagonal operators that include contributions from the diffusion of w i , and the  i are weighted finite-difference gradient operators. Using Equations 46, 48 and 49 to eliminate ′ i , ′ i and ′ i from Equation 45 gives an equation of the form where the  i are modified finite-difference divergence operators. The right-hand side terms  w i and  Π are functions of the residuals Q i , Q i , Q w i . In the absence of the diffusion terms, the tridiagonal  i operator in Equation 50 reduces to multiplication by a coefficent; in that case the w ′ i can be eliminated from Equations 50, 51 to leave an equation for the single unknown Π ′ that has the same structure as a one-dimensional version of a standard semi-implicit Helmholtz problem (e.g. Wood et al., 2014). In discrete form this takes the form of a tridiagonal system, for which Gaussian elimination gives an efficient numerical solution. In our case, however, the w-diffusion terms must be included in the linearization. Then no further analytical elimination is possible, and the coupled system Equations 50, 51 must be solved numerically. With the unknowns w ′ i and Π ′ ordered according to their height, Equations 50 and 51 form a heptadiagonal system, which again can be solved directly by Gaussian elimination.
Once w ′ i and Π ′ are found, the other increments are found by back-substitution. All of the estimates for the step n + 1 prognostic variables are then updated: for a generic variable , Provided the Newton iterations converge, the updated variables approach the solution of the nonlinear system Equations 33-35 and 6. Note that all variables or coefficients labelled as being at time level n + 1, including K,  and d, are updated at each iteration; there are no "lagged" or "frozen" coefficients. Four Newton iterations are used for all the results shown below. Some other relevant numerical details are as follows. Seeking quasi-steady solutions of the continuous governing equations for z 0 ≪ z ≪ z * yields the well-known convective boundary-layer results that w i should scale like z 1∕3 while i and 2 − 1 should scale like z −1∕3 . To improve the accuracy of the advection of these fields, the semi-Lagrangian interpolation of R w i is done in a stretched coordinate z 1 = z 1∕3 while the semi-Lagrangian interpolation of R i is done in a stretched coordinate z 2 = z −1∕3 0 − (z + z 0 ) −1∕3 . Cubic Lagrange interpolation is used, dropping to linear interpolation for departure points in the lowest and highest model layers. The semi-Lagrangian departure point calculation is also done in terms of a stretched coordinate z 3 = z 2∕3 , and vertical interpolation of w 2 to compute the entrainment rate is done in the z 1 coordinate.
Conservation of the mass integral of is crucial for the correct growth of the boundary layer. The diffusion of is discretized in space in a way that conserves the integral. However, the semi-Lagrangian advection is not conservative, so a simple post hoc conservation fixer is applied.
In order to ensure a smooth time evolution of z * and the quantities that depend on it, z * must be allowed to take values in between w-levels. Similar to the method used by Siebesma et al. (2007), we take z * to be the height at which a linear extrapolation of w 2 2 goes to zero.

Physical aspects
Figures 1d-f, 3, 4 and 5 show some example results from the two-fluid single-column model. The test case is again based on case 1 of Siebesma et al. (2007), with the same initial potential temperature profile and surface potential temperature flux. The governing equations are now fully compressible rather than Boussinesq or anelastic. The domain is 2400 m deep discretized on a uniform 20 m grid with a 6 s time step. Figure 3 shows a sequence of vertical profiles of the mean potential temperature * at the initial time and after 1.5, 3.5, 5.5 and 7.5 hr from the single-column model and, for comparison, from the LES. The boundary-layer depth grows at a realistic rate. There is a superadiabatic layer near the surface, a strong inversion at the boundary-layer top, and a well-mixed layer in between with a potential temperature minimum near the middle of the boundary layer. The single-column results agree well with those from EDMF schemes (e.g. Siebesma et al., 2007). They are also similar to those from LES, except that the boundary-layer top inversion is unrealistically sharp; this excessively sharp inversion is common in EDMF schemes too (e.g. Siebesma et al., 2007;Witek et al., 2011b) and suggests that improvements to parametrized terms near the boundary-layer top may be needed. Figure 4 shows the vertical profiles of * and 2 after 8 hr, relative to the minimum value of * and scaled using * ≡ Q * ∕w * . Again the results are very similar to those of EDMF schemes (e.g. Siebesma et al., 2007). The updraught is positively buoyant up to z ≈ 0.9z * ; above this the updraught becomes negatively buoyant, contributing to the boundary-layer top entrainment flux. The updraught potential temperature profile is almost identical at 20 m resolution and 5 m resolution. However, the height of the minimum in the mean potential temperature profile is lower at finer resolution.
Figure 1d-f shows scaled vertical profiles of 2 , w 2 , and 2 w 2 after 8 hr. The updraught fraction 2 is fairly uniform at around 0.15, decreasing quickly to zero at the boundary-layer top. The updraught vertical velocity and mass flux profiles are very similar to the LES results shown in panels (b)-(c). The mass flux profile is also very close to the prescribed profile used by Siebesma et al. (2007) (dashed curve in (f)). Figure 5 shows two contributions to the potential temperature flux: the local eddy diffusive contribution ∑ i i i K i i ∕ z, and a resolved advective contribution , along with their sum. The sum decreases linearly with height, consistent with a uniform heating rate, from Q * at the surface to zero at a height of about 0.9z * . In the lowest part of the boundary layer, the resolved transport contribution increases with finer vertical resolution, and this increase is compensated by a reduced eddy diffusive contribution. Above 0.9z * the net transport reaches a minimum -(the negative of) the entrainment flux -of about −0.1Q * before returning to zero at z = z * . This entrainment flux is consistent with values seen in LES. Note, however, that its magnitude fluctuates considerably as the boundary-layer top moves between model levels (Figure 5c) -another There is one further contribution to the potential temperature transport not shown in Figure 5, namely that due to the mean ascent: . In Boussinesq and anelastic models this contribution vanishes identically. However, for the fully compressible system used here, heating of the boundary layer causes expansion leading to a mean ascent that peaks at the boundary-layer top. The resulting potential temperature transport has a peak value of about 0.5Q * at the boundary-layer top.

Numerical aspects
For the run presented above, the advective, acoustic, and diffusive Courant numbers at the end of the simulation are respectively c ad = max(w 2 )Δt∕Δz = 0.47, c ac = c s Δt∕Δz ≈ 100, and c diff = max(K)Δt∕Δz 2 = 3.5, where c s ≈ 350 m/s is the sound speed. For the convergence test discussed below c ad and c ac are kept roughly independent of resolution while c diff roughly doubles with each halving of Δz, reaching around 30 on the finest grid. Thus, the numerical method presented in Section 4 can stably handle large acoustic and diffusive Courant numbers. However, despite the use of a semi-Lagrangian advection scheme, we found that the model became unstable when we attempted to increase the advective Courant number beyond 1. The symptoms of the instability suggest that it may be an acoustic mode, whose propagation is drastically slowed by the semi-implicit time stepping, coupled to advection by w 1 and w 2 . However, a more complete understanding and a solution have so far eluded us. A numerical method stable for large advective Courant numbers is desirable in order that a model time step is not constrained by the presence of strong updraughts. Despite solving a fully compressible equation set at very high acoustic Courant number, we found that the solution ) is the height of the minimum in * . t(ent flx) is the time at which a non-zero boundary-layer-top entrainment flux is first established. All diagnostics except t(ent flx) are taken at the end of an 8 hr run.
above the top of the boundary layer remains remarkably quiescent and noise-free, even as z * approaches the model lid. We suspect that this good behaviour is related to the form and discretization of the entrainment and detrainment terms, which conserve the mass integral of in each model layer and so preserve the pressure there.
Because of the sharp curvature of the potential temperature profile near the ground, the advective potential temperature tendencies at the first level above the ground are very large and very sensitive to the details of the advection scheme. The resulting errors have a significant effect on the updraught properties and hence on the solution throughout the boundary layer. For example, we experimented with linear, cubic Lagrange, and cubic spline interpolation, each using either z or z 2 as the coordinate, and found differences in updraught fraction, updraught w, and normalized mean potential temperature of order 10%.
On a closely related point, the growth of the boundary layer is linked to the budget of potential temperature, so accurate conservation of potential temperature is important. However, the semi-Lagrangian advection used for R i is not inherently conservative, so a numerical fixer is used to restore the correct mass integral after advection. Numerous conservation fixers have been proposed in the literature, including simple globally uniform additive or multiplicative factors, as well as more sophisticated schemes (e.g. Zerroukat and Allen, 2015). We found that our numerical solutions were very sensitive to the exact form of the fixer used, often giving spurious kinks in the profile of 2 or potential temperature flux. We obtained best results when the conservation fixer was applied entirely to the first level above the ground. This is consistent with the fact, noted above, that the errors in advection of are dominated by the first level above the ground. Including this conservation fixer also greatly reduced the sensitivity of the solution to the choice of advection scheme for .
To test the well-posedness of the two-fluid model and the ability of the numerical methods to solve it, we investigated the convergence of the solution with increasing resolution. Table 1 shows the time and space resolutions used along with a selection of key diagnostics. The finest resolution run was also used as a reference solution allowing errors to be estimated at coarser resolutions ( Figure 6; also Figures 4 and 5). Many diagnostics, such as z * , minimum 2 , maximum w 2 , and the root mean square errors, show a clear convergence with increasing resolution. However, the maximum in 2 w 2 and the height of the minimum in * appear to be converging rather slowly, if at all. These two quantities are related, since a larger mass flux produces a larger advective potential temperature transport and leads to a minimum * at a lower height. We have found the 2 profile, and hence 2 w 2 and z(min * ), to be sensitive to both the parametrizations and the details of the numerics near the bottom boundary. The slow convergence of these diagnostics, despite the modifications to numerics already introduced in Section 4, suggests the need for an improved and perhaps quite different numerical handling of the near-singularities at the bottom boundary.
For comparison, in the LES the normalized height of the mean potential temperature minimum z(min * )∕z * is quite variable in time, for example varying between 0.28 and 0.44 over a couple of hours. Because the middle part of the boundary layer is very well mixed, very small changes in the mean potential temperature are enough to significantly change the height of its minimum. It is not clear whether this variability indicates sampling errors due to the finite LES domain size, or a genuine physical variability. In either case, the LES results for this diagnostic do not provide a robust reference value for the single-column model results.
Some diagnostics do not show a convergence with increasing resolution, particularly the maximum errors in w 2 and 2 , and the time at which a clear entrainment flux is first established t(ent flx). The maximum errors in w 2 and 2 reflect the fact that the solutions for w 2 and 2 are discontinuous at z = z * , so that even a small error in z * implies a finite error in w 2 and 2 . The initial penetration of the updraught above its level of neutral buoyancy to establish an entrainment flux appears to be very sensitive to the details of the numerical handling of advection, diffusion, and detrainment near z = z * . Small changes to the profile of diffusivity, such as allowing it to vanish just above z * rather than exactly at z * , or small changes to the factors 1.001 in Equation 21 or 0.01 in the definition ofw, can lead to quite different behaviour. A further issue is that simply diagnosing z * becomes problematic if the w 2 profile is noisy. Our preferred approach for future work would be to avoid parametrizations that depend explicitly on a diagnosed boundary-layer height z * . We anticipate that an eddy diffusivity derived from a TKE (e.g. Angevine, 2005;Witek et al., 2011a;Sušelj et al., 2012;Tan et al., 2018) will lead to better convergence with resolution.

CONCLUSIONS AND DISCUSSION
We have presented a two-fluid single-column model for the dry, shear-free, convective boundary layer in which non-local transports by coherent structures such as thermals are represented by the partitioning of the fluid into two components, updraught and environment, each with a full set of prognostic dynamical equations. Local eddy diffusive transport and entrainment and detrainment are represented by parametrizations similar to those used in EDMF schemes. A semi-implicit semi-Lagrangian numerical solution method is presented and shown to provide stable solutions for large acoustic and diffusive Courant numbers, though it becomes unstable for large advective Courant numbers. In many aspects the solutions obtained are similar to those obtained with EDMF schemes. In particular, they are able to capture the countergradient potential temperature transport in the upper half of the boundary layer, the location of the mean potential temperature minimum in the middle of the boundary layer, and the occurrence of a boundary-layer top entrainment flux of order 0.1 times the surface potential temperature flux.
The work has highlighted several valuable lessons to be borne in mind for future development. One lesson concerns the need to include vertical diffusion of w in order to control an inherent Kelvin-Helmholtz-like instability of the two-fluid governing equations. It would be interesting to explore whether this instability could usefully be linked to the TKE budget and hence to the parametrization of eddy diffusivity and perhaps to entrainment and detrainment. Another lesson concerns the importance of conservation of mass and potential temperature for obtaining accurate numerical solutions.
We strongly advocate testing the convergence with resolution of any proposed parametrization (as well as testing at the proposed operational resolution) as a means of disentangling the properties of the underlying mathematical model from those of the particular numerical implementation (and any numerical "fixes" needed to make it work). Our results, especially those of the convergence test, have highlighted several areas where further understanding or development is needed.
First, there are several issues related to the use of parametrizations that depend explicitly on a diagnosed boundary-layer height z * . The method used here to diagnose z * can be vulnerable to noise in the w 2 profile. Moreover, the parametrizations described in Section 2 result in the diffusivity being identically zero above z * and the detrainment rate blowing up at z * , while the advective form of the w and equations means that their numerical advective tendencies are identically zero at any level where w = 0. As a consequence, the solutions for w 2 and 2 are discontinuous at z = z * and the inversion in * is unrealistically sharp, while the initial penetration of the updraught above its level of neutral buoyancy to establish an entrainment flux is very sensitive to resolution. A further point is that the parametrizations used here would be unsuitable for the moist case, in which some fraction of the updraught should be permitted to penetrate above the boundary-layer top to form cumulus cloud. For future work we propose to investigate parametrizations that do not depend explicitly on z * and which produce smoother, more realistic solutions.
Second, the solutions for both w i and i vary sharply near the ground. This impacts the numerical solution in several ways, including the calculation of semi-Lagrangian departure points and the interpolation of advected fields to those departure points, the calculation of gradients in diffusion terms, and the vertical interpolation of fields between w-levels and -levels. Unlike the rather artificial singular behaviour near the top of the boundary layer, which might be ameliorated by alternative parametrizations, this near-singular behaviour near the ground appears to be inherent to the mathematical description of the convecting boundary layer (e.g. Holtslag and Nieuwstadt, 1986;Troen and Mahrt, 1986;Holtslag and Moeng, 1991). Numerical methods can often be adapted to better handle singularities if the form of the singularity is known. The well-known convective boundary-layer scaling of w and near the ground motivated our introduction of the stretched coordinates z 1 , z 2 and z 3 in Section 4. This introduction led to some improvements, and we were able to obtain w i profiles that scale like z 1∕3 . However, we were not able to obtain the correct scaling for i . A complicating factor is that the entrainment source term added to i before advection scales in a different way from i itself. A related point is that, while semi-Lagrangian schemes give excellent accuracy when the Lagrangian time-scale is long, as is typical of the free atmosphere, they can be much less accurate when the Eulerian time-scale is long and a quasi-steady balance between strong source terms and advection is required, as here. Since these numerical errors near the ground affect the updraught properties and hence the solution throughout the boundary layer, it would be desirable to develop improved numerical schemes that can accurately capture the solution near the ground.
Finally, as mentioned above, despite the use of semi-implicit semi-Lagrangian numerics, the solution becomes unstable for large advective Courant numbers. Because updraught vertical velocities are typically much larger than grid-cell mean vertical velocities, using the multi-fluid approach in a weather or climate model with the numerical method described here would require a significant reduction in time-step size compared with current approaches. Therefore it would be desirable to understand the mechanism of this instability and to develop an improved numerical method that is more stable.
Despite these areas in need of further work, the results indicate that the multi-fluid approach can capture some of the essential physics of a convecting fluid. Therefore we are encouraged to extend the approach to two and three dimensions and to include moist processes in order to investigate whether the potential benefits of the approach mentioned in Section 1 can indeed be realized.

This work was funded by the Natural Environment Research
Council under grant NE/N013123/1 as part of the ParaCon programme. We are grateful to two anonymous reviewers, whose constructive comments helped improve this article.

A. TWO-FLUID KELVIN-HELMHOLTZ-LIKE INSTABILITY
It has proved challenging to obtain numerical solutions of the multi-fluid equations. Such problems are familiar from attempts to solve similar equation sets in engineering applications, and have been linked to the fact that (in the absence of their right-hand sides) the equations are ill-posed in the sense of possessing complex characteristic speeds (e.g. Stewart and Wendroff, 1984).
Another manifestation of this "ill-posedness" is that a simple basic state with uniform but different velocities in the two fluids is unstable. Here we demonstrate this instability in the simplest possible scenario: a one-dimensional, two-fluid, incompressible flow. The governing equations are then ∑ Introduce a steady basic state with constant volume fractions 0 i and velocities W i , where 0 1 + 0 2 = 1. Small perturbations to this basic state satisfy the linearized equations ∑ Seeking solutions proportional to e i(mz− t) and eliminating unknowns leads to the dispersion relation The frequency is real if one of 0 i vanishes or if W 1 = W 2 ; in all other cases the basic state is unstable. The growth rate is proportional to the wavenumber m and to the velocity difference |W 2 − W 1 |.
Recall now that Equations A1-A3 describe the dynamics of a conditionally filtered system. There are infinitely many possible pre-filtered states that would give the same conditionally filtered basic state. One obvious choice is a two-dimensional (x, z) basic state (where the filter is an average in x) in which u = 0 and w is a function only of x, taking only the values W 1 and W 2 in alternating layers of width D 1 and D 2 , respectively, where D 1 ∕D 2 = 0 1 ∕ 0 2 . For such a state, all of the right-hand side terms in the conditionally filtered governing equations would indeed be zero, at least initially. It is well known that such a state would be unstable to Kelvin-Helmholtz instability. It is reasonable to conjecture, therefore, that the instability found in the two-fluid model is a conditionally filtered representation of this Kelvin-Helmholtz instability.
To support this conjecture we carried out a stability analysis of this two-dimensional state. The analysis follows standard methods (e.g. Drazin and Reid, 1981) so the details are omitted. We consider the limit mD i → 0; this is the relevant limit for comparison with the two-fluid model, since the wavelength 2 ∕m is resolved but the width D i is subfilter-scale. Then, identifying D i ∕(D 1 + D 2 ) with 0 i , the dispersion relation is found to be (A8) Both pairs of solutions have the same growth rates as predicted by the two-fluid model, in particular the same dependence on the shear |W 1 − W 2 |, the volume fractions 0 i , and the wavenumber m. The second pair of solutions also has the same propagation speed as predicted by the two-fluid model. The two-fluid model does not have enough degrees of freedom to capture the first pair of solutions.
Interestingly, the solutions captured by the two-fluid model have the property that, when one of the volume fractions is small, the propagation speed approaches the basic state velocity of the smaller fluid fraction. Thus, the two-fluid model is remarkably successful at capturing some key aspects of this Kelvin-Helmholtz instability.
For the two-dimensional instability considered here, the flow would quickly become turbulent; then the subfilter-scale fluxes and pressure terms on the right-hand side of Equations 1-4, which were omitted in Equations A1-A3, would no longer be negligible.
As a thought experiment, one can also consider a dry convective boundary layer that is statistically homogeneous in the horizontal. A conditionally filtered description of such a flow (where the filter could be a horizontal mean) would evolve slowly, and not show the kind of dramatic transience predicted by Equations A1-A3 and A7. Clearly the subfilter-scale terms, and perhaps the relabelling terms, must be playing a crucial role in stabilizing the conditionally filtered equations. We therefore argue that the subfilter-scale terms are an indispensible part of the multi-fluid equations, and that useful solutions cannot be expected if they are omitted.
A mutual drag between fluid components could be expected to provide some damping of the two-fluid instability. However, the form given by Equations 31-32, which is independent of the vertical scale of any perturbations, is unable to control the instability, whose growth rate is proportional to vertical wavenumber. Similarly, entrainment of or w as in Section 2.3 provides only scale-independent damping and cannot control the instability. Subfilter-scale fluxes of potential temperature, parametrized, for example, as a vertical diffusion, are also unable to control the instability. (Consider the case of uniform potential temperature, when such fluxes would have no effect.) However, subfilter-scale fluxes of vertical velocity, parametrized as a vertical diffusion, would preferentially damp the smallest vertical scales. Moreover, a comparison of the damping rate based on Equation 7 with the growth rate given by Equation A7 indicates that such a diffusion of vertical velocity would be sufficient to control the instability on scales smaller than the boundary-layer depth z * . This argument motivates our inclusion of a vertical diffusion on w i in Section 2.2.

B. APPROXIMATE LINEARIZATION OF EQUATIONS 33-35
A full linearization of Equations 33-35, with the left-and right-hand sides given by Equations 36-41, would include a number of less stiff terms that are not crucial to the convergence of the Newton iterations, but which would make analytical elimination of unknowns infeasible. Therefore we neglect several of the less stiff terms. Note that R i , R i , and R w i depend only on the known time step n fields, so the only increments involved in the linearization of the right-hand sides of Equations 33-35 are the w ′ i arising through the semi-Lagrangian departure point calculation.
For the left-hand sides of Equations 45-47, we use