Stochastic Galerkin method for cloud simulation
Abstract.
We develop a stochastic Galerkin method for a coupled NavierStokescloud system that models dynamics of warm clouds. Our goal is to explicitly describe the evolution of uncertainties that arise due to unknown input data, such as model parameters and initial or boundary conditions. The developed stochastic Galerkin method combines the spacetime approximation obtained by a suitable finite volumefinite method with a spectraltype approximation based on the generalized polynomial chaos expansion in the stochastic space. The resulting numerical scheme yields a secondorder accurate approximation in both space and time and exponential convergence in the stochastic space. Our numerical results demonstrate the reliability and robustness of the stochastic Galerkin method for some typical atmospheric scenarios.
1. Introduction
Clouds constitute one of the most important component in the Earthatmosphere system. They influence the hydrological cycle and by interacting with radiation they control the energy budget of the system. However, clouds are one of the most uncertain components, which, unlike the atmospheric flows, cannot be modeled using first principles of physics.
Clouds are composed by myriads of water particles in different phases (liquid and solid), and thus they need to be described by a large ensemble in a statistical sense. A common way of obtaining such an ensemble is by using a mass or size distribution, which would lead to a Boltzmanntype evolution equation. Although there are some approaches available in literature to formulate cloud models in such a way [beheng2010, khvorostyanov1995, khain_etal2000], a complete and consistent description is missing. Since measurements of size distributions of cloud particles are difficult, we are often restricted to averaged quantities such as, for example, mass of water per dry air (mass concentrations). Therefore, models are often formulated in terms of socalled bulk quantities, that is, mass and number concentrations of the respective water species. Many cloud processes are necessary to describe the time evolution of the cloud as a statistical ensemble, that is, particle formation or annihilation, growth/evaporation of particles, collision processes, and sedimentation due to gravity. For each of the processes, we have to formulate a representative mathematical term in the sense of a rate equation. Although for some processes the physical mechanisms are quite understood, the formulation of the process rates usually contain uncertain parameters, thus cloud models come with inherent uncertainty. On the other hand, the initial conditions for atmospheric flows and the embedded clouds are also not perfectly constrained, leading to uncertainties in the environmental conditions. It is wellknown from former studies that uncertainties in cloud processes and in environmental conditions can lead to drastic changes in simulations, thus these uncertainties influence predictability of moist atmospheric flows, clouds and precipitation in a crucial way; for instance, the distribution of latent heat is changed, which in turn can influence frontogenesis [igel_vandenheever2014] or convection [grabowski2015, marinescu_etal2017].
For investigations of the impact of these uncertain cloud model parameters as well as the impact of variations in environmental conditions on atmospheric flows, sensitivity studies are usually carried out. Since one or more parameters are (randomly) varied, the Monte Carlo approach can be used. This, however, requires a large ensemble of simulations to be conducted, which makes Monte Carlo methods computationally expensive and requires a very fine sampling of the parameter space and possible environmental conditions. We therefore choose a different way of representing random variations by using spectral expansions in the stochastic space. This approach enables us to investigate the impact of variations in cloud model parameters and initial conditions on the evolution of moist flows with embedded clouds.
We consider a mathematical model of cloud dynamics that consists of the NavierStokes equations coupled with the cloud evolution equations for the water vapor, cloud water and rain. In this model developed in [cloudpaper, porz18] and presented in Section 2, the NavierStokes equations describe weakly compressible flows with viscous and heat conductivity effects, while microscale cloud physics is modeled by the system of advectiondiffusionreaction equations.
In this paper, we study a stochastic version of the coupled NavierStokescloud model in order to account for uncertainties in input quantities, such as model parameters, initial and boundary conditions, etc. Our main goal is to design an efficient numerical method for quantifying uncertainties in solutions of the studied system. In recent years, a wide variety of uncertainty quantification methods has been proposed and investigated in the context of physical and engineering applications. These methods include stochastic Galerkin methods based on generalized polynomial chaos (gPC) [XK02, XG08, XP12, CGH, PDL, TLE12, TLNE10a, TLNE10b, LSK04, WK06, HSX, JXZ, DPL13], stochastic collocation methods [MZ09, WLB09, XH05, TZ10], and multilevel Monte Carlo methods [SMS13, MSS13, MSS12a, MSS12b, MS12]. Each of these groups of methods has its own pros and cons. While results obtained by the Monte Carlo simulations are generally good, the approach is not very efficient due to a large number of realizations required. Stochastic collocation methods are typically more efficient than the Monte Carlo ones, since they only require solving the underlying deterministic system at the certain quadrature nodes in the stochastic space. These data are then used to reconstruct the gPC expansion using an appropriate set of orthogonal polynomials. Stochastic Galerkin methods offer an alternative approach for computing the gPC expansion. In general, they are more rigorous and efficient than the Monte Carlo and collocation ones; see, e.g., [EMPT11]. On the other hand, application of the stochastic Galerkin method requires changes in the underlying code, because in this approach one needs to solve a system of PDEs for the gPC expansion coefficients.
We develop a new stochastic Galerkin method for the coupled NavierStokescloud system. We restrict our consideration to the case in which the uncertainties are only in the cloud dynamics; extension to the full stochastic NavierStokescloud model is left to future studies. Thus, we need to solve the deterministic NavierStokes equations coupled with the PDE system for the gPC expansion coefficients for the cloud variables. Our numerical method is an extension of the approach proposed in [cloudpaper] for the purely deterministic version of the coupled NavierStokescloud system. This method is based on the operator splitting approach, in which the system is split into the macroscopic NavierStokes equations and microscopic cloud model with random inputs. The NavierStokes equations are then solved by an implicitexplicit (IMEX) finitevolume method, while for the cloud equations we develop a stochastic Galerkin method based on the gPC. The resulting gPC coefficient system is numerically solved by a finitevolume method combined with an explicit RungeKutta method with an enlarged stability region [Dumka3].
The paper is organized as follows. We start in Section 2 with the description of the deterministic NavierStokescloud model. We then continue in Section 3 with the presentation of the stochastic model. Sections 4 and 5 are devoted to the numerical method for the deterministic and stochastic models, respectively. Finally, in Sections 6 and 7, we report on numerical experiments for wellknown meteorological benchmarks—rising warm bubble and RayleighBénard convection—for the deterministic and stochastic models, respectively. Our numerical results clearly demonstrate that the proposed stochastic Galerkin method is capable of quantifying the uncertainties of complex atmospheric flows.
2. Deterministic mathematical model
We study a mathematical model of cloud dynamics, which is based on the compressible nonhydrostatic NavierStokes equations for moist atmosphere (that is, mixture of ideal gases dry air and water vapor),
(2.1)  
and evolution equations for cloud variables,
(2.2)  
Here, is the density, is the velocity vector, is the moist potential temperature, is the pressure, is the acceleration due to gravity, is the dynamic viscosity, the thermal conductivity, and the cloud diffusivity. We denote by the time variable and by the space vector; in the threedimensional (3D) and in the twodimensional (2D) cases. Furthermore, and in the 3D and 2D cases, respectively. We set and . Note that the systems (2) and (2) are coupled through the source term , which represents the impact of phase changes and will be defined below, see (2.6). The temperature can be obtained from the moist adiabatic ideal gas equation
(2.3) 
where is the reference pressure at sea level. In addition to the usual definition of a potential temperature, we use with the ideal gas constant of dry air , the gas constant of water vapor and the specific heat capacity of dry air for constant pressure . In order to close the system, we determine the pressure from the equation of state that includes moisture
(2.4) 
We note that in the dry case reduces to , and the moist ideal gas equation as well as the moist equation of state become their dry analogon.
In this paper, we restrict our investigations to clouds in the lower part of the troposphere, that is, to clouds consisting of liquid droplets exclusively. All of the processes involving ice particles are left for future research. For the representation of liquid clouds in our model we use the socalled single moment scheme, that is, equations for the bulk quantities of mass concentrations of different water phases. For the representation of the relevant cloud processes we adapt a recently developed cloud model [porz18]. Note that for bulk models, the process rates cannot be derived completely from first principles. Consequently, some uncertain parameters show up naturally. This underlies the need of a rigorous sensitivity study which is the goal of the present paper.
Generally, we follow the standard approach in cloud physics modeling for separating hydrometeors of different sizes, as firstly introduced in [kessler1969]. This relies on the observations that small droplets have a negligible falling velocity. In addition, measurements indicate two different modes of droplets in the size distribution, which can be associated to small cloud droplets and large rain drops [warner1969]. Thus, we use the cloud variables and indicating mass concentration of (spatially stationary) cloud droplets and (falling) rain drops, respectively, and the water vapor concentration , that is,
The rest of this section is devoted to a description of the different terms on the righthand side (RHS) of (2), which represents the following relevant cloud processes.

Condensation/evaporation of cloud droplets
Cloud droplets can be formed by the activation of socalled cloud condensation nuclei. Liquid aerosol particles can grow by water vapor uptake to larger sizes; this effect can be described by the Köhler theory; see, e.g., [koehler1936, petters_kreidenweis2007]. As described in detail in [porz18], we represent the cloud droplet number concentration by a nonlinear relationship
Growth and evaporation of small cloud particles are dominated by diffusion processes. If the water vapor concentration is larger than the equilibrium water vapor concentration , water molecules diffuse to the water droplet and thus the cloud particle is growing. For , the water droplet is evaporating. These effects are represented in the terms and . In particular,
(2.5) where
Note that since cloud droplets are activated for water vapor concentrations larger than the thermodynamic equilibrium () the term is added as a source of liquid water in (2.5).
We introduce an additional closure for the number concentration of rain drops , which is explicitly used in [porz18]. Under the assumption that the size of rain drops is distributed to an exponential law; see, e.g., [marshall_palmer1948], we obtain the exponent . Note that this relation will be inserted in any formulation of cloud process rates, involving . Finally, the evaporation of rain drops is changed by hydrodynamic effects of air motions around the drops. This is corrected by an additional empirical relationship. The final formulation of the evaporation rate is given by
where

Autoconversion: Collision of cloud droplets, forming rain drops
The growth of cloud droplets to larger sizes is dominated by collision processes. The collision of two cloud droplets leading to a larger rain drop is called autoconversion; see, e.g., [khain2015]. This rate can be modeled as
Note that the coefficient cannot be measured or derived from the first principles. It is a free parameter, which must be fixed using parameter estimations. Thus, the impact of the uncertainty of this parameter is of high interest. In our deterministic experiments, we choose .

Accretion: Collection of cloud droplets by rain drops
Falling rain drops can also collect smaller cloud droplets. This process is called accretion and can be modelled as
Again, the parameters and cannot be derived from the first principles and the impact of their uncertainty is of high interest as well. In our deterministic experiments, we chose and .

Sedimentation of rain drops
Large rain drops are accelerated by gravity force. Frictional forces balance gravity, thus we can assume that a rain drop falls with a terminal velocity, which depends only on the mass of the drop and the density of air. The terminal velocity is given by
We have to introduce an additional hyperbolic term into the equation for the evolution of , that is, the term is included.
Note that the condensation and evaporation processes are formulated explicitly, in contrast to the usual approach of saturation adjustment (see, e.g., [LV11]), which is less accurate, but commonly used in operational weather forecast models. This explicit formulation introduces stiffness caused by modeling cloud processes on the RHS of the cloud equations with fractional exponents between and . We handle this stiffness by replacing terms like , , with
Due to the condensation and evaporation processes latent heat is released or absorbed. These processes are modelled by the source term in (2):
(2.6) 
where is the specific latent heat of vaporization.
Solving the NavierStokes equations (2) in a weakly compressible regime is known to cause numerical instabilities due to the multiscale effects. We follow the approach typically used in meteorological models, where the dynamics of interest is described by a perturbation of a background state, which is the hydrostatic equilibrium. The latter expresses a balance between the gravity and pressure forces. Denoting by , , , and the respective background state, the hydrostatic equilibrium satisfies
where is obtained from the equation of state (2.4)
(2.7) 
Let , , , and stand for the corresponding perturbations of the equilibrium state, then
The pressure perturbation is derived from (2.4) and (2.7) using the following Taylor expansion
which results in
The perturbation formulation of the NavierStokes equations (2) then reads as
(2.8)  
Note that though the systems (2) and (2) are equivalent, the perturbation formulation (2) is preferable for the development of a numerical scheme. For alternative representations of cloud dynamics and their numerical investigations, we refer the reader to [kroner, COSMO] and references therein.
3. Stochastic mathematical model
In the meteorological model described in Section 2, some data or parameters may contain uncertainty. In this paper, we consider the case, where the uncertainty arises from the initial data or some coefficients in the microphysical cloud parameterizations. In order to mathematically describe the uncertainty, we introduce a random variable . We assume that either the initial data or some wellchosen model parameters depend on , that is,
or
Consequently, the solution at later time will also depend on , that is, for , and the system (2) for cloud variables will read as
(3.1)  
From now on we will stress the dependence on , but we will omit the dependence on and to simplify the notation. We would like to point out that the solution of the NavierStokes equations (2) will also depend on , because of the source term . In this paper, we will consider a simplified situation by replacing
in (2) by which only depends on the expected values of the cloud variables
This ensures that all of the fluid variables, , and , remain deterministic.
4. Numerical scheme for the deterministic model
The numerical approximation of the coupled model (2), (2) is based on the secondorder Strang operator splitting. Therefore, we split the whole system into the macroscopic NavierStokes flow equations and the microscopic cloud equations. The NavierStokes equations (2) are approximated by an IMEX finitevolume method and the cloud equations (2) are approximated by a finitevolume method in space and an explicit RungeKutta method with an enlarged stability region in time.
4.1. Operator form
Let and denote the solution vectors of (2) and (2), respectively. Then, the coupled system can be written as
where and are advection fluxes and , and , denote the diffusion and reaction operators of the respective systems. They are given by
(4.1)  
In order to derive an asymptotically stable, accurate and computational efficient scheme for the NavierStokes equations, we first split the equations into linear and nonlinear parts; see [bispen1, cloudpaper] and references therein. Consequently, we introduce

with and ;

with

with and .
We would like to point out that the choice of the linear and nonlinear operators is crucial. We choose the linear part to model linear acoustic and gravitational waves as well as linear viscous fluxes. The nonlinear part describes nonlinear advective effects together with the remaining nonlinear viscous fluxes and the influence of the latent heat. We will use the following notation:
4.2. Discretization in space
The spatial discretization is realized by a finitevolume method. We take a cuboid computational domain , which is divided into uniform Cartesian cells. The cells are labelled in a certain order using a singleindex notation. For simplicity of notation, we assume that the cells are cubes with the sides of size so that . We also introduce the notation for the set of all neighbouring cells of cell , .
We assume that at a certain time the approximate solution is realized in terms of its cell averages
In order to simplify the notation, we will now omit the time dependence of and . Next, we introduce the notation and and consider the following approximation of the advection, diffusion and reaction operators:
Analogous notation will be used for the approximations , and of the cloud operators.
4.2.1. Advection
The advection terms are discretized using flux functions as follows:
where the numerical fluxes , and approximate the corresponding fluxes between the computational cells and , and denotes the th component of the outer normal unit vector of cell in the direction of cell . We use the Rusanov numerical flux for and and the central flux for . For and a secondorder discretization is obtained via a MUSCLtype approach using piecewise linear reconstructions with the minmod limiter. The numerical fluxes are then given by
(4.2)  
Here, , and , denote the corresponding interface values, which are computed using a piecewise linear reconstruction so that
where the slopes are computed by the mindmod limiter,
applied in a componentwise manner. Here,
and and are obtained similarly. Thereby is the other neighboring cell of in the opposite direction from . Finally, the values and are given by
where denotes the spectral radius of the corresponding Jacobians.
Remark 4.1.
Note that in the computation of in (4.2.1), we use the cell averages rather than the point values at the cell interfaces for the following two reasons. First, the flux is secondorder accurate. Second, in Section 4.3, we will treat the linear part of the flux implicitly and this is much easier to do when the numerical flux is linear as well.
4.2.2. Diffusion
The components of the discrete diffusion operators are discretized in a straightforward manner using secondorder central differences.
4.2.3. Reaction
The reaction terms are discretized by a direct evaluation of the reaction operators at the cell centres:
After the spatial discretization, we obtain the following system of timedependent ODEs:
(4.3)  
(4.4) 
This system has to be solved using an appropriate ODE solver as discussed in Section 4.3.
4.3. Discretization in time
Let and denote the numerical approximation of the solutions and at the discrete time level . We evolve the solution to the next time level , where is the size of the Strang operator splitting time step. In the operator splitting approach, we first numerically solve the ODE system (4.3) with , we then numerically integrate the ODE system (4.4) with and finally we solve the system (4.3) again with .
Notice that the system (4.3) may be very stiff as the NavierStokes equations are in the weakly compressible regime. We therefore follow the approach in [bispen1] (see also [diss]), and employ the secondorder ARS(2,2,2) IMEX method from [ARS222]:
(4.5)  
where , , , , and satisfies the following CFL condition:
For solving the linear systems arising in (4.5), we use the generalized minimal residual (GMRES) method combined with a preconditioner, the incomplete LU factorization (ILU). As it was shown in [bispen1] (see also [diss]), the resulting method is both accurate and efficient in the weakly compressible regime.
The ODE system (4.4) is also stiff, but its stiffness only comes from the diffusion and powerlawtype source terms. We therefore efficiently solve it using the large stability domain thirdorder RungeKutta method from [Dumka3]. We have utilized the ODE solver DUMKA3, which is a free software that can be found in [Dumka3_code]. We note that DUMKA3 selects time steps automatically, but in order to improve its efficiency, one needs to provide the code with a time step stability restriction for the forward Euler method; see [Dumka3_code, Dumka3]. This bound is obtained by , where satisfies the following CFL condition for the cloud system:
5. Numerical scheme for the stochastic model
In this section, we describe a generalized polynomial chaos stochastic Galerkin (gPCSG) method for the system of cloud equations (3). Such method belongs to the class of intrusive methods and the use of the Galerkin expansion leads to a system of deterministic equations for the expansion coefficients. In the gPCSG method, the solution is sought in the form of a polynomial expansion
(5.1) 
where , , are polynomials of th degree that are orthogonal with respect to the probability density function . More precisely, the polynomials satisfy
(5.2) 
where is the Kronecker symbol and is the sample space. The choice of the orthogonal polynomials depends on the distribution of . In our case, we use a uniformly distributed , which defines the Legendre polynomials. We use the same expansion for the uncertain coefficients,
(5.3) 
for the source terms on the RHS of (3),
(5.4)  
as well as for the raindrop fall velocity,
(5.5) 
Since , we also obtain
(5.6) 
We note that if is very small, the computation of the coefficients should be desingularized; see, e.g.,[KurganovPetrova], where several desingularization strategies were discussed.