Numerical methods for entrainment and detrainment in the multi‐fluid Euler equations for convection

Convection schemes are a large source of error in global weather and climate models, and modern resolutions are often too fine to parametrize 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 subgrid‐scale convection, including net mass transport by convection and non‐equilibrium dynamics. The technique involves splitting the atmosphere into multiple fluids. For example, the atmosphere could be divided into buoyant updraughts and stable regions. The fluids interact through a common pressure, drag and mass transfers (entrainment and detrainment). Little is known about the numerical properties of mass transfer terms between the fluids. We derive mass transfer terms which relabel the fluids and derive numerical properties of the transfer schemes, including boundedness, momentum conservation and energy conservation on a co‐located grid. Numerical simulations of the multi‐fluid Euler equations using a C‐grid are presented using stable and unstable treatments of the transfers on a well‐resolved two‐fluid dry convection test case. We find two schemes which are conservative, stable and bounded for large time steps, and maintain their numerical properties on staggered grids.


INTRODUCTION
The modelling of atmospheric convection is at the forefront of current meteorological research due to the large errors caused by convection schemes in atmospheric models (e.g., Arakawa, 2004;Yano et al., 2004;Lean et al., 2008;Yano 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 convection schemes assume the scales of convection are small relative to the dynamical flow and that there is no net mass transport due to convection (e.g., Arakawa and Schubert, 1974). This assumption is increasingly unrealistic at finer resolutions (Kwon and Hong, 2017), since mass transport by convection could be larger than other mass fluxes once the size of a convective cell is close to the grid-scale. Convection schemes such as Arakawa and Schubert (1974), Gregory and Rowntree (1990) and Lappen and Randall (2001) ignore net mass transport by convection. Convection schemes also rarely use non-equilibrium dynamics, but instead commonly assume quasi-equilibrium (e.g., Arakawa and Schubert, 1974;Kain and Fritsch, 1990) in which convection is assumed to be in equilibrium with the large-scale forcing. Although some convection schemes incorporate these aspects (such as Gerard and Geleyn, 2005;Kuell and Bott, 2008), few of them offer a consistent treatment of resolved and sub-grid convection. The conditional filtering (or conditional averaging) technique has been proposed for convection modelling due to the possibility of modelling convection over all resolutions, while also representing net mass transport by convection and non-equilibrium dynamics Weller and McIntyre, 2019).
Conditional filtering involves dividing space into various fluids. In the case of convection, one could label fluid 0 as regions of neutrally buoyant air, fluid 1 as convective updraughts and fluid 2 as downdraught regions. One then calculates a convolution with a space-dependent filter, which may be a volume average or a more complicated filter such as a Guassian. Each filtered fluid has its own properties and prognostic variables such as volume fraction, density, temperature and velocity , and the equations of motion are solved for each fluid individually. As the scheme allows for the advection of any fluid to neighbouring 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 is used in other fields of science and engineering (Dopazo, 1977;Baer and Nunziato, 1986;Méchitoua et al., 2003;Guelfi et al., 2007). Lappen and Randall (2001) used the technique to model cumulus convection, using an updraught fluid and a stable fluid. However, the stable fluid in that study was assumed to subside in the same column as the updraught, meaning the scheme does not incorporate net vertical mass transport by convection. More recently, Yano (2014) and Tan et al. (2018) have described how to conditionally filter the anelastic model with the aim of representing subgrid-scale convection, and  have done the same for the fully compressible Euler equations. Subsequent studies have built upon these foundations including  who investigate the conservation properties and normal modes of the equations, Thuburn et al. (2019) who use the method for a two-fluid single-column convective boundary-layer scheme, and Weller and McIntyre (2019) who formulate a numerical solution of the multi-fluid compressible Euler equations. Thus far, little is known about the numerical properties of solutions to the multi-fluid equations. Stewart and Wendroff (1984) and Thuburn et al. (2019) note that the multi-fluid Euler equations are ill-posed when using a single pressure for all fluids and when sub-filter terms are ignored. This property is confirmed by Weller and McIntyre (2019), as drag or mixing between the fluids is necessary to prevent the fluid properties unphysically diverging from each other. However, the numerical properties of transfer terms between fluids has received little attention.
Entrainment and detrainment are key in the interaction of convective clouds and their environment (De Rooy et al., 2013) but many existing mass flux schemes do not conserve momentum in their steady state equations which may cause model inaccuracies (Arakawa, 2004). Transfer terms that exchange mass and other properties between fluids are also crucial for formulating a parametrization of convection using conditional filtering. These transfer terms will be equivalent to entrainment (including cloud-base entrainment) and detrainment which may be adapted from existing frameworks such as Arakawa and Schubert (1974), Betts and Miller (1993), Neggers et al. (2002) and Siebesma et al. (2007). Fluid transfers are given in  and Weller and McIntyre (2019) in terms of transfer rates. Weller and McIntyre (2019) also propose 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 analyse the numerical properties of the transfer terms between fluids for the multi-fluid compressible Euler equations (defined in Section 2). We formulate 20 possible numerical schemes and analyse their properties including conservation, boundedness and stability in Section 3. We then apply the transfer terms to well-resolved two-fluid dry convection test cases in Section 4.

GOVERNING EQUATIONS
The multi-fluid compressible Euler equations are derived in  and we will be using the notation convection from Weller and McIntyre (2019). We have three equations for each fluid including the continuity equation, the potential temperature equation, and momentum equation, where z is the height coordinate and c v = c p ∕ is the heat capacity of dry air at constant volume and is the heat capacity ratio. These energies will be used to assess the numerical stability of the transfer schemes.

FLUID TRANSFER SCHEMES
In convection modelling, entrainment and detrainment are both mass exchanges between the updraught and surrounding environment. It is therefore important that the numerical implementation of these transfer terms for a multi-fluid system have accurate conservation properties and that they do not produce new extrema. For accuracy, mass and momentum should be conserved. For stability, the fluid mass ( i ) must remain positive and velocities should be bounded. For accuracy and stability, potential and internal energy should also be conserved. In addition, the kinetic energy should not increase as resolved kinetic energy decreases when two fluids of differing velocity mix. In reality, any sink of grid-scale energy should be a source of turbulent kinetic energy which, in turn, may be a source of internal energy; as we are not modelling sub-filter-scale variability in this study, we will be ignoring these effects. In this section, we demonstrate the conservation and boundedness properties of the mass transfer terms for the multi-fluid equations and present solutions of the multi-fluid Euler equations with these transfers.

Notation and numerics
In our governing equations, we have assumed that 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 will not form part of this study. For a mean fluid property i ∈ [ i , u i ], the governing equations can be generalised as where F i contains right-hand-side terms such as the pressure gradient term. Applying this to the temperature and momentum equations we get: We will assume that the transfer terms are operator split such that other processes (advection and F i ) act on the prognostic variables first, followed by the transfers: where n is the time-level (t = nΔt), m is the time-level after applying the advection and F i terms, and is the Crank-Nicolson off-centring coefficient. The numerical treatment of the transfer terms will be discussed in Section 3.3. The mass transfers are then based on the most up-to-date states (m) rather than the previous time-level (n). This allows the transfer terms to be independent of the numerical properties of the advection and F i terms. The total momentum should be equal before (m) and after (n + 1) mass transfer such that: The internal energy should also be conserved by the temperature equation transfers: As the transfer terms involve a division by i , we will replace 1∕ i with 1∕max( i , 10 −16 kg ⋅ m −3 ) in any numerical method to avoid issues when one fluid vanishes.

Transferring mass
When transferring mass, we must ensure that mass is conserved and all i remain positive. The mass in each fluid at the end of the timestep ( n+1 i ) is given by the mass after advection ( m i ) plus the discretised transfer term integrated over time Δt: where C determines the numerical treatment of transfer terms in the continuity equation: transfers are conducted explicitly if C = 0 and implicitly if C = 1. Re-arranging these equations for n+1 where and A is a label used to identify the coefficients and .
We use A = C in the continuity equation (mass transfers) and we will later use A = M and A = T for momentum and temperature transfers respectively. The total mass is clearly conserved as ∑ i n+1 i = ∑ i n i and the total potential energy is also conserved. ij is between 0 and 1 for all A when ΔtS ij ≤ 1, meaning 0 and 1 remain positive. When A = 1, any positive ΔtS ij > 0 may be used.

Transferring fluid properties: Method 1
We must also model the transfer of velocity and temperature associated with the re-labelling of mass between fluids described in Section 3.2. The new value of the variable i ∈ [ i , u i ] should be bounded by the old values of fluids i and j (at time level m) so that new extrema are not generated. Also, momentum should be conserved and energy should not increase. Assuming operator-split transfers, the new fluid properties for fluids 0 and 1 are written as where m i are the values after advection. If A = 0, then 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 n+1 where q,r With two degrees of freedom in each of C , A , q and r, a total of 16 different transfer schemes exist using this method. In Appendix A.1, we derive the momentum/internal energy conservation properties for the following four schemes: The other 12 schemes do not conserve momentum/internal energy; we prove this by presenting numerical analysis of the conservation properties. The relative momentum changes (ΔF REL ≡ (F n+1 − F m )∕F 0 ) due to the transfer schemes are calculated using initial conditions which cover a large parameter range, including conditions observed in convective clouds. The transfer schemes were initiated with m ,150] m⋅s −1 and S 10 in the range [0, 1] s −1 , each uniformly discretised 50 times. For a given time step, 1.25 × 10 5 transfers are therefore tested and the range of the relative momentum change for each scheme is plotted. This is shown in Figure 1. These results confirm the momentum conservation analysis of schemes 1-4 (the relative momentum change for these schemes is always zero). The other 12 schemes do not conserve momentum (or internal energy) and will not be analysed further.
Boundedness is also an important property for a stable numerical scheme. will remain bounded if 2), meaning schemes 2 and 4 are bounded. A special case also exists for scheme 1 which ensures values of remain bounded for ΔtS ij ≤ 1. Figure 2a shows the relative energy changes (ΔE REL ≡ (E n+1 − E m )∕E 0 ) of schemes 1-4 over the same range of parameter space used in Figure 1. Schemes 2 (blue) and 4 (black) do not increase the total kinetic energy of the system for any ΔtS ij > 0 and scheme 1 (grey) for ΔtS ij ≤ 1. Scheme 3 (red) may produce large energy increases for any ΔtS ij due to unbounded velocities which cause increases in kinetic energy. The energy analysis comprehensively samples the parameter space which is useful for convection modelling, but these results do not concretely prove that schemes 2 and 4 are always energy diminishing. We should therefore also consider other transfer methods with known energy properties.

Transferring fluid properties: Method 2 (Mass-weighted transfers)
Transfer terms can also be obtained by considering the flux form equations, which are obtained by combining the continuity Equations (1) and (9) with the chain rule. These transfers unconditionally guarantee the conservation of i i (momentum and internal energy) as with the mass transfers seen in Section 3.2. By defining our mass-weighted quantity as Φ i ≡ i i , we get: The equation takes a similar form to Equation (12), meaning we get the solution We have proposed this alternative method as we can demonstrate that the total kinetic energy of the system never increases when C = A (Appendix A.3). In Appendix A.2, we also show that i is bounded when mass and momentum transfers are treated consistently ( C = A ). We will therefore not consider schemes where C ≠ A . Using purely explicit or purely implicit treatments, we therefore have two more viable transfer schemes for the multi-fluid equations: 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 time step. 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. Figure 2b shows the relative energy changes of schemes 5 and 6 over the same parameter space range used for method 1 schemes. The energy changes are consistent with the analysis in Appendix A.3, whereby scheme 5 never increases in energy for ΔtS ij ≤ 1 and scheme 6 for ΔtS ij > 0.

Transfers on a staggered grid
So far, we have assumed that our mass transfers are conducted in the same location, i.e. on a co-located grid (A-grid). But how do the numerical methods change when using a staggered grid? Following the 2D C-grid  (2019), we keep our prognostic mass and temperature defined at cell centres and define our velocities on cell faces. Henceforth, a cell-centred variable ( ) which is linearly interpolated onto cell faces will be denoted by [ ] f and a variable defined on cell faces will be denoted by [ ] c when it is interpolated onto the cell centres. The numerical transfer schemes for the mass and potential temperature remain the same, but some adjustments must be made for the velocity transfers (method 1): , for example.

F I G U R E 2
The range of relative energy changes of (a) schemes 1-4 and (b) schemes 5 and 6. The minimum and maximum energy changes are calculated using the same parameter-space range as Figure 1. Schemes are energy-diminishing if the relative energy change due to the transfer is never above zero, indicated by the horizontal dotted white line. Energy diminishing and energy producing schemes are shown with solid and dashed lines respectively. Scheme 2 (a, blue), scheme 4 (a, black) and scheme 6 (b, black) are energy diminishing. Scheme 1 (a, grey) and scheme 5 (b, grey) also have energy-diminishing properties for ΔtS ij ≤ 1. Although schemes 2 and 4 behave similarly for ΔtS ij ≤ 1, the profiles are not identical. Square-root scales are used for both axes We will use , once again following Weller and McIntyre (2019). The momentum transfers for method 2 become N i is the fluid mass calculated by conducting the mass transfers on the cell faces; this aids in a consistent and accurate conversion of the mass flux to the fluid velocity.
With velocities defined on cell faces, the kinetic energy is calculated on the faces and then interpolated back onto the cell centres: This interpolation method ensures kinetic energy is conserved when converting to the cell centre values (Ringler et al., 2010).

Summary of proposed transfer terms
We have presented six numerical transfer schemes which maintain positivity of mass and conserve mass, momentum, potential energy and internal energy for ΔtS ij ≤ 1. These schemes are presented in Table 1. Schemes 2, 4, 5 and 6 keep the fluid temperatures and velocities bounded, although scheme 5 only does this for ΔtS ij ≤ 1. Only schemes 2, 4 and 6 are kinetic-energy-diminishing for all time steps, meaning schemes 1, 3 and 5 can cause numerical instabilities if ΔtS ij is large. From our analysis, we recommend scheme 4, Momentum and E I conserved?

E K decreases?
Method 1 Note: Schemes 4 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 ΔtS ij ≤ 1, double ticks show that the property also occurs for all ΔtS ij > 0, and crosses mean the property is not fulfilled. Properties which have not been proven, but have been numerically verified (in Figure 2) are shown in brackets.
as they are fully implicit and fulfil all the numerical criteria we have set. Scheme 2 is also a viable scheme if the transfer rate is limited to ΔtS ij ≤ 1 to maintain positive mass. Section 4 will test these schemes on 2D staggered grids.

RISING BUBBLE TEST CASES
In order to test the properties of the various transfer schemes on a staggered grid, we have implemented them into the multi-fluid fully compressible Euler equation solver from Weller and McIntyre (2019) using operator splitting. We will run test cases adapted from the single-fluid rising bubble test case (defined in Bryan and Fritsch, 2002) where an initially stationary temperature anomaly rises and generates resolved circulations ( Figure 3). 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 (where zero-gradient fields are imposed and no fluxes are perpendicular to the boundaries). A uniform potential temperature field of = 300 K is initially chosen with the system in hydrostatic balance and zero velocity. A positive temperature perturbation is then applied at t = 0 s: The perturbation is only applied for L ≤ 1 where x c = 10 km, z c = 2 km and x r = z r = 2 km. For the twofluid experiments the warm anomaly will be applied to 1 only, whereas fluid 0 will remain initialised as 0 = 300 K. We use a 2D C-grid with i and i defined at cell centres and the normal component of u i defined at cell faces. The time-stepping is centred Crank-Nicolson with a time step of Δt = 2 s. A van Leer advection scheme is chosen to maintain positivity of the mass of each fluid. All details of the numerical set-up and the numerical solvers used are described in Weller and McIntyre (2019), with the exception of a numerical adjustment which must be made for an operator-split Crank-Nicolson multi-fluid scheme (described in Appendix B).

Full bubble test case
The first test case is initialised with all mass in fluid 1: 0 = 0, 1 = 1. The transfer rate is chosen to transfer a large quantity of fluid 1 to fluid 0: where min = 0.1. This means that the explicit schemes will transfer 10% of the mass in the first time step (and none thereafter). As fluid 0 initially has no mass, it should inherit the properties of fluid 1 when mass is transferred. We therefore expect the solution to be the same as the single-fluid test case shown in Figure 3. The test case is run for all 20 transfer schemes, including the non-conservative schemes. For each scheme we calculate the relative energy change from the single fluid test case: where E n SF and E n MF are the total energies at time step n for the single-fluid and multi-fluid simulations respectively. With a float precision up to 16 decimal places, we expect fluid 0 to inherit the density, velocity and temperature of fluid 1 to machine precision. A relative energy change of ΔE n RSF ∼ 10 −15 is therefore reasonable for an energy-conserving scheme.
The energy changes for all schemes are shown in Table 2. All six conservative schemes produce small energy decreases after one time step with the implicit mass schemes (3, 4 and 6) producing the smallest energy changes with ΔE 1 RSF = −1.18 × 10 −15 as they have not transferred the full 10% of the mass in the first time step (unlike schemes 1, 2 and 5). The relative energy changes of schemes 1-6 remain of the order 10 −15 by t = 1000 s with schemes 5 and 6 having the smallest errors. Many of the non-conservative schemes produce large energy increases due to lack of internal energy conservation and unbounded velocities. Some of these schemes become unstable before the end of the test case at t = 1000 s. Note that two of the non-conservative schemes behave similarly to schemes 1-6, but internal energy and momentum are not conserved exactly.

Half-bubble test case
We have already shown solutions for transfers to an empty fluid. But how do the schemes behave when transferring between fluids with comparable mass and different properties? For this we use a two-fluid test case from Weller and McIntyre (2019), where half the mass is initialised with the warm anomaly (fluid 1) and the other half without (fluid 0):

F I G U R E 3
The temperature profile evolution of the full bubble test case, which has the same analytical solution as the single-fluid test case from Bryan and Fritsch (2002). The warm anomaly rises and induces large-scale resolved circulations. The black arrows give the relative magnitudes and directions of the velocity vectors TA B L E 2 The relative energy changes of the two-fluid rising bubble test case, relative to the single-fluid test case  The two-fluid equations with different fluid properties require some form of stabilisation (Thuburn et al., 2019;Weller and McIntyre, 2019). We will use a diffusive mass transfer to couple the fluids: where K = 200 m 2 ⋅s −1 is a large-enough diffusion coefficient to maintain numerical stability for this test case (as shown in Weller and McIntyre, 2019). The temperature and volume fraction distributions for this test case are shown in Figure 4; slower circulations form than in the full bubble test case due to the lower mean temperature anomaly. The energy changes of all schemes relative to the initial conditions are shown in Figure 5. Dashed lines represent negative energy changes and solid lines show positive energy changes. Schemes 1-6 follow similar energy evolutions, where energy decreases relative to the initial conditions. The non-conservative schemes (light grey) exhibit various behaviours; many blow up within the first time steps and produce large energy increases, whereas some schemes (which use implicit transfers) follow similar energy evolutions to the conservative schemes. Note that energy changes due to other aspects of the numerical scheme (e.g. the non-energy-conserving van Leer advection scheme) are far larger than changes due to the transfer schemes; hence all conservative schemes appear to behave similarly.
These simulations on a staggered grid are consistent with the analysis of the transfer schemes in Section 3 as schemes 1-6 conserve mass, momentum, internal energy and potential energy and are stable for the given test cases with no positive increases in total energy.

CONCLUSIONS
Transfer terms between the fluid components in the multi-fluid equations can be used to couple the fluids, represent physical exchanges and stabilise the equations, but the numerical treatment of the transfer terms must F I G U R E 4 The (a) temperature and (b) volume fraction profiles of the half-bubble test case using scheme 2. The black arrows give the relative magnitudes and directions of the velocity vectors. As the warm anomaly is initially only in half the fluid, the distributions differ from the single-fluid case in Figure 3, including a slower circulation over the domain F I G U R E 5 The relative energy change from the initial energy for all of the transfer schemes for the two-fluid rising bubble test case. The solid and dashed lines represent positive and negative energy differences respectively. Dotted lines represent the transition between positive and negative energy changes. Schemes 1-6 are given by the black lines, while the remaining non-conservative schemes are given by the light grey lines be stable. We have presented various numerical methods for treating the transfer terms between the fluids. These schemes 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 and internal energy, and are kinetic energy diminishing. These properties help to keep the overall numerical scheme accurate and stable on co-located and staggered grids. Of the six conservative schemes, we have shown that scheme 3 can produce energy increases (Figure 2). We have also shown that schemes 5 and 6 do not increase energy for ΔtS ij ≤ 1 and ΔtS ij > 0 respectively. We have not proved this for schemes 1, 2 and 4, but have thoroughly explored the relevant parameter space for convection and have not found any instances of kinetic energy increases (other than scheme 1 for ΔtS ij > 1). The fully implicit schemes (schemes 4 and 6) automatically handle large mass transfers and scheme 6 also produces the smallest energy changes in the full bubble test case. As the energy properties are exactly known, scheme 6 has the most desirable numerical properties. Our numerical analysis suggests that scheme 4 is also an ideal transfer scheme and that scheme 2 is useful if the transfer is limited so that mass positivity is maintained. By using any of the fully implicit schemes, the entrainment and detrainment in the multi-fluid equations can be conducted in a numerically stable manner. The physical form of the entrainment and detrainment transfer terms should be the focus of future studies so that convective processes can be accurately represented in the multi-fluid equations.

A. NUMERICAL PROPERTIES OF THE TRANSFER SCHEMES
Here we derive the conservation, boundedness and energy properties of the transfer schemes which are important for a stable numerical method.

A.1. Conservation properties
By multiplying Equation (13) by Equation (16) and summing over both fluid components, we get the total mass-weighted properties: For conservation of internal energy/momentum, we require that N 01 m 0 + N 10 m 1 = 0. As m 0 and m 1 are independent, we additionally impose that N ij = 0 for conservation.
Setting A = 0 and r = n + 1: By using A = 0 and r = n + 1, Setting C = 0 and q = m means that Cij =ΔtS ij which results in N ij = 0, meaning scheme 1 is conservative. If instead we use C = 1 and q = n + 1 (scheme 3), we get Setting A = 1 and r = m: By instead using A = 1 and r = m (and expanding the n+1 i terms in Equation (A2) to get m i terms), we get The terms in the parentheses are the same as Equation (A3). This means that N ij = 0 if C = 0 and q = m (scheme 2) or if C = 1 and q = n + 1 (scheme 4).

A.2. Boundedness properties
For a two-fluid system, bounded velocity transfer terms can be generalised by where 0 ≤ ij ≤ 1 ensures the new velocities are bounded. Method 1 has ij = q,r Mij (defined in Equation (17)).

q,r
Mij is clearly positive given positive mass and transfer rates. To investigate whether q,r Mij ≤ 1, we make the denominator small (the worst case scenario) such that S ji = 0s −1 . This gives j } which is bounded if 0 ≤ Cij ≤ 1. This is always true for ΔtS ij ≤ 1, although boundedness is also guaranteed for ΔtS ij > 1 when C = M = 1.

A.3. Energy properties
For method 2, the total momentum is conserved if the new momenta (F n+1 i ) satisfy where F m i ≡ m i u m i . The new kinetic energy after transfers have been applied is found by multiplying Equations (A6) and (A8) and summing over all fluids: When 01 = 10 and 01 ≥ 0, the kinetic energy will never increase. For method 2 (and when C = M ), ij 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 than Equation (A8); 0 u 1 and 1 u 0 terms are also present with method 1. Instead, the energy changes are calculated over a large range of parameter space and are shown in Figure 2a.

B. NUMERICAL ADJUSTMENTS FOR AN OPERATOR-SPLIT CRANK-NICOLSON MULTI-FLUID SCHEME
The numerical multi-fluid scheme used for this study follows the implementation by Weller and McIntyre (2019), with the exception of operator-split transfers. For a fluid property such as temperature or velocity ( ), the solution for the Crank-Nicolson scheme (before transfers) is given by where is the off-centring coefficient and the i ∕ t terms represent the advection and F in Equation (10). As ( i ∕ t) n is stored from the previous time step, we must ensure that it remains consistent with the fluid properties when transfers are made. This is done by computing for method 1 schemes and for method 2 schemes. ( i ∕ t) n+1 is therefore stored for the next time step (and is not used to calculate n+1 i ). Absence of these terms leads to errors in the numerical solution when using operator-split transfers, especially when a fluid has a small volume fraction or if large transfers are conducted. These terms are not necessary if the Crank-Nicolson off-centring coefficient is set to = 1, but the scheme will be limited to first-order accuracy in time.