The representation of subgrid‐scale convection is a weak aspect of weather and climate prediction models and the assumption that no net mass is transported by convection in parametrizations is increasingly unrealistic as models enter the grey zone, partially resolving convection. The solution of conditionally averaged equations of motion (multifluid equations) is proposed in order to avoid this assumption. Separate continuity, temperature, and momentum equations are solved for inside and outside convective plumes, which interact via mass‐transfer terms, drag, and by a common pressure. This is not a convection scheme that can be used with an existing dynamical core—this requires a whole new model. This article presents stable numerical methods for solving the multifluid equations, including large transfer terms between the environment and plume fluids. Without transfer terms, the two fluids are not sufficiently coupled and solutions diverge. Two transfer terms are presented, which couple the fluids together in order to stabilize the model: diffusion of mass between the fluids (similar to turbulent entrainment) and drag between the fluids. Transfer terms are also proposed to move buoyant air into the plume fluid and vice versa as would be needed to represent initialization and termination of subgrid‐scale convection. The transfer terms are limited (clipped in size) and solved implicitly in order to achieve bounded, stable solutions. Results are presented for a well‐resolved warm bubble with rising air being transferred to the plume fluid. For stability, equations are formulated in advective rather than flux form and solved using bounded finite‐volume methods. Discretization choices are made to preserve boundedness and conservation of momentum and energy when mass is transferred between fluids. The formulation of transfer terms in order to represent subgrid convection is the subject of future work.

Traditional parameterizations of convection are a large source of error in weather and climate prediction models, and the assumptions behind them become worse as resolution increases. Multifluid modeling is a promising new method of representing subgrid‐scale and near‐grid‐scale convection allowing for net mass transport by convection and nonequilibrium dynamics. The air is partitioned into two or more fluids, which may represent, for example, updrafts and the nonupdraft environment. Each fluid has its own velocity, temperature, and constituents with separate equations of motion. This paper presents two‐fluid Boussinesq equations for representing subgrid‐scale dry convection with sinking and w = 0 air in Fluid 0 and rising air in Fluid 1. Two vertical slice test cases are developed to tune parameters and to evaluate the two‐fluid equations: a buoyant rising bubble and radiative convective equilibrium. These are first simulated at high resolution with a single‐fluid model and conditionally averaged based on the sign of the vertical velocity. The test cases are next simulated with the two‐fluid model in one column. A model for entrainment and detrainment based on divergence leads to excellent representation of the convective area fraction. Previous multifluid modeling of convection has used the same pressure for both fluids. This is shown to be a bad approximation, and a model for the pressure difference between the fluids based on divergence is presented.

The multifluid equations are derived from the compressible Euler equations (or any of the usual approximate equation sets used in meteorology) by conditional filtering. They have the potential to provide the basis for an improved representation of cumulus convection and its coupling to the boundary layer and larger scale flow in numerical models. The present article derives the prognostic equations for subfilter-scale turbulent second moments in the multifluid framework, along with certain systematic simplifications of them, thus providing a multifluid analogue of the well-known Mellor and Yamada hierarchy of turbulence closures. As well as enabling a more accurate calculation of subfilter-scale fluxes and the effects of subfilter-scale variability on cloud fraction, liquid water, and buoyancy, the second moment information can be used to obtain a more accurate parameterization of entrainment and detrainment. A subset of the turbulence equations derived here is employed in the two-fluid single-column model described in Part 2 and applied to the simulation of shallow cumulus cases in Part 3.

Two-fluid modelling has recently emerged as a promising approach to representing cumulus convection in weather and climate models. This study applies the two-fluid model described in Part II (Thuburn et al., 2022b) to a shallow cumulus convection case study over land (ARM). Large eddy simulation data is used to tune the majority of the closures that determine the properties of entrained and detrained air. The two-fluid model is generally able to reproduce the profiles of the mean and turbulent quantities over all stages of the diurnal cycle. As such, the initiation of shallow convection and the evolution of the cloud layer are well captured. The robustness of the two-fluid model is further verified using a steady-state test case (BOMEX), in which the cloud properties are also well modelled.

Multi‐fluid models have recently been proposed as an approach to improving the representation of convection in weather and climate models. This is an attractive framework as it is fundamentally dynamical, removing some of the assumptions of mass‐flux convection schemes which are invalid at current model resolutions. However, it is still not understood how best to close the multi‐fluid equations for atmospheric convection. In this paper we develop a simple two‐fluid, single‐column model with one rising and one falling fluid. No further modelling of sub‐filter variability is included. We then apply this model to Rayleigh–Bénard convection, showing that, with minimal closures, the correct scaling of the heat flux (Nu) is predicted over six orders of magnitude of buoyancy forcing (Ra). This suggests that even a very simple two‐fluid model can accurately capture the dominant coherent overturning structures of convection.