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

TL;DR

This paper investigates numerical methods for modeling convection in climate models using multi-fluid Euler equations, focusing on entrainment and detrainment schemes that improve stability and conservation properties.

Contribution

It derives and analyzes mass transfer schemes between fluids in multi-fluid Euler models, demonstrating stable, conservative, and bounded numerical solutions for convection.

Findings

Identified schemes that are conservative, stable, and bounded for large timesteps.

Demonstrated numerical properties on a well-resolved convection test case.

Maintained properties on staggered grids.

Abstract

Convection schemes are a large source of error in global weather and climate models, and modern resolutions are often too fine to parameterise 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 multiple fluids. For example, the atmosphere could be divided into buoyant updrafts 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. Numerical solutions 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 timesteps, and maintain their numerical properties on staggered grids.

Figures

Tables

Table 1. Table 1: The transfer properties of the transfer schemes from methods 1 and 2. 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 Δ ​ t ​ S i ​ j ≤ 1 Δ 𝑡 subscript 𝑆 𝑖 𝑗 1 \Delta tS_{ij}\leq 1 , double ticks show that the property also occurs for all Δ ​ t ​ S i ​ j > 0 Δ 𝑡 subscript 𝑆 𝑖 𝑗 0 \Delta tS_{ij}>0 and crosses mean the property is not fulfilled.

	Scheme	Positive ${\eta }_i$?	Bounded ${\theta }_i$ & $u_i$?	\LongunderstackMomentum &
IE conserved?	KE decreases?
	Method 1
1	${\alpha }_C=0$, ${\alpha }_A=0$, $q=m$, $r=n+1$	✓	✗	✓✓	✓
2	${\alpha }_C=0$, ${\alpha }_A=1$, $q=m$, $r=m$	✓	✓✓	✓✓	✓✓
3	${\alpha }_C=1$, ${\alpha }_A=0$, $q=n+1$, $r=n+1$	✓✓	✗	✓✓	✗
4	${\alpha }_C=1$, ${\alpha }_A=1$, $q=n+1$, $r=m$	✓✓	✓✓	✓✓	✓✓
	Other schemes	✓	✓ for ${\alpha }_A=1$	✗	✗
	Method 2 (Mass-weighted transfers)
5	${\alpha }_C=0$, ${\alpha }_A=0$	✓	✓	✓✓	✓
6	${\alpha }_C=1$, ${\alpha }_A=1$	✓✓	✓✓	✓✓	✓✓
	Other schemes	✓	✓ for ${\alpha }_C={\alpha }_A$	✓✓	✓ for ${\alpha }_C={\alpha }_A$

Table 2. Table 2: The relative energy changes of the 2-fluid rising bubble test case, relative to the single-fluid test case. Energy changes are shown for the first timestep ( n = 1 𝑛 1 n=1 ) where the transfers are largest and after 1000 1000 1000~{} s ( n = 500 𝑛 500 n=500 ). Given the machine precision, a relative energy change of the order 10 − 15 superscript 10 15 10^{-15} is expected for energy conservation - results around this range are indicated in bold.

Name	Method	${\alpha }_C$	${\alpha }_A$	$q$	$r$	$\Delta E_{RSF}^1$	$\Delta E_{RSF}^{500}$
Scheme 1	1	$0$	$0$	$m$	$n+1$	$-1.18\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟏𝟒}}$	$-4.71\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟏𝟓}}$
$-$	1	$0$	$1$	$m$	$n+1$	$+6.77\times {10}^{-7}$	$-1.34\times {10}^{-6}$
$-$	1	$1$	$0$	$m$	$n+1$	$-1.23\times {10}^{-7}$	$-1.97\times {10}^{-6}$
$-$	1	$1$	$1$	$m$	$n+1$	$+5.85\times {10}^{-7}$	$-2.35\times {10}^{-6}$
$-$	1	$0$	$0$	$m$	$m$	$+1.05\times {10}^{11}$	$-$
Scheme 2	1	$0$	$1$	$m$	$m$	$-1.18\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟏𝟒}}$	$-7.07\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟏𝟓}}$
$-$	1	$1$	$0$	$m$	$m$	$+9.51\times {10}^{10}$	$-$
$-$	1	$1$	$1$	$m$	$m$	$-1.18\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟏𝟓}}$	$-6.73\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟏𝟓}}$
$-$	1	$0$	$0$	$n+1$	$n+1$	$+1.35\times {10}^{-7}$	$-6.53\times {10}^{-6}$
$-$	1	$0$	$1$	$n+1$	$n+1$	$+7.12\times {10}^{-7}$	$-1.07\times {10}^{-6}$
Scheme 3	1	$1$	$0$	$n+1$	$n+1$	$-1.18\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟏𝟓}}$	$-6.56\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟏𝟓}}$
$-$	1	$1$	$1$	$n+1$	$n+1$	$+6.15\times {10}^{-7}$	$-2.00\times {10}^{-6}$
$-$	1	$0$	$0$	$n+1$	$m$	$+8.47\times {10}^{10}$	$-$
$-$	1	$0$	$1$	$n+1$	$m$	$-1.18\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟏𝟒}}$	$-4.71\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟏𝟓}}$
$-$	1	$1$	$0$	$n+1$	$m$	$+7.86\times {10}^{10}$	$-$
Scheme 4	1	$1$	$1$	$n+1$	$m$	$-1.18\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟏𝟓}}$	$-8.75\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟏𝟓}}$
Scheme 5	2	$0$	$0$	$-$	$-$	$-1.18\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟏𝟒}}$	$-3.87\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟏𝟓}}$
$-$	2	$0$	$1$	$-$	$-$	$+5.13\times {10}^{-14}$	$-$
$-$	2	$1$	$0$	$-$	$-$	$+4.58\times {10}^{-14}$	$-$
Scheme 6	2	$1$	$1$	$-$	$-$	$-1.18\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟏𝟓}}$	$-3.19\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟏𝟓}}$

Equations

$$\frac{\partial {\eta }_i}{\partial t}+\nabla .({\eta }_i{\mathbf{u}}_i)=\underset{\text{Mass transfers }}{\underbrace{\sum _{j\ne i}({\eta }_j{S}_{ji}-{\eta }_i{S}_{ij})}},$$

$$\frac{\partial {\theta }_i}{\partial t}+{\mathbf{u}}_i.\nabla {\theta }_i=\underset{\text{Transfer mean temp. }}{\underbrace{\sum _{j\ne i}\left(\frac{{\eta }_j}{{\eta }_i}{S}_{ji}({\theta }_j-{\theta }_i)\right)}}-\underset{\text{Heat transfer }}{\underbrace{\sum _{j\ne i}{H}_{ij}}},$$

$$\frac{\partial {\mathbf{u}}_i}{\partial t}+{\mathbf{u}}_i.\nabla {\mathbf{u}}_i=\mathbf{g}-\underset{\text{Pressure gradient }}{\underbrace{c_p{\theta }_i\nabla \pi }}+\underset{\text{Transfer mean velocity. }}{\underbrace{\sum _{j\ne i}\left(\frac{{\eta }_j}{{\eta }_i}{S}_{ji}({\mathbf{u}}_j-{\mathbf{u}}_i)\right)}}-\underset{\text{Drag }}{\underbrace{\sum _{j\ne i}{\mathbf{D}}_{ij}}},$$

$$p_0{\pi }^{\frac{1-\kappa }{\kappa }}=R\sum _i{\eta }_i{\theta }_i.$$

$$E=E_P+E_I+E_K$$

$$E_P$$

$$E_I$$

$$E_K$$

$$\frac{\partial {\varphi }_i}{\partial t}+{\mathbf{u}}_i.\nabla {\varphi }_i=F_i+\sum _{j\ne i}{S}_{ji}\frac{{\eta }_j}{{\eta }_i}\left[{\varphi }_j-{\varphi }_i\right],$$

$${\varphi }_i^m$$

$${\varphi }_i^{n+1}$$

$$\begin{array}{rl}{\eta }_0^{n+1}& ={\eta }_0^m-\Delta t\left[(1-{\alpha }_C){\eta }_0^m+{\alpha }_C{\eta }_0^{n+1}\right]S_{01}+\Delta t\left[(1-{\alpha }_C){\eta }_1^m+{\alpha }_C{\eta }_1^{n+1}\right]S_{10},\\ {\eta }_1^{n+1}& ={\eta }_1^m-\Delta t\left[(1-{\alpha }_C){\eta }_1^m+{\alpha }_C{\eta }_1^{n+1}\right]S_{10}+\Delta t\left[(1-{\alpha }_C){\eta }_0^m+{\alpha }_C{\eta }_0^{n+1}\right]S_{01},\end{array}$$

$$\begin{array}{rl}{\eta }_0^{n+1}& =(1-{\lambda }_{C01}){\eta }_0^m+{\lambda }_{C10}{\eta }_1^m,\\ {\eta }_1^{n+1}& =(1-{\lambda }_{C10}){\eta }_1^m+{\lambda }_{C01}{\eta }_0^m,\end{array}$$

$$\begin{array}{rl}{\lambda }_{Aij}& \equiv \frac{\Delta t{S}_{ij}}{1+{\alpha }_A\Delta t(S_{ji}+S_{ij})}\end{array}$$

$$\begin{array}{rl}{\varphi }_0^{n+1}& ={\varphi }_0^m-(1-{\alpha }_A)\Delta t\frac{{\eta }_1^q}{{\eta }_0^r}{S}_{10}({\varphi }_0^m-{\varphi }_1^m)-{\alpha }_A\Delta t\frac{{\eta }_1^q}{{\eta }_0^r}{S}_{10}({\varphi }_0^{n+1}-{\varphi }_1^{n+1}),\\ {\varphi }_1^{n+1}& ={\varphi }_1^m-(1-{\alpha }_A)\Delta t\frac{{\eta }_0^q}{{\eta }_1^r}{S}_{01}({\varphi }_1^m-{\varphi }_0^m)-{\alpha }_A\Delta t\frac{{\eta }_0^q}{{\eta }_1^r}{S}_{01}({\varphi }_1^{n+1}-{\varphi }_0^{n+1}),\end{array}$$

$$\begin{array}{rl}{\varphi }_0^{n+1}& =(1-{\nu }_{A10}^{q,r}){\varphi }_0^m+{\nu }_{A10}^{q,r}{\varphi }_1^m,\\ {\varphi }_1^{n+1}& =(1-{\nu }_{A01}^{q,r}){\varphi }_1^m+{\nu }_{A01}^{q,r}{\varphi }_0^m,\end{array}$$

$${\nu }_{Aij}^{q,r}=\frac{\Delta t{S}_{ij}\frac{{\eta }_i^q}{{\eta }_j^r}}{1+{\alpha }_A\Delta t\left[S_{ij}\frac{{\eta }_i^q}{{\eta }_j^r}+S_{ji}\frac{{\eta }_j^q}{{\eta }_i^r}\right]}.$$

$$\frac{\partial ({\eta }_i{\varphi }_i)}{\partial t}+\nabla .({\eta }_i{\varphi }_i{\mathbf{u}}_i)={\eta }_i{F}_i+\sum _j\left[S_{ji}{\eta }_j{\varphi }_j-S_{ij}{\eta }_i{\varphi }_i\right],$$

$$\begin{array}{rl}{\Phi }_0^{n+1}& ={\Phi }_0^m-\Delta t\left[(1-{\alpha }_A){\Phi }_0^m+{\alpha }_A{\Phi }_0^{n+1}\right]S_{01}+\Delta t\left[(1-{\alpha }_A){\Phi }_1^m+{\alpha }_A{\Phi }_1^{n+1}\right]S_{10},\\ {\Phi }_1^{n+1}& ={\Phi }_1^m-\Delta t\left[(1-{\alpha }_A){\Phi }_1^m+{\alpha }_A{\Phi }_1^{n+1}\right]S_{10}+\Delta t\left[(1-{\alpha }_A){\Phi }_0^m+{\alpha }_A{\Phi }_0^{n+1}\right]S_{01}.\end{array}$$

$$\begin{array}{rl}{\Phi }_i^m& \equiv {\eta }_i^m{\varphi }_i^m\\ {\Phi }_0^{n+1}& =(1-{\lambda }_{A01}){\Phi }_0^m+{\lambda }_{A10}{\Phi }_1^m,\\ {\Phi }_1^{n+1}& =(1-{\lambda }_{A10}){\Phi }_1^m+{\lambda }_{A01}{\Phi }_0^m,\\ {\varphi }_i^{n+1}& =\frac{{\Phi }_i^{n+1}}{{\eta }_i^{n+1}}.\end{array}$$

$$\begin{array}{rl}{w}_0^{n+1}& =\left(1-[{\nu }_{M10}^{q,r}{]}_f\right)w_0^m+[{\nu }_{M10}^{q,r}{]}_f{w}_1^m,\\ w_1^{n+1}& =\left(1-[{\nu }_{M01}^{q,r}{]}_f\right)w_1^m+[{\nu }_{M01}^{q,r}{]}_f{w}_0^m,\end{array}$$

$$\begin{array}{rl}{N}_0^{n+1}& \equiv \left(1-[{\lambda }_{C01}{]}_f\right)[{\eta }_0^m{]}_f+[{\lambda }_{C10}{]}_f[{\eta }_1^m{]}_f,\\ N_1^{n+1}& \equiv \left(1-[{\lambda }_{C10}{]}_f\right)[{\eta }_1^m{]}_f+[{\lambda }_{C01}{]}_f[{\eta }_0^m{]}_f,\\ {\mathbf{F}}_i^m& \equiv [{\eta }_i^m{]}_f{\mathbf{u}}_i^m\\ {\mathbf{F}}_0^{n+1}& =\left(1-[{\lambda }_{M01}{]}_f\right){\mathbf{F}}_0^m+[{\lambda }_{M10}{]}_f{\mathbf{F}}_1^m,\\ {\mathbf{F}}_1^{n+1}& =\left(1-[{\lambda }_{M10}{]}_f\right){\mathbf{F}}_1^m+[{\lambda }_{M01}{]}_f{\mathbf{F}}_0^m.\\ {\mathbf{u}}_i^{n+1}& =\frac{{\mathbf{F}}_i^{n+1}}{N_i^{n+1}},\end{array}$$

$$E_K^{n+1}=\sum _i\frac{1}{2}{\left[N_i^{n+1}(w_i^{n+1}{)}^2\right]}_c.$$

$${\theta }^{\prime }=2{\cos }^2\left(\frac{\pi }{2}L\right).$$

$$S_{01}$$

$$S_{10}$$

$$\Delta E_{RSF}^n=\frac{E_{MF}^n-E_{SF}^n}{E_{SF}^0},$$

$$\begin{array}{rl}{\sigma }_1& =\{\begin{array}{ll}0.5& \text{for }L<1,\\ 0& \text{otherwise},\end{array}\\ {\sigma }_0& =1-{\sigma }_1.\end{array}$$

$$S_{ij}=\frac{1}{2}\frac{K_{\sigma }}{{\eta }_i}\text{max}\left(0,{\nabla }^2({\eta }_j-{\eta }_i)\right),$$

$$\begin{array}{rl}{\mathbf{u}}_0^{n+1}& =(1-{\beta }_{10}){\mathbf{u}}_0^m+{\beta }_{10}{\mathbf{u}}_1^m,\\ {\mathbf{u}}_1^{n+1}& =(1-{\beta }_{01}){\mathbf{u}}_1^m+{\beta }_{01}{\mathbf{u}}_0^m,\end{array}$$

$$\begin{array}{rl}{\nu }_{Mij}^{q,r}=\frac{\Delta t{S}_{ij}\frac{{\eta }_i^q}{{\eta }_j^r}}{1+{\alpha }_M\Delta t{S}_{ij}\frac{{\eta }_i^q}{{\eta }_j^r}}& \le 1\\ \Delta t{S}_{ij}\frac{{\eta }_i^q}{{\eta }_j^r}& \le \frac{1}{1-{\alpha }_M}.\end{array}$$

$$\begin{array}{rl}{\mathbf{F}}_0^{n+1}& =(1-{\lambda }_{M01}){\mathbf{F}}_0^m+{\lambda }_{M10}{\mathbf{F}}_1^m,\\ {\mathbf{F}}_1^{n+1}& =(1-{\lambda }_{M10}){\mathbf{F}}_1^m+{\lambda }_{M01}{\mathbf{F}}_0^m,\end{array}$$

Full text

1 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 [eg 1, 15, 26, 27]. We have reached the “grey zone” in which improved computational power allows convection to be partially resolved but cannot yet be explicitly simulated [8]. 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 [eg 2]. This assumption is increasingly unrealistic at finer resolutions [13], since mass transport by convection could be larger than other mass fluxes once a convective cell is close to the grid-scale. This possibility is commonly ignored in convection schemes [such as 2, 9, 14] and so is non-equilibrium dynamics [for example 11]. Although many convection schemes incorporate some of these aspects [such as 7, 12], 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, whilst also representing net mass transport by convection and non-equilibrium dynamics [24, 25].

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 updrafts and fluid 2 as downdraft 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 [24], 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 [6, 3, 16, 10]. [14] 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 the same column as the updraft, meaning the scheme does not incorporate net vertical mass transport by convection. More recently, [24] and [21] 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 [23] who investigate the conservation properties and normal modes of the equations, [22] who use the method for a 2-fluid single-column convective boundary layer scheme, and [25] 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. [20] and [22] note that the multi-fluid Euler equations are ill-posed when sub-filter terms are ignored. This property is confirmed by [25], 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.

Transfer terms that exchange mass and other properties between fluids are crucial for formulating a parameterisation 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 [2], [4], [17] and [19]. Fluid transfer terms are given in [23] and [25] in terms of transfer rates. [25] 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.

2 Governing equations

The multi-fluid compressible Euler equations are derived in [24] and we will be using the notation convection from [25]. We have three equations for each fluid including the continuity equation,

[TABLE]

the potential temperature equation,

[TABLE]

and momentum equation,

[TABLE]

where $i$ is the label for fluid $i$. Our prognostic variables are the fluid mass per unit volume ($\eta_{i}$), the fluid potential temperature ($\theta_{i}$) and the fluid velocity ($\bm{u}_{i}$). We have defined $\eta_{i}\equiv\sigma_{i}\rho_{i}$, where $\sigma_{i}$ is the fluid volume fraction and $\rho_{i}$ is the density of fluid $i$. $S_{ij}$ is the unidirectional mass transfer rate from fluid $i$ to fluid $j$ ($S_{ij}\geq 0$). $c_{p}$ is the heat capacity of dry air at constant pressure and $\bm{g}$ is the gravitational acceleration. $\pi\equiv p^{\kappa}/p_{0}^{\kappa}$ is the Exner pressure where $p$ is the pressure, $p_{0}$ is a reference pressure, $\kappa=R/c_{p}$ and $R$ is the gas constant of dry air. $H_{ij}$ is the heat transfer between fluids $i$ and $j$ (which does not exchange mass between fluids) and $\bm{D}_{ij}$ is the drag between fluids $i$ and $j$ - these exchange terms will not be used in this study. Additionally, the equation of state for dry air is used to relate the pressure and fluid temperatures:

[TABLE]

The total energy of the multi-fluid system is given by

[TABLE]

where $E_{P}$ is the potential energy, $E_{I}$ is the internal energy and $E_{K}$ is the kinetic energy, defined respectively as:

[TABLE]

where $z$ is the height coordinate and $c_{v}=\frac{c_{p}}{\gamma}$ is the heat capacity of dry air at constant volume and $\gamma$ is the heat capacity ratio. These energies will be used to assess the numerical stability of the transfer schemes.

3 Fluid transfer schemes

In convection modelling, entrainment and detrainment are both mass exchanges between the updraft 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 ($\eta_{i}$) must remain positive and velocities should be bounded. For accuracy and stability, potential and internal energy should also be conserved and the kinetic energy should not increase (resolved kinetic energy decreases when two fluids of differing velocity mix). 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.

3.1 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 $\phi_{i}\in[\theta_{i},\bm{u}_{i}]$, the governing equations can be generalised as

[TABLE]

where $F_{i}$ contains right-hand-side terms such as the pressure gradient term. Applying this to the temperature and momentum equations we get:

• •

Momentum equation: $\bm{\phi}_{i}=\bm{u}_{i}$, $\bm{F}_{\bm{u}i}=-c_{p}\theta_{i}\bm{\nabla}\pi+\bm{g}$.

• •

Temperature equation: $\phi_{i}=\theta_{i}$, $\bm{F}_{\theta i}=0$.

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:

[TABLE]

where $n$ is the time-level ($t=n\Delta t$), $m$ is the time-level after applying the advection and $F_{i}$ terms and $\alpha$ is the Crank-Nicolson off-centering coefficient. 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. For each momentum equation, the total momentum should be equal before ($m$) and after ($n+1$) mass transfer such that: $\sum_{i}\eta_{i}^{n+1}\bm{u}_{i}^{n+1}=\sum_{i}\eta_{i}^{m}\bm{u}_{i}^{m}$. The internal energy should also be conserved by the temperature equation transfers: $c_{v}\pi\sum_{i}\eta_{i}^{n+1}\theta_{i}^{n+1}=c_{v}\pi\sum_{i}\eta_{i}^{m}\theta_{i}^{m}$.

3.2 Transferring mass

When transferring mass, we must ensure that mass is conserved and all $\eta_{i}$ remain positive. The mass in each fluid at the end of the timestep ($\eta_{i}^{n+1}$) is given by the mass after advection ($\eta_{i}^{m}$) plus the discretised transfer term integrated over time $\Delta t$:

[TABLE]

where $\alpha_{C}$ determines whether the transfer terms in the continuity equation are numerically treated explicitly ($\alpha_{C}=0$) or implicitly ($\alpha_{C}=1$). Re-arranging these equations for $\eta_{i}^{n+1}$, we get

[TABLE]

where

[TABLE]

and $A$ is a label used to identify the transfer coefficients. 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 $\sum_{i}\eta_{i}^{n+1}=\sum_{i}\eta_{i}^{n}$ and the total potential energy is also conserved. $\lambda_{ij}$ is between 0 and 1 for all $\alpha_{A}$ when $\Delta tS_{ij}\leq 1$, meaning $\eta_{0}$ and $\eta_{1}$ remain positive. When $\alpha_{A}=1$, any positive $\Delta tS_{ij}>0$ may be used.

3.3 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 $\phi_{i}\in[\theta_{i},\bm{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 [math] and $1$ are written as

[TABLE]

where $\phi_{i}^{m}$ are the values after advection. If $\alpha_{A}=0$, then $\phi_{i}$ is treated explicitly and $\alpha_{A}=1$ means $\phi_{i}$ is treated implicitly. Note that these equations have additional degrees of freedom in the time-level choice for $\eta_{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 $\phi_{i}^{n+1}$, we obtain

[TABLE]

where

[TABLE]

With two degrees of freedom in each of $\alpha_{C}$, $\alpha_{A}$, $q$ and $r$, a total of 16 different transfer schemes exist using this method. It is trivial to derive the conservation of momentum/internal energy for the following four schemes:

$\alpha_{C}=0$, $\alpha_{M}=\alpha_{T}=0$ with $q=m$, $r=n+1$. 2. 2.

$\alpha_{C}=0$, $\alpha_{M}=\alpha_{T}=1$ with $q=m$, $r=m$. 3. 3.

$\alpha_{C}=1$, $\alpha_{M}=\alpha_{T}=0$ with $q=n+1$, $r=n+1$. 4. 4.

$\alpha_{C}=1$, $\alpha_{M}=\alpha_{T}=1$ with $q=n+1$, $r=m$.

The other 12 schemes do not conserve momentum/internal energy - proof of this is intractable so we instead present numerical analysis of the conservation properties. The relative momentum changes ($\Delta F_{\text{REL}}\equiv(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 $\eta_{0}^{m}=1~{}$kg m*-3*, $u_{0}^{m}=1~{}$ms*-1*, $\theta_{0}^{m}=300~{}$K, $\theta_{1}^{m}=301~{}$K, $S_{01}=1~{}$s*-1*. We also use $\Delta t$ in the range $[0,5]~{}$s, $\eta_{1}^{m}$ in the range $[10^{-8},2]~{}$kg m*-3*, $u_{1}^{m}$ in the range $[-150,150]~{}$ms*-1* and $S_{10}$ in the range $[0,1]~{}$s*-1* - each uniformly discretised 50 times. For a given timestep, $1.25\times 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.

For schemes 1-4, $\phi$ will remain bounded if $\nu_{ij}^{q,r}\in[0,1]$, but this can only be guaranteed when $\alpha_{A}=1$ (we prove this in appendix section 6.1) meaning schemes 1 and 3 can produce unbounded velocities and temperatures. Figure 2.a shows the relative energy changes ($\Delta E_{\text{REL}}\equiv(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 $\Delta tS_{ij}>0$ and scheme 1 (grey) for $\Delta tS_{ij}\leq 1$. Scheme 3 (red) may produce large energy increases for any $\Delta tS_{ij}$. 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 schemes with known energy properties.

3.4 Transferring fluid properties - Method 2 (Mass-weighted transfers)

Transfer terms can also be obtained by considering the flux form equations,

[TABLE]

which are obtained by combining the continuity equation 1 and equation 9 with the chain rule. These transfers unconditionally guarantee the conservation of $\eta_{i}\phi_{i}$ (momentum and internal energy) as with the mass transfers seen in section 3.2. By defining our mass-weighted quantity as $\Phi_{i}\equiv\eta_{i}\phi_{i}$, we get:

[TABLE]

The equation takes a similar form to (12), meaning we get the solution

[TABLE]

We have proposed this alternative method as we can demonstrate that the total kinetic energy of the system never increases when $\alpha_{C}=\alpha_{A}$ (see appendix 6.2). In appendix 6.1, we also show that $\phi_{i}$ is bounded when mass and momentum transfers are treated consistently ($\alpha_{C}=\alpha_{A}$) - we will therefore not consider schemes where $\alpha_{C}\neq\alpha_{A}$. Using purely explicit or purely implicit treatments, we therefore have 2 more viable transfer schemes for the multi-fluid equations:

$\alpha_{C}=0$, $\alpha_{A}=\alpha_{M}=\alpha_{T}=0$. 2. 6.

$\alpha_{C}=1$, $\alpha_{A}=\alpha_{M}=\alpha_{T}=1$.

Note that schemes such as $[\alpha_{C}=0.5,\ \alpha_{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. Figure 2.b 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 6.2, whereby scheme 5 never increases in energy for $\Delta tS_{ij}\leq 1$ and scheme 6 for $\Delta tS_{ij}>0$.

3.5 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 staggered grid? Following the C-grid setup used in [25], we keep our prognostic mass and temperature defined at cell centres and define our velocities on cell faces. Henceforth, a cell-centred variable ($\mu$) which is linearly-interpolated onto cell faces will be denoted by $[\mu]_{f}$ and a variable defined on cell faces will be denoted by $[\mu]_{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):

[TABLE]

where $w$ is the vertical velocity. $[\nu_{M10}^{q,r}]_{f}$ has various degrees of freedom in the choice of interpolations, such as $[S_{ij}]_{f}\frac{[\eta_{i}^{q}]_{f}}{[\eta_{j}^{r}]_{f}}$ or $\left[S_{ij}\frac{\eta_{i}^{q}}{\eta_{j}^{r}}\right]_{f}$, for example. We will use $\frac{[S_{ij}\ \eta_{i}^{q}]_{f}}{[\eta_{j}^{r}]_{f}}$, once again following [25]. The momentum transfers for method 2 become

[TABLE]

where $[\lambda_{Aij}]_{f}=\frac{\Delta t[S_{ij}]_{f}}{1+\alpha_{M}\Delta t([S_{ji}]_{f}+[S_{ij}]_{f})}$. $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:

[TABLE]

This interpolation method ensures kinetic energy is conserved when converting to the cell centre values [18].

3.6 Summary of proposed transfer terms

We have presented 6 numerical transfer schemes which maintain positivity of mass and conserve mass, momentum, potential energy and internal energy for $\Delta tS_{ij}\leq 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 $\Delta tS_{ij}\leq 1$. Only schemes 2, 4 and 6 are kinetic-energy-diminishing for all timesteps meaning schemes 1, 3 and 5 can cause numerical instabilities if $\Delta tS_{ij}$ is large. From our analysis, we recommend schemes 4 and 6 as they fulfil all the numerical criteria we have set. Scheme 2 is also a viable scheme if the transfer rate is limited to $\Delta tS_{ij}\leq 1$ to maintain positive mass. Section 4 will test these schemes on 2D staggered grids.

4 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 [25] using operator splitting. We will run test cases adapted from the single-fluid rising bubble test case [defined in 5] where an initially stationary temperature anomaly rises and generates resolved circulations (see figure 3). The domain extends to $x\in[-10,10]~{}$km and $z\in[0,10]~{}$km with uniform grid spacings $\Delta x=\Delta z=100~{}$m and wall boundaries on all sides (where zero-gradient fields are imposed and no fluxes perpendicular to the boundaries). A uniform potential temperature field of $\theta=300~{}$K is initially chosen with the system in hydrostatic balance and zero velocity. A warm temperature perturbation is then applied at $t=0~{}$s:

[TABLE]

The perturbation is only applied for $L\leq 1$ where $L\equiv\sqrt{\frac{x-x_{c}}{x_{r}}+\frac{z-z_{c}}{z_{r}}}$, $x_{c}=10~{}$km, $z_{c}=2~{}$km and $x_{r}=z_{r}=2~{}$km. For the 2-fluid experiments the warm anomaly will be applied to $\theta_{1}$ only, whereas fluid 0 will remain initialised as $\theta_{0}=300~{}$K.

We use a 2D C-grid with $\eta_{i}$ and $\theta_{i}$ defined at cell centres and the normal component of $\bm{u}_{i}$ defined at cell faces. The time-stepping is centred Crank-Nicolson with a timestep of $\Delta t=2~{}$s. A van-Leer advection scheme is chosen to maintain positivity of the mass of each fluid. All details of the numerical setup and the numerical solvers used are described in [25], with the exception of a numerical adjustment which must be made for an operator-split Crank-Nicolson multi-fluid scheme (described in appendix 6.3).

4.1 Full bubble test case

The first test case is initialised with all mass in fluid 1: $\sigma_{0}=0$, $\sigma_{1}=1$. The transfer rate is chosen to transfer a large quantity of fluid 1 to fluid 0:

[TABLE]

where $\sigma_{\text{min}}=0.1$. This means that the explicit schemes will transfer $10\%$ of the mass in the first timestep (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:

[TABLE]

where $E_{SF}^{n}$ and $E_{MF}^{n}$ are the total energies at timestep $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 $\Delta E_{RSF}^{n}\sim 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 timestep with the implicit mass schemes (3, 4 and 6) producing the smallest energy changes with $\Delta E_{RSF}^{1}=-1.18\times 10^{-15}$ as they have not transferred the full $10\%$ of the mass in the first timestep (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 unconservative schemes behave similarly to Schemes 1-6, but internal energy and momentum are not conserved exactly in these schemes.

4.2 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 2-fluid test case from [25], where half the mass is initialised with the warm anomaly (fluid 1) and the other half without (fluid 0):

[TABLE]

The 2-fluid equations with different fluid properties require some form of stabilization [25, 22]. We will use a diffusive mass transfer to couple the fluids:

[TABLE]

where $K_{\sigma}=200~{}$m${}^{2}~{}$s*-1* is a large-enough diffusion coefficient to maintain numerical stability for this test case [as shown in 25].

The temperature and volume fraction distributions for this test case are shown in figure 4 - slower circulations form compared with 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 timesteps and produce large energy increases whereas some schemes (which use implicit transfers) follow similar energy evolutions to the conservative schemes. Note that the energy changes due to the numerics of the Crank-Nicolson 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 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.

5 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 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 & 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 6 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 $\Delta tS_{ij}\leq 1$ and $\Delta 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 $\Delta 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 but scheme 2 (when enforcing mass positivity) and scheme 4 are also good candidates. By using any of these three 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 in the multi-fluid equations should be the focus of future studies so that convective processes can be accurately represented.

6 Appendix

6.1 Boundedness properties

For a 2-fluid system, bounded velocity transfer terms can be generalised by

[TABLE]

where $0\leq\beta_{ij}\leq 1$ ensures the new velocities are bounded. Method 1 has $\beta_{ij}=\nu_{Mij}^{q,r}$ (defined in equation 17). $\nu_{Mij}^{q,r}$ is clearly positive given positive mass and transfer rates. To investigate whether $\nu_{Mij}^{q,r}\leq 1$, we make the denominator small (the worst case scenario) such that $S_{ji}=0~{}$s*-1*. This gives

[TABLE]

This boundedness condition is only guaranteed if $\alpha_{M}=1$.

When $\alpha_{M}=\alpha_{C}$, method 2 has $\beta_{ij}=\frac{\lambda_{Cij}\eta_{i}^{m}}{\lambda_{Cij}\eta_{i}^{m}+(1-\lambda_{Cji})\eta_{j}^{m}}$ which is bounded if $0\leq\lambda_{Cij}\leq 1$. This is always true for $\Delta tS_{ij}\leq 1$, although boundedness is also guaranteed for $\Delta tS_{ij}>1$ when $\alpha_{C}=\alpha_{M}=1$.

6.2 Energy properties

For method 2, the total momentum is conserved if the new momenta ($\bm{F}_{i}^{n+1}$) satisfy

[TABLE]

where $\bm{F}_{i}^{m}\equiv\eta_{i}^{m}\bm{u}_{i}^{m}$. The new kinetic energy after transfers have been applied is

[TABLE]

where $\Delta K\equiv\frac{1}{2}\left(\bm{u}_{0}^{m}-\bm{u}_{1}^{m}\right).\left(\mu_{01}\bm{u}_{0}^{m}-\mu_{10}\bm{u}_{1}^{m}\right)$ and $\mu_{ij}\equiv\left[\lambda_{Mij}(1-\beta_{ij}-\beta_{ji})+\beta_{ji}\right]\eta_{i}^{m}$. When $\mu_{01}=\mu_{10}$ and $\mu_{01}\geq 0$, the kinetic energy will never increase. For method 2 (and when $\alpha_{C}=\alpha_{M}$), $\mu_{ij}$ is given by

[TABLE]

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 32 - $\eta_{0}u_{1}$ and $\eta_{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 2.a.

6.3 Numerical adjustments for an operator-split Crank-Nicolson multi-fluid scheme

The numerical multi-fluid scheme used for this study follows the implementation by [25], with the exception of operator-split transfers. For a fluid property such as temperature or velocity ($\phi$) the solution for the Crank-Nicolson scheme (before transfers) is given by

[TABLE]

where $\alpha$ is the off-centering coefficient. As $\left(\frac{\partial\phi_{i}}{\partial t}\right)^{n}$ is stored from the previous timestep, we must ensure that it remains consistent with the fluid properties when transfers are made. This is done by computing

[TABLE]

for method 1 schemes and

[TABLE]

for method 2 schemes. Absence of these terms lead 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-centering coefficient is set to $\alpha=1$ but the scheme will be limited to first-order accuracy in time.

acknowledgements

The authors acknowledge funding from the NERC RevCon project NE/N013735/1 lead by Bob Plant. RevCon is part of the ParaCon project lead by Alison Stirling at the UK Met Office.

\printendnotes