{'fraction_non_alphanumeric': 0.035449975187175006, 'fraction_numerical': 0.01637646449608389, 'mean_word_length': 4.373058196151171, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 16, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Numerical analysis of the multi-fluid equations with applications for convection modelling

Convection schemes are a large source of error in global weather and climate models, and modern resolutions are often too fine to parameterize convection but are still too coarse to fully resolve it. Recently, numerical solutions of multi-fluid equations have been proposed for a more flexible and consistent treatment of sub-grid scale convection, including net mass transport by convection and non-equilibrium dynamics. The technique involves splitting the atmosphere into various fluid components which are defined by their physical properties, such as the stable environment or buoyant updrafts. The fluids interact through a common pressure, but such equations are ill-posed, meaning extra terms such as drag or mass transfers are needed to stabilise them. Relatively little is known about the numerical stability of the multi-fluid equations or the numerical properties of mass transfer terms between the fluids. This study aims to provide an insight for these two topics: 1. We analyse the physical and numerical stability of the multi-fluid shallow water equations with a common height gradient. We show that the 2-fluid linearized equations are less stable than the 1-fluid case. However, the full, non-linear equations have a greater numerical stability due to a direct energy cascade which removes energy from the most unstable modes. 2. We derive mass transfer terms which relabel the fluids and derive numerical properties of the transfer schemes including boundedness, momentum conservation and energy properties. Numerical solutions of the multi-fluid Euler equations are presented using stable and unstable treatments of the transfers. We find three transfer schemes which possess desirable numerical properties, two of which are suitable for use on staggered grids.

into various fluid components which are defined by their physical properties, such as the stable environment or buoyant updrafts. The fluids interact through a common pressure, but such equations are ill-posed, meaning extra terms such as drag or mass transfers are needed to stabilise them. Relatively little is known about the numerical stability of the multi-fluid equations or the numerical properties of mass transfer terms between the fluids. This study aims to provide an insight for these two topics: 1. We analyse the physical and numerical stability of the multi-fluid shallow water equations with a common height gradient. We show that the 2-fluid linearized equations are less stable than the 1-fluid case. However, the full, non-linear equations have a greater numerical stability due to a direct energy cascade which removes energy from the most unstable modes. 1

| INTRODUCTION
The modelling of atmospheric convection is at the forefront of current meteorological research due to the large errors caused by current convection schemes in atmospheric models (eg Arakawa, 2004;Lean et al., 2008;Yano et al., 2004Yano et al., , 2018. We have reached the "grey zone" in which improved computational power allows convection to be partially resolved but cannot yet be explicitly simulated (Gerard et al., 2009). Many assumptions designed for convection schemes on large-scale grids are increasingly unrealistic at finer resolutions due to the increased significance of convective processes at smaller scales (Kwon and Hong, 2017). Processes such as net mass flux by convection are commonly ignored in convection schemes (such as Arakawa and Schubert, 1974;Gregory and Rowntree, 1990;Lappen and Randall, 2001) as well as non-equilibrium dynamics (for example Kain and Fritsch, 1990). Although many convection schemes incorporate some of these aspects (such as Gerard and Geleyn, 2005;Kuell and Bott, 2008), the conditional filtering (or conditional averaging) technique has been proposed for convection modelling due to the possibility of consistent treatment of resolved and sub-grid convection, whilst also representing net mass transport by convection and non-equilibrium dynamics Weller and McIntyre, 2019).
Conditional filtering involves dividing space into various defined fluids, each with their own physical properties . In the case of convection, one could use fluid 0 as the neutrally buoyant air, fluid 1 as convective updraft regions and fluid 2 as downdraft regions. Each fluid has its own properties and prognostic variables such as volume fraction, density, temperature and velocity. The equations of motion are solved for each fluid individually. As the scheme allows for the advection of any fluid to neighboring cells, net mass transport by convection can occur in which the properties of the convective mass are transported based on the local dynamics.
Conditional filtering has already been implemented in various fields of science and engineering (Dopazo, 1977;Méchitoua et al., 2003;Guelfi et al., 2007). Moreover, Lappen and Randall (2001) used the technique to model cumulus convection, using an updraft fluid and a stable fluid. However, the stable fluid in that study was assumed to subside in order to counteract any mass transfer due to the convective updraft, meaning the scheme does not incorporate net vertical mass transport by convection. More recently,  and Tan et al. (2018) have described how to conditionally filter the fully compressible Euler equations with the aim of representing sub-grid scale convection. Subsequent studies have built upon these foundations including  who investigate the conservation properties and normal modes of the equations, Thuburn et al. (pre-release: 2019) who use the method for a 2-fluid single-column convective boundary layer scheme, and Weller and McIntyre (2019) who formulate a numerical solution of the multi-fluid compressible Euler equations.
Little is known about the numerical properties of solutions to the multi-fluid equations. Stewart and Wendroff (1984) and Thuburn et al. (pre-release: 2019) note that the multi-fluid Euler equations are ill-posed when sub-filter terms are ignored. This property is confirmed by Weller and McIntyre (2019), as drag or diffusion between the fluids is necessary to prevent the fluids from diverging. Weller and McIntyre (2019) also explain that positivity of mass terms is crucial for numerical stability. However, none of these studies comprehensively analyse the differences in numerical behaviour of the 1-fluid and multi-fluid equations. It is also unclear what effect the non-linear advection terms have of the system behaviour, as Thuburn et al. (pre-release: 2019) analyse a linearized equation set. In addition, the drag and diffusion used in Weller and McIntyre (2019) are manually chosen to stabilise the system and are not derived through analysis. Numerical analysis of the multi-fluid equations is therefore important for informing scheme choices for multi-fluid convection modelling. This motivates our investigation of the multi-fluid shallow water equations which provide an idealised framework without a heat equation or equation of state. When advection terms are small and the shallow water equation set is linearized, the system reduces to acoustic wave equations, also present in the multi-fluid fully compressible Euler equations. Investigating the core instabilities associated with the multi-fluid shallow water equations will therefore inform us of some properties present in the full equation set for convection modelling.
A parameterisation of convection using conditional filtering also relies on formulating transfer terms that exchange mass and other properties between fluids. These transfer terms will be equivalent to entrainment, detrainment and the triggering of convection which may be adapted from existing frameworks including Betts andMiller (1993), Neggers et al. (2002) and Siebesma et al. (2007), for example. General fluid transfer terms (without physics) are given in  and Weller and McIntyre (2019). Weller and McIntyre (2019) also create a numerical scheme for the fluid transfers, but only one transfer scheme is considered in which the numerical treatment of the mass transfer is explicit (and the momentum/temperature transfer is treated implicitly) which may not be suitable for all transfer rates. This motivates us to present more mass transfer schemes and analyse their numerical properties to obtain the most desirable numerical solutions.
In this study, we derive the multi-fluid shallow water equations (section 2) and conduct stability analysis on the linearized equation and subsequently investigate the properties of the non-linear equation set (section 3). We also analyse the numerical properties of the transfer terms between fluids, applicable to both the multi-fluid shallow water and multifluid Euler equations. We will discuss the schemes with the most desirable numerical properties including conservation, boundedness and stability, and present numerical solutions of the multi-fluid Euler equations using these transfer terms (section 4).

| GOVERNING EQUATIONS
The multi-fluid fully compressible Euler equations are derived in . In this section, we derive the inviscid, non-rotating multi-fluid shallow water equations using the methodology from . The single-fluid equations are given by the continuity equation, and the momentum equation which may be represented in flux form by u is the velocity of the fluid, h is the fluid height and g is the acceleration due to gravity. In order to convert equations (1) and (2) into multi-fluid form, we must first consider how to distinguish between fluids at any scale using the approach outlined by . At any infinitesimal point in space, we assume only one fluid (i ) will occupy that position. We label this point using the indicator I i = 1. We also label the absence of any other fluids at this point using The conditionally-filtered distribution of some function F (x) is given by the convolution: where G (x) provides the domain-weighting to F (x ) (such as volume averaging) with the condition that ∫ D G (x) d x = 1. By choosing F (x) to be proportional to I i , we make F (x) the volume-weighted distribution of some property of fluid i .
For example, the volume fraction of a fluid is obtained by σ i = I i . By applying the conditional filter to the governing equations, assuming the filter commutes with temporal and spacial operators, and assuming each fluid shares the same height (h i = h), we get the continuity equation, and the momentum equation in flux form, where u i is the velocity of fluid i , F u i S F is the sub-filter scale flux and D i j is the drag between fluids i and j . In this study, we are only considering the numerics of resolved states (as in Weller and McIntyre, 2019) and thus omit the sub-filter flux terms. We will also ignore the drag terms in order to investigate the core numerical properties of the equation set.
Our idealised momentum equation in flux form is therefore which, in the equivalent advective form, is given by Note that our assumed shared height between fluids results in a single height-gradient acting on all of the fluids. This is similar to the shared pressure assumption made by Stewart and Wendroff (1984),

| STABILITY ANALYSIS
Various studies on the physical and numerical stability of the shallow water equations have been conducted (for example Stelling, 1983;Rebollo et al., 2003;Anmala and Mohtar, 2011) but this is the first study of the multi-fluid shallow water equations. Studies on the stability of other multi-fluid equation sets include Stewart and Wendroff (1984), who show that the two-phase Euler equations with a common pressure are ill-posed and thus inherently unstable. Thuburn et al. (pre-release: 2019) showed that there exists a Kelvin-Helmholtz-like instability for the linearised 1D two-fluid incompressible equations of motion unless 1 fluid vanishes or their velocities are equal. In this section, we will investigate whether the same instabilities exist in the linearised 1D 2-fluid shallow water equations and assess the implications for the non-linear equation set.

Linearized equations
In order to linearize the multi-fluid shallow water equations (necessary for linear stability analysis), we assume that each of our prognostic variables (η i ≡ σ i h and u i ) is equal to the sum of a mean state (denoted by the overbar) and a perturbation (denoted by the prime): We assume that the mean state variables do not vary in time or space. Substituting into our governing equations and neglecting any second-order perturbation terms, we get the linearized multi-fluid shallow water equations,

Physical stability of linearized equations
In order to diagnose the physical stability of our 1D 2-fluid equation set, we assume wave-like solutions of the form where H i and U i are constants, M is the wavenumber, ω is the angular frequency and u i is the x-component of u i . For oscillatory behaviour (and therefore physical stability of our system), we expect M , ω and K ≡ ω M to be real. Substituting the solutions into 10 and 11 for a 2-fluid system with i = {0, 1}, eliminating unknowns and re-arranging for K , we get the quartic equation, which may be solved numerically or analytically to find stable values of the mean heights and velocities. Note that for the special cases where η 0 = 0 (or η 1 = 0 -only 1 fluid present) or g = 0 ms −2 (fluids no longer coupled -2 1-fluid equations), K can be solved using the quadratic equation, giving the 1-fluid dispersion relation K = u i ± g (η 0 + η 1 ) which is real and therefore stable for all g , η i ≥ 0. We also get the same relation when u 0 = u 1 where the mean fluid properties are the same but the perturbations may differ; a quasi-one-fluid state.
We assess the physical stability by numerically solving the two-fluid dispersion relation (12). If any of the four solutions has a positive imaginary component, the system is physically unstable. Figure 1, shows various examples of the stability regions for c 1 = u 1 ∆t ∆x ∈ [−1.1, 1.1] (220 increments), c g = g (η 0 + η 1 ) ∆t ∆x ∈ [0, 1.1] (110 increments). Parameter-space with physical stability is given by the red-toned and white regions. Regions of physical instability are given by the blue and grey-toned regions. We have also included points which are unstable but within the numerical round-off error of the stable/unstable boundary (denoted by ∼ unst abl e, the lighter/darker shades of grey, red and blue). This error estimate is given by the Bauer-Fike theorem (Bauer and Fike, 1960) which measures the error propagation through eigenvalue calculations. Using 64-bit floats, the error ranges are typically seen to be E ≈ 10 −8 . For physical stability, we therefore use this error cut-off when the imaginary part of K is less than 16E . The choice of 16E stems from the assumption that each element of four eigenvectors may have its own rounding error.
For the single fluid cases (figure 1.a), we verify that the linear shallow water equations are stable for g (η 0 + η 1 ) ≥ 0, as can be shown analytically. However, the single fluid stability regions do not extend to the two-fluid case where the velocities of the fluids differ, as shown in figures 1.b-d. Instead, we observe physical stability when the two fluid velocities differ substantially, as well as only an infinitesimal stable region where the velocities are equal. At its smallest, the (red/white) physical stability region is bounded by c g ≤ 1 2 |c 0 − c 1 | (to within machine accuracy) which occurs when η 0 = η 1 = 0.5 m (figure 1.b). The stable parameter space extends to c g ≤ |c 0 − c 1 | as η 0 η 1 → 0 (see figures 1.c-d) but returns to the single-fluid distribution once η 0 η 1 = 0. The volume of the physically stable 2-fluid parameter-space is, at best, less than half that of the 1-fluid system, indicating that the 2-fluid shallow water model also introduces the instabilities described by Stewart and Wendroff (1984).

Numerical stability of linearized equations
We choose to analyse the numerical stability of the linearized equation set using first-order upwind advection, forwardbackward timestepping and a C-grid. Weller and McIntyre (2019) showed that total variation diminishing (TVD) schemes are necessary to keep our prognostic mass variables positive (which maintains numerical stability) and firstorder upwind is the only TVD scheme amenable to linear stability analysis. Moreover, using a staggered rather than co-located grid helps to prevent grid-scale oscillations caused by computational modes. Numerically discretizing (10 and 11), we getη where we have definedη i ≡ as our dimensionless prognostic variables. ∆t is the timestep, ∆x is the spacial discretization and variables attached to the vertical bars indicate the numerical labels: the superscript gives the time level and the subscript the spacial level so that x = k ∆x and t = n∆t . Note that we have used a form of the first order upwind interpolation valid for both positive and negative velocities (used by Rebollo et al., 2003). Due to |u i |, the advection term reduces to the familiar first order advection scheme for both u i > 0 and u i < 0.
We define the vector v T | n k ≡ η 0 | n k ,η 1 | n k ,ũ 0 | n k + 1 2 ,ũ 1 | n k + 1 2 T and conduct Von-Neumann stability analysis by searching is an integer wavenumber and A is the amplification matrix. In order to find A such that v | n+1 where B, C, D, E, F and G are also matrices, obtained using equations (13) and (14). We then substitute v | n+1 where I is the identity matrix. For numerical stability, we require that ρ(A) (the spectral radius of the amplification matrix) is less than or equal to 1 for all wavenumbers. To assess the numerical stability, we calculate the spectral radius of the amplification matrix (A) for various wavenumbers, M ∆x ∈ [0, 2π] discretised by π/20 -this is equivalent to a spacial domain with 40 cells, which can support 20 wavemodes. If any eigenvalue for any value of M ∆x is greater than 1, then the numerical scheme is unstable for the given mean values. Note that numerical stability without physical stability should not occur for consistent numerical methods according to the Lax equivalence theorem (Lax and Richtmyer, 1956). However, as we are using a non-infinite wavenumber resolution, it is possible that some false-positive numerically stable results will occur. We compute the spectral radius for all the same parameter values of {η i , u i } as the physical stability analysis, also shown in figure 1.
Regions of numerical stability are given by the white and blue-toned colours. We have once again included points which are unstable but within the numerical round-off error, when the spectral radius is less than 1 + 16E , indicated by the ∼ unst abl e tag. In the 2-fluid cases (figures 1.b-d), we cannot expect better stability than the one-fluid stability relation: This maximum stability region is shown by the dotted lines in figures 1.b-d. Our analysis shows that the numerical stability largely coincides with the overlap between the 1-fluid stability triangle and the physical stability zones. This pattern is most clear in figure 1.c.ii as we are "close" to a one-fluid system with η 1 η 0 . However, the distributions deviate from this trend when either velocity approaches 0 ms −1 , as can be seen in all 2-fluid figures for c 0 → 0 (figures . In all of these cases, the numerical stability-range is weakened with respect to the physical system (red). The amplification matrix becomes undefined as u i → 0 which indicates the scheme is numerically unstable when either fluid velocity is zero. Where we have η 0 η 1 (figures 1.d.i and 1.d.ii), we also see different behaviour for the numerical stability which is no longer confined within the 1-fluid stability triangle.
This is due to the small quantity of η 1 present, meaning the total system stability is more dependent on fluid 0. Instead, the fluid 1 stability regions can be though of as a "slice" of the fluid 0 stability triangle.
The final discrepancy occurs for c 0 = 0.505 where the numerical scheme is stable when the system is physically unstable, as indicated by the blue regions. However, these blue regions become narrower and eventually disappear when using higher resolutions of M ∆x . It is always the smallest wave modes (largest wavelengths) that are unstable (see figure 2) and so sufficient resolution is needed to resolve them. The numerical stability of the 2-fluid linearized shallow water equations is therefore consistent with the physical stability of the equations. There is clear reduction in numerical stability relative to the 1-fluid system, but this may not be fully representative of the numerical properties of the non-linear equation set (equations 4,7), as will be discussed is section 3.4.

Numerical simulations of the linear and non-linear equations
The implications of the linear equations' stability on the non-linear multi-fluid shallow water equations are not clear: 1. Should the non-linear equations be assumed unstable because the linear equations are not stable in the given parameter-space?
2. What do we infer from the observation that only the long wavelengths are unstable?
3. How do the non-linear advection terms influence the unstable large-scale modes?
To investigate this further, we can observe the energy evolution of the discretized linear and non-linear systems as a numerical method will be unstable if the total energy increases with time. The total energy of the non-linear shallow water equations (per unit density-area) is given by the sum of the total potential energy (P ) and total kinetic energy (K ), where and respectively. Following the single-fluid methods in Ringler et al. (2010), one can show that the multi-fluid shallow water equations conserve the total energy via where the divergence term cancels when integrated over the entire domain using the divergence theorem and assuming periodic boundary conditions. The 2-fluid linearized shallow water equations do not conserve the same energy quantity.
Instead, we can derive that is a conserved quantity via where k i ≡ u i .u i and k i ≡ 1 2 u i .u i . The curl terms also disappear when integrating over the domain (assuming periodic boundary conditions) using an alternative form of Gauss' divergence theorem: where V is an arbitrary volume and S is the surface of that volume. This means that we can use E lin as our conserved "energy" quantity for the linearized equations. Note that any first-order perturbation terms can also be added to E lin without affecting conservation.
The numerical scheme we choose for the non-linear equations also uses forward-backward timestepping, first order upwind advection and a staggered grid, as given by: where and For our numerical simulations, we choose a 1D domain in the x direction with x ∈ [0., 1.] m, ∆x = 1 40 m and periodic boundary conditions as well as a timestep of ∆t = 0.01 s. For any mean parameter values Q (Q ∈ [η i , u i ]), we define initial square-wave perturbations: Square wave perturbations are chosen as they encompass all wavemodes. We run both the linear and non-linear schemes for 5 seconds and compute the total energy quantities, given by equations (18,22), for each timestep. Note that, even though the physical equations conserve energy, the numerical schemes are not formulated to conserve E lin and E non . Figure 3 shows snapshots of these simulations for {η 0 = 0.5m, η 1 = 0.5m, c 0 = 0.005, c 1 = 0.1, c g = 0.1}. A scheme is labelled unstable if the total energy at the end of the simulation is greater than the initial energy. Otherwise the scheme is classified as "stable", although this may cause false positive results if instabilities are present but have not yet become prevalent in the 5 second simulation. Figure 1 shows the stability contours for the linear (cyan) and non-linear (orange) numerical simulations. The areas within these contours are classified as stable. We immediately see strong concordance in stability with the analysis and simulations of the linearized system. This is particularly evident in all figures (except for 1.c.i) where the cyan contour encompasses the white and blue predicted stability regions, including the stability deficits as either fluid velocity tends towards 0 ms −1 . There are some regions where the numerical simulation is stable but the analysis suggests that the scheme should be unstable. In these cases, the linear simulations have not become unstable yet according to our detection criteria. We observe energy increases for a few selected cases of these simulations when run for 50 s rather than 5 s (not shown). Unlike the linear scheme, the non-linear simulation stability region does not reduce with 50 s simulation times, suggesting that the orange stability contour is representative of the non-linear stability.
The non-linear simulations are energetically stable for a greater range of parameter-space relative to the linearized equations (see orange lines in figure 1). Perhaps through coincidence, the stable region is approximately concordant where Z 0 = |c 0 | + σ 0 |c g |, Z 1 = |c 1 | + σ 1 |c g | and Z 2 = |σ 0 c 0 + σ 1 c 1 | + |c g |. We do not mathematically justify the use of Z non-linear in this study but instead outline that it closely resembles the non-linear stability region. However, we could reason that Z 0 and Z 1 represent the stability conditions for fluids 0 and 1 respectively, where each fluid is only influenced by a fraction of the height gradient term proportional to how much of the fluid is present. Z 2 is a statement of the 1-fluid stability relation, where the grid-box mean velocity is used -this once again implies that the gravity waves act on the overall fluid mass. This reasoning is only speculation and it can not explain the regions where simulation stability is not achieved. Further understanding of the stable non-linear regions should be the subject of future work.
The enhanced stability of the non-linear scheme relative to the linear scheme suggests that the non-linear advection terms are removing energy from the system. This allows for greater numerical stability, as displayed by the energy time series in figure 3 for simulation L1. Energy depleting effects are noted in Leonard (1975), for example, where the non-linear advection term for resolved scales is responsible for an energy cascade from large scales (low wavenumbers) to small scales (large wavenumbers). At small scales, energy is subsequently dispersed if viscosity is present or lost numerically if the scales are smaller than the grid scale. We know from figure 2 that the low wavenumbers are unstable (in the solid line case where max [|c 0 |, |c 1 |] + |c g | ≤ 1). Therefore, a direct energy cascade would explain the improvement in numerical stability. Note that turbulent flow in two dimensions (as is approximately the case in the atmosphere at scales greater than ∼ 10 km) is known to have a net inverse energy cascade to the larger scales (Boffetta, 2007;Yano et al., 2018). However, 3-dimensional turbulence (as present within clouds) produces a net energy cascade to the smaller scales, like our 1-dimensional case (Tennekes, 1978;Yano et al., 2018). We must therefore be cautious in extrapolating these results to the 2D and 3D cases. However, we can generally expect the non-linear multi-fluid shallow water equation scheme to be more stable than the linear scheme due to the the non-linear advection terms which cause a direct energy cascade to the small scales where modes smaller than the grid scales are destroyed.

Form of the mass transfer terms between fluids
If S i j is the transfer rate from fluid i to fluid j , then the continuity equation for each fluid becomes For the multi-fluid shallow water equations, we use η i = σ i h. For the fully compressible multi-fluid Euler equations we instead use η i = σ i ρ i as the prognostic mass quantity, where ρ i is the density of fluid i (Weller and McIntyre, 2019). As we are considering resolved systems, we assume the mass transferred between fluids will take its associated mean properties from the original fluid. Refinement of the transfer terms to incorporate sub-filter-scale variation can later be incorporated and does not form part of this study. For a mean prognostic variable φ i , the governing equation is given by in flux form and in advective form. F i contains additional right-hand-side terms such as the height gradient term in the multi-fluid SWE momentum equation. We can use this template to construct the following governing equations: • Momentum equation for the multi-fluid shallow water equations: η i = σ i h, φ i = u i , F ui = −g η i h. This gives equations 6 and 7 with velocity/momentum transfers.
• Momentum equation for the multi-fluid Euler equations: . Ω is the angular velocity vector of the domain, c p is the heat capacity, θ i is the temperature of fluid i , π is the Exner pressure, g is the gravitational acceleration and d i j is the drag between fluids i and j .
• Temperature equation for the multi-fluid Euler equations: transfer between fluids i and j .
We will assume that the transfer terms are operator split such that other processes (such as advection) act on the prognostic variables first and the mass transfers are then based on their most up-to-date states (which we will denote by the time-level m). For each momentum equation, the total momentum should be equal before (m) and after (n + 1) mass transfer such that: i η i | n+1 u i | n+1 = i η i | m u i | m . The internal energy should also be conserved by the temperature equation transfers: c v π i η i | n+1 θ i | n+1 = c v π i η i | m θ i | m , where c v is the heat capacity at constant volume.

Mass transfers for the continuity equation
When transferring mass, we must ensure that all η i remain positive and that the total energy does not increase (we will consider other conservation properties in sections 4.3 and 4.4). The mass in each fluid at the end of the timestep, η i | n+1 , is given by the mass after advection (η i | m ) plus the discretised transfer term integrated over time ∆t : where η i | m are the values of η i after advection and η i | n+1 are the values after mass transfer. α C determines whether the transfer terms in the continuity equation are numerically treated explicitly (α C = 0) or implicitly (α C = 1). Re-arranging these equations for η i | n+1 , we get where and A is a label used to identify between transfers in the continuity (C), momentum (M) and temperature (T) equations.
The total mass is clearly conserved as i η i | n+1 = i η i | n . Therefore, the total potential energy (19) is conserved. As λ i j is between the values of 0 and 1 for ∆t S ≤ 1, η 0 and η 1 remain positive, ensuring numerical stability of the scheme.

φ transfers -Method 1
We must also model the transfer of velocity and temperature associated with the re-labelling of mass between fluids described in section 4.2. The new value of some variable φ i (velocity or temperature of the fluid, for example) should be bounded by the old values so that new extrema are not generated . Also, momentum should be conserved for best accuracy and energy should not increase for numerical stability. We assumed in equation 33 that any mass transferred takes with it the grid-box-mean variable from the fluid it originated from. Numerically, this equates to where α A = 0 means φ i is treated explicitly and α A = 1 means φ i is treated implicitly. Note that these equations have additional degrees of freedom in the time-level choice for η i , where q and r are the time level choices for the numerator and denominator respectively. A is the label for each governing equation. For the momentum and temperature equations, we will use A = M and A = T respectively. Rearranging for φ i | n+1 , we obtain where ν Ai j | q,r = ∆t S i j With this transfer scheme, momentum and internal energy are conserved when: 1. α C = 0, α M = α T = 0 with q = m, r = n + 1.
If ν i j | q,r ∈ [0, 1], then the new values of φ will remain bounded. This can only be guaranteed when α A = 1 (we prove this in appendix section 6.1). Moreover, appendix 6.2 describes tests which imply that schemes 1, 2 and 4 do not increase the total kinetic energy of the system for ∆t S i j ≤ 1. Although schemes 2 and 4 fulfil all of the desired numerical properties of our transfer terms, we have not explicitly proven that energy decreases, meaning alternative schemes should also be considered if they have energy-diminishing properties.

φ transfers -Method 2
An alternative approach takes advantage of the flux form equations (32) which uses the mass-weighted variable as the prognostic quantity (η i φ i ). These transfers unconditionally guarantee the conservation of η i φ i (such as momentum and internal energy) as whatever leaves one fluid enters another. Starting with our prognostic variable φ i , we define our mass-weighted quantity as Φ i ≡ η i φ i ) and apply the transfers according to equation 32: The equation takes a similar form to (34), meaning we get the solution In order to diagnose the value of φ i , we have chosen to divide Φ i by the new value of the effective mass. We have proposed this alternative method as we can demonstrate that φ i is bounded for any α C = α A (see appendix 6.1). We also show that this transfer scheme never increases the total kinetic energy of the system when mass and momentum transfers are treated consistently (see appendix 6.2). Note that, when using different combinations of α C and α A , the 4.5 SUMMARY OF PROPOSED TRANSFER TERMS total energy increases and unbounded values can arise. Using purely explicit or purely implicit treatments, we therefore have 2 more viable transfer schemes for the multi-fluid equations: 5. α C = 0, α A = α M = α T = 0.
Note that schemes such as [α C = 0.5, α A = 0.5] can also be used but there is no increase in order of accuracy as the scheme is operator-split and thus the time level m is not that of the previous timestep. Using time level n instead of m introduces instabilities into the numerical method as updates from the prognostic equations such as advection will be ignored in the transfer scheme.

Summary of proposed transfer terms
There exist 6 numerical transfer schemes which maintain positivity of mass and conserve mass, momentum, potential energy and internal energy. These schemes are presented in table 1. Some of these schemes also keep the prognostic variables bounded and never increase the total energy of the system (as we have shown that the kinetic energy can only decrease).
Only schemes 2, 4 and 6 are kinetic-energy-diminishing for all parameters meaning schemes 1, 3, 4 and 5 can cause numerical instabilities if ∆t S i j is large. We should also note that additional degrees of freedom arise when using a staggered, rather than co-located, grid. For example, if we use a staggered grid, with η i defined at cell centres and {φ i , Φ i } at cell faces, we could either conduct all our transfer terms at cell centres (by interpolating φ i or Φ i ), at the faces (by interpolating η i ) or instead by interpolating the transfer term S i j . This first two examples can cause large energy increases using any transfer scheme due to inconsistencies in the transfer rates which arise from interpolating S i j (not shown). The latter case is investigated in section 4.6.

Rising bubble test case
In order to test the properties of the various transfer schemes on a staggered grid, we have implemented them in a 2-fluid test case solving the multi-fluid fully compressible Euler equations, described in Weller and McIntyre (2019).
The test case is an adaption of the single-fluid rising bubble test case outlined in Bryan and Fritsch (2002). We use a 2D C-grid with η i and θ i defined at cell centres and u i defined at cell faces. The domain extends to x ∈ [−10, 10] km and z ∈ [0, 10] km with uniform grid spacings ∆x = ∆z = 100 m and wall boundaries on all sides. A uniform potential temperature field of θ 0 = θ 1 = 300 K is chosen with the system in hydrostatic balance and zero velocity. A warm temperature perturbation is applied to fluid 0 at t = 0 s: The perturbation is only applied for L ≤ 1 where L ≡ x −xc xr + z −zc zr , x c = 10 km, z c = 2 km and x r = z r = 2 km. A van-Leer advection scheme is chosen to maintain positivity of the mass. All details of the numerical setup and the numerical solvers used are described in Weller and McIntyre (2019). The expected evolution of the 1-fluid test case is presented in figure 4.a, where the warm bubble rises and generates symmetric convective rotors.
For the 2-fluid test case, both fluids initially have identical properties with the exception that no fluid 1 is present. The magnitude of the transfer terms is also borrowed from Weller and McIntyre (2019): where S d i j diffuses the two fluids together via where K σ = 100 m 2 s −1 is an arbitrarily-chosen diffusion coefficient in order to keep the simulation numerically stable.
S bi j transfers warm perturbations from the "neutral" fluid 0 to the "buoyant" fluid 1 and vice versa for cold perturbations in fluid 0: where K θ = 10 6 m 2 s −1 is a diffusion coefficient. For the velocity transfers, S i j and η i are linearly interpolated onto the cell faces where the transfer coefficients ν M i j and λ M i j are calculated.
The test case is run for all transfer schemes, including the 14 non-momentum-conserving schemes. The temperature anomalies of schemes 2 and 6 relative to the 1-fluid test case are presented in figure 4. The distribution of the temperature anomalies indicates that kinetic energy has been lost due to the transfer schemes, as the rising bubble lags behind the 1-fluid test case. This is also shown by the energy anomalies in figure 4. As the properties of the two fluids are initially identical and the transfer terms (analytically) only re-label the fluid, we expect the fluid properties to be identical to those of the 1-fluid test case. No temperature or energy anomalies should therefore exist, other than those introduced by numerical errors. The momentum-conserving method 1 schemes (schemes 1-4) perform best and closely to one-another. They produce the smallest temperature and energy anomalies despite the schemes having different boundedness and energy properties. Schemes 5 and 6 from method 2 initially perform similarly to schemes 1-4 but diverge after 150s with larger temperature and energy anomalies. The energy errors are comparable to the worst non-momentum-conserving method 1 schemes. It should be noted that all schemes suffer from some energy increases due to the interpolation errors on the C-grid.
Our simulations suggest that transfer schemes 1-4 give the most desirable numerical properties, however, the test case does not expose potential stability issues such as unbounded velocities/temperatures in schemes 1 and 3, as well as possible energy increases in scheme 3. The simulations show a poor performance of the method 2 schemes, but this is due to the fact that the method uses more interpolated variables relative to method 1 which introduce more interpolation errors. Had we used a co-located grid, the method 2 schemes would be as viable as the method 1 schemes.
Based on our analysis and results, schemes 2 and 4 provide the best overall performance as they possess all of the desired numerical properties which we outlined and the best accuracy in our 2-fluid staggered-grid test case.
We have shown that the multi-fluid shallow water equations with a single height gradient (and without drag or diffusion between fluids) can be stable when the fluid velocities differ. Our results are not inconsistent with those of Stewart and Wendroff (1984) and Thuburn et al. (pre-release: 2019) but provide further insight into the stability of the multi-fluid equations. While the stability of the 2-fluid linearised equations is reduced with respect to the 1-fluid systems, we do at least obtain some physical and numerical stability when we have differing fluid velocities and when using one common pressure gradient, despite no additional mixing terms such as drag and diffusion between the fluids (as were necessary in Weller and McIntyre, 2019). The non-linear equations are more stable than the linearized equations and we have no reason to believe that they are less stable than the 1-fluid non-linear equations. This numerical stability is due to the non-linear advection terms taking energy away from the large-scale modes where instabilities occur. We are of course analysing an idealised system, but the stability properties uncovered here may be applicable to other multi-fluid equation sets.
We have also presented various numerical methods for treating the transfer terms between the fluids. These terms are applicable to any multi-fluid equation set where the mean properties of a fluid are transferred with its mass. We have shown that some of the transfer schemes maintain positive mass, keep prognostic variables bounded, conserve momentum, potential & internal energy and are kinetic energy diminishing. These properties help to keep the overall numerical scheme accurate and stable. However, some of these numerical properties are lost when using a staggered grid due to interpolation errors. Transfer schemes 2 and 4 yielded the most desirable properties on a staggered grid but scheme 6 should also be considered on co-located grids.
More work must be done to understand the stability of the multi-fluid compressible Euler equations. The inherent need for drag or diffusion in multi-fluid schemes (Weller and McIntyre, 2019) means that various arbitrary mixing coefficients are needed to prevent instability. The analysis in this study could perhaps be used to analytically find a drag/diffusion coefficient to ensure total stability of the linearized system which could then be applied to the fully compressible Euler equations for convection modelling. Further studies should also focus on the effects that staggered grids have on the transfer schemes and whether some numerical properties present on co-located transfer schemes can be revived.

Boundedness properties of transfer terms
For a 2-fluid system, bounded velocity transfer terms can be generalised by where β i j ∈ [0, 1] ensures the new velocities do not generate new extrema. Method 1 has β i j = ν M i j | q,r (defined in equation 39). ν M i j | q,r is clearly positive given positive mass and transfer rates. To investigate whether ν M i j | q,r ≤ 1, we make the denominator small (the worst case scenario) such that S j i = 0 s −1 . This gives us This boundedness condition is only guaranteed if α M = 1.
When α M = α C , method 2 has β i j = . This is always true for ∆t S i j ≤ 1, although boundedness is also guaranteed for ∆t S i j > 1 when α C = α M = 1.

Energy properties of transfer terms
For method 2, the total momentum is conserved if the new momenta (F i | n+1 ) satisfy The new kinetic energy after transfers have been applied is where ∆K ≡ 1 2 (u 0 | m − u 1 | m ) . (µ 01 u 0 | m − µ 10 u 1 | m ) and µ i j ≡ λ M i j (1 − β i j − β j i ) + β j i η i | m . When µ 01 = µ 10 and µ 01 ≥ 0, the kinetic energy will never increase. For method 2 (and when α C = α M ), µ i j is given by which is symmetric and always positive meaning this scheme never produces positive energy changes.
Such a proof is less trivial for method 1 as these schemes conserve momentum differently compared to equation 48; η i u j terms are also present in this case. Instead, we have calculated the energy change due to these transfers for a wide range of parameter-space: η i ∈ [10 −16 , 10 5 ], u i ∈ [10 −16 , 10 5 ], S i j ∈ [0, 1], ∆t ∈ [0, 5] -all of which are uniformly discretized 50 times. These 50 7 transfer simulations cover a far wider parameter-space than would be expected in the atmosphere (η i ∼ 1, u i ∼ 0 − 100). Thus if the schemes behave well in the full parameter space, it is reasonable to assume the transfer schemes are useful for convection modelling.
In the 50 7 simulations, scheme 1 was energy-diminishing for ∆t S i j ≤ 1 and schemes 2 & 4 were energy-diminishing for all ∆t S i j > 0. Scheme 3 produced energy increases for ∆t S i j < 1 and could therefore produce numerical instabilities if used.

A C K N O W L E D G E M E N T S
The transfer properties of the transfer schemes from methods 1 and 2. Schemes 2, 5 and 6 have all of the ideal transfer properties. In this study, we have not shown that schemes 1 or 4 always decrease energy but we have observed no energy increases in idealised test cases. Ticks indicate that the scheme fulfils the given property for ∆t S i j ≤ 1, double ticks show that the property also occurs for all ∆t S i j > 0 and crosses mean the property is not fulfilled.

2-fluid cases
F I G U R E 1 The physical and numerical stability regions of the 2-fluid linearized shallow water equations, where (a) shows one-fluid (or quasi-one-fluid) bases and (b-d) show 2-fluid cases. White regions indicate a physically and numerically stable system whereas dark grey regions mean both are unstable. Red regions occur where we have physical stability but numerical instability and the opposite is true in the blue regions. Also included are unstable regions which are within rounding-error of the stability condition (denoted by ∼ unst abl e) and are given by the different shades of red or blue. The stable regions for the numerical simulations are contained within the cyan (linear) and orange (non-linear) contours. L1, L2 and L3 are labels.
a.i) L1 simulation initial conditions b) L1 energy timeseries a.ii) L1 linear simulation, t = 5 s a.iii) L1 non-linear simulation, t = 5 s F I G U R E 3 Snapshots of the simulation initial conditions (a.i) as well as the 5-second profiles of the linear scheme (a.ii) and non-linear scheme (a.iii) for {η 0 = 0.5m, η 1 = 0.5m, c 0 = 0.005, c 1 = 0.1, c g = 0.1} from label L1 in figure 1.b.i. Note that the mass in the linear scheme has become negative and that large-scale modes are beginning to dominate the system, implying numerical instabilities. This is confirmed in the relative energy change timeseries (b) where the linear scheme experiences an exponential increase in energy. Note that the non-linear scheme remains numerically stable for the same initial conditions.