Stochastic Galerkin method for cloud simulation

We develop a stochastic Galerkin method for a coupled Navier-Stokes-cloud 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 space-time approximation obtained by a suitable finite volume method with a spectral-type approximation based on the generalized polynomial chaos expansion in the stochastic space. The resulting numerical scheme yields a second-order 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. We also use the proposed method to study the behavior of clouds in certain perturbed scenarios, for examples, the ones leading to changes in macroscopic cloud pattern as a shift from hexagonal to rectangular structures.


Introduction
Clouds constitute one of the most important component in the Earth-atmosphere 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 Boltzmann-type evolution equation. Although there are some approaches available in literature to formulate cloud models in such a way [4,17,18], 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 so-called 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 well-known 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 [15] or convection [14,24].
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. In most practical studies, a much smaller set of ensembles (with about 10 − 100 samples only) is used.
In order to improve both the efficiency and accuracy of a numerical method, we 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 physics that consists of the Navier-Stokes equations coupled with the cloud evolution equations for the water vapor, cloud water and rain. In this model developed in [22,34] and presented in Section 2, the Navier-Stokes equations describe weakly compressible flows with viscous and heat conductivity effects, while microscale cloud physics is modeled by the system of advection-diffusionreaction equations.
Meteorological applications typically inherit several sources of uncertainties, such as model parameters, initial and boundary conditions. Consequently, purely deterministic models are insufficient in such situations and more sophisticated methods need to be applied to analyze the influence of uncertainties on numerical solutions. In this paper, we study a stochastic version of the coupled Navier-Stokes-cloud model in order to account for uncertainties in input quantities. 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) [10,33,40,41,46], stochastic collocation methods [23,44,45], and multilevel Monte Carlo methods [28,29,37]. 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. The Monte Carlo as well as the stochastic collocation method fall into a class of the non-intrusive methods. 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., [12]. The stochastic Galerkin method is an intrusive method since it requires changes in the underlying code. In fact, one needs to solve a system of PDEs for the gPC expansion coefficients.
We develop a new stochastic Galerkin method for the coupled Navier-Stokes-cloud system. As it has already been mentioned above, the largest source of uncertainties is cloud physics. Therefore, we restrict our consideration to the case in which the uncertainties are only in the cloud physics representation; extension to full stochastic Navier-Stokes-cloud model is left to future studies. Thus, we need to solve the deterministic Navier-Stokes 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 [22] for the purely deterministic version of the coupled Navier-Stokes-cloud system. This method is based on the operator splitting approach, in which the system is split into the macroscopic Navier-Stokes equations and microscopic cloud model with random inputs. The Navier-Stokes equations are then solved by an implicit-explicit (IMEX) finite-volume 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 finite-volume method combined with an explicit Runge-Kutta method with an enlarged stability region [27].
The paper is organized as follows. We start in Section 2 with the description of the deterministic Navier-Stokes-cloud model. The numerical method for the deterministic model is presented in Sections 3. In Sections 4, we report on numerical experiments for well-known meteorological benchmarks-rising warm bubble and Rayleigh-Bénard convection-for the deterministic model. We then continue in Section 5 with the presentation of the stochastic model, which is followed by the description of the numerical method (Section 6) and presentation of the numerical experiments (Section 7) for the stochastic model. Our numerical results clearly demonstrate that the proposed stochastic Galerkin method is capable of quantifying the uncertainties of complex atmospheric flows.

Deterministic mathematical model
We study a mathematical model of cloud dynamics, which is based on the compressible nonhydrostatic Navier-Stokes equations for moist atmosphere (that is, mixture of ideal gases dry air and water vapor), and evolution equations for cloud variables, Here, ρ is the density, u = (u 1 , u 2 , u 3 ) is the velocity vector, θ is the moist potential temperature, p is the pressure, g is the acceleration due to gravity, µ m is the dynamic viscosity, µ h is the thermal conductivity, and µ q is the cloud diffusivity. The cloud variables representing the mass concentration of water vapor, cloud droplets and rain drops, q v , q c , q r , respectively, as well as the right-hand side (RHS) terms E, C, A 1 , A 2 will be defined below. We denote by t the time variable and by x the space vector; x = (x 1 , x 2 , x 3 ) in the three-dimensional (3-D) and x = (x 1 , x 3 ) in the two-dimensional (2-D) cases. Furthermore, e 3 = (0, 0, 1) and e 3 = (0, 1) in the 3-D and 2-D cases, respectively. We set µ m = 10 −3 and µ h = 10 −2 = µ q . Note that the systems (2.1) and (2.2) are coupled through the source term S θ , which represents the impact of phase changes and will be defined below, see (2.6). The temperature T can be obtained from the moist adiabatic ideal gas equation where p 0 = 10 5 Pa is the reference pressure at sea level. In addition to the usual definition of a potential temperature, we use R m = (1−q v −q c −q r )R+q v R v with the ideal gas constant of dry air R = 287.05 J /(kg·K), the gas constant of water vapor R v = 461.51 J /(kg·K) and the specific heat capacity of dry air for constant pressure c p = 1005 J /(kg·K). In order to close the system, we determine the pressure from the equation of state that includes moisture We note that in the dry case R m reduces to R, S θ = 0 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 so-called 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 [34]. 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 [16]. 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 [42]. Thus, we use the cloud variables q c and q r indicating mass concentration of (spatially stationary) cloud droplets and (falling) rain drops, respectively, and the water vapor concentration q v , that is, q = mass of the respective phase mass of dry air for ∈ {v, c, r}.
The rest of this section is devoted to a description of the different terms on the RHS of (2.2), which represents the following relevant cloud processes, see [34].
• General properties of a single water particle As we exclusively investigate water clouds, we can assume a spherical shape of water particles. For small cloud droplets this is a very good approximation, while for large rain drops drag effects change their shape [38,39]. However, for our investigations of ensembles of rain drops, the spherical shape approximation is appropriate. Thus, mass and radius of droplets are related by the usual equation where the Schmidt and Reynolds numbers are defined as (2.5) respectively. In (2.5), µ is the dynamic viscosity of air, which is expressed according to [11] by For cloud droplets, we neglect the ventilation correction, thus the mass rate of diffusion for a cloud droplet of mass m c can be expressed as For rain drops, growth due to the diffusion is negligible, and thus we obtain the mass rate for rain drops of mass m r as Here, (·) + := max(·, 0) denotes the positive part.
• Collision of rain drops with cloud droplets: Accretion A spherical rain drop of mass m r (radius r) falls with terminal velocity v τ (m r ) through a volume V = πr 2 v τ (m r )∆t (during a time interval ∆t) and collects cloud droplets of total mass M c = V ρq c . Thus, the corresponding growth rate of the rain drop is given by r with an efficiency k 2 > 0.

2.2.
Ensemble/collective properties. For the description of clouds as an ensemble of water particles, we would have to introduce such averaged quantities as mass concentrations (as described above, that is, q c and q r ) as well as number concentrations of cloud droplets, n c , and rain drops, n r . Since we do not extend the systems of equations for these two quantities, we introduce relations between mass and number concentrations in order to keep the main effects in a simplified way.
• Formation of cloud droplets: Activation Cloud droplets can be formed by the activation of so-called cloud condensation nuclei (CCN). 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., [19,32]. As described in detail in [34], we represent the cloud droplet number concentration n c by a nonlinear relation Here, N ∞ denotes the maximum number of CCN (depending on environmental conditions, e.g. clean or polluted air), m 0 can be interpreted as the activation mass of cloud droplets, and N 0 is the approximated number of activated droplets at q v = q * . In our investigations, we set these three parameters to the following values: For the initialization of the cloud droplet production, we introduce an additional factor in case of supersaturation • Relation between number and mass concentration for rain drops In contrast to the formulation in [34], we do not include another equation for the number concentration of rain drops. In a similar way as for cloud activation, we use a relation between n r and q r , that is, a closure of the form n r = f (q r , c). Since we implicitly assume that the rain drops are distributed according to their size, this approach should be used for mimicking the shape of the distribution in a proper way. We propose the (non)linear relation n r = c r q γ r , 0 < γ ≤ 1. Assuming a constant mean mass of rain drops m r , we can determine the constants as c r = m −1 r and γ = 1. This approach would be meaningful for the case of a symmetric size distribution of rain droplets, centered around the mean mass. However, it is well-known that size distributions of rain are usually skew to larger sizes, thus a linear relation is not appropriate. For sizes of rain drops often an exponential distribution is assumed, this leads to an exponent γ = • Rates for diffusion processes For the mass mixing ratios q c and q r , the rates for the diffusion processes are given by multiplication of the single particle rates by the number concentration of the respective particles, namely: Applying the relations between the mass and number concentrations as stated above, we obtain Note, that the rates for activation and diffusion growth of cloud droplets are combined in the model formulation, that is, C = C 1 + C act .

• Rate for accretion
For the rate of accretion of rain water, we obtain Note that for compensating effects of the averaging the parameter k 2 can be adjusted (such that k 2 = 0.8 kg = k 2 ) and the impact of the uncertainty of this parameter is of high interest, since it cannot be measured or derived from the first principles.
• Collision of cloud droplets, forming a rain drop: Autoconversion Beside the growth of an existing rain drop by collecting cloud droplets, a rain drop can be formed by the collision of two cloud droplets. According to [34], we formulate the growth rate similarly to population models, namely: Note that the coefficient k 1 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 k 1 = 0.0041 kg s −1 , as indicated in [34].
• Sedimentation of rain mass mixing ratio We have introduced an additional convection term into the equation for the evolution of q r in (2.2), that is, The parameter α can be derived empirically, but the influence of uncertainty in α is of high interest.
Note that activation and diffusion processes are formulated explicitly, in contrast to the usual approach of saturation adjustment (see, e.g., [21]), 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 −1 and 1. In order to avoid the root evaluation for an argument that is close to zero, we introduce a cut-off function and replace ζ ξ , ξ ∈ (−1, 1), with Due to phase changes (activation and growth/evaporation of water particles) latent heat is released or absorbed. These processes are modeled by the source term in (2.1): Solving the Navier-Stokes equations (2.1) 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 byp,ρ, u = 0,θ and ρθ the respective background state, the hydrostatic equilibrium satisfies wherep is obtained from the equation of state (2.4) Let p , ρ , u , θ and (ρθ) stand for the corresponding perturbations of the equilibrium state, then The pressure perturbation p is derived from (2.4) and (2.7) using the following Taylor expansion which results in Consequently, a physical motivation from the hydrostatic balance state leads to a numerically preferable scaling and the perturbation formulation of the Navier-Stokes equations (2.1) then reads as For alternative representations of cloud dynamics and their numerical investigations, we refer the reader to [3,36] and references therein.

Numerical scheme for the deterministic model
The numerical approximation of the coupled model (2.8), (2.2) is based on the second-order Strang operator splitting. Therefore, we split the whole system into the macroscopic Navier-Stokes flow equations and the microscopic cloud equations. The Navier-Stokes equations (2.8) are approximated by an IMEX finite-volume method and the cloud equations (2.2) are approximated by a finite-volume method in space and an explicit Runge-Kutta method with an enlarged stability region in time.
3.1. Operator form. Let w := (ρ , ρu, (ρθ) ) and w q := (ρq v , ρq c , ρq r ) denote the solution vectors of (2.8) and (2.2), respectively. Then, the coupled system can be written as where F and F q are advection fluxes and D, R and D q , R q denote the diffusion and reaction operators of the respective systems. They are given by In order to derive an asymptotically stable, accurate and computational efficient scheme for the Navier-Stokes equations, we first split the equations into linear and nonlinear parts; see [6,22] and references therein. Consequently, we introduce 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:

Discretization in space.
The spatial discretization is realized by a finite-volume method. We take a cuboid computational domain Ω ⊂ R d , which is divided into N uniform Cartesian cells. The cells are labeled in a certain order using a single-index notation. For simplicity of notation, we assume that the cells are cubes with the sides of size h so that |C i | = h d . We also introduce the notation S(i) for the set of all neighboring cells of cell C i , i = 1, . . . , N .
We assume that at a certain time t 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 w i (t) and (w q ) i (t). Next, we introduce the notation and consider the following approximation of the advection, diffusion and reaction operators: Analogous notation will be used for the approximations ( 3.2.1. Advection. The advection terms are discretized using flux functions as follows: where the numerical fluxes H L ij , H N ij and (H q ) ij approximate the corresponding fluxes between the computational cells C i and C j , and n a discretization is obtained via a MUSCL-type approach using piecewise linear reconstructions with the minmod limiter. It is well-known this approach yields an approximation, which is almost second-order accurate as its accuracy deteriorates at extrema. The numerical fluxes are given by Here, w − ij , w + ij and (w q ) − ij , (w q ) + ij denote the corresponding interface values, which are computed using a piecewise linear reconstruction so that where the slopes s ij are computed by the minmod limiter, applied in a component-wise manner. Here, and (w q ) − ij and (w q ) + ij are obtained similarly. Thereby C j * is the other neighboring cell of C i in the opposite direction from C j . Finally, the values λ ij and (λ q ) ij are given by where σ denotes the spectral radius of the corresponding Jacobians.
Remark 3.1. Note that in the computation of H L ij in (3.2), we use the cell averages rather than the point values at the cell interfaces for the following two reasons. First, the flux is second-order accurate. Second, in Section 3.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.

3.2.2.
Diffusion. The components of the discrete diffusion operators are discretized in a straightforward manner using second-order central differences.

3.2.3.
Reaction. The reaction terms are discretized by a direct evaluation of the reaction operators at the cell centers: After the spatial discretization, we obtain the following system of time-dependent ODEs: This system has to be solved using an appropriate ODE solver as discussed in Section 3.3.

Discretization in time.
Let w n h and (w q ) n h denote the numerical approximation of the solutions w h (t) and (w q ) h (t) at the discrete time level t n . We evolve the solution to the next time level t n+1 = t n + ∆t n , where ∆t n is the size of the Strang operator splitting time step. In the operator splitting approach, we first numerically solve the ODE system (3.3) with ∆t n NS = ∆t n /2, we then numerically integrate the ODE system (3.4) with ∆t n and finally we solve the system (3.3) again with ∆t n NS . Notice that the system (3.3) may be very stiff as the Navier-Stokes equations are in the weakly compressible regime. We therefore follow the approach in [6] (see also [5]), and employ the second-order ARS(2,2,2) IMEX method from [2]: where α = 1 − 1 / √ 2, δ = 1 − 1 /2β, t n+ 1 2 = t n + ∆t n NS , t n+ 1 4 = t n + ∆t n NS/2, and ∆t n NS satisfies the following CFL condition: For solving the linear systems arising in (3.5), we use the generalized minimal residual (GMRES) method combined with a preconditioner, the incomplete LU factorization (ILU). As it was shown in [6] (see also [5]), the resulting method is both accurate and efficient in the weakly compressible regime.
The ODE system (3.4) is also stiff, but its stiffness only comes from the diffusion and power-law-type source terms. We therefore efficiently solve it using the large stability domain third-order Runge-Kutta method from [27]. We have utilized the ODE solver DUMKA3, which is a free software that can be found in [26]. 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 [26,27]. This bound is obtained by min{∆t n , ∆t n cloud }, where ∆t n cloud satisfies the following CFL condition for the cloud system:

Deterministic numerical experiments
In this section, we test the numerical method described in Section 3. The experimental order of convergence is computed for the so-called free convection of a moist warm air bubble and the structure formation in cloud dynamics is shown in the Rayleigh-Bénard convection. The latter will be simulated in both the 2-D and 3-D cases.
4.1. Free convection of a moist warm air bubble in 2-D. We start with the well-known meteorological benchmark describing the free convection of a smooth warm air bubble; see, e.g., [7,9].
Example 1. In this experiment, the warm bubble rises and deforms axisymmetrically due to the shear friction with the surrounding air at the warm/cold air interface, gradually forming a mushroom-like shape. The warm air bubble is placed at (3500 m, 2000 m) with the initial perturbation: For the cloud variables we choose the following initial conditions: Furthermore, we apply the no-flux boundary conditions u · n = 0, ∇ρ · n = 0, ∇(ρθ) · n = 0, ∇(ρq ) · n = 0, ∈ {v, c, r}.
In Figure 1, we show the potential temperature θ and cloud variables q v , q c and q r , computed on a 160 × 160 mesh at t = 150 and 200s. One can clearly observe condensation taking place on the interface between cold and warm air and leading to cloud formation in this region. In consequence, rain is formed in the clouds and falls towards the surface. Note that the order of magnitude of the different water mass concentrations is very different, that is, q v q c q r , as expected. The experimental convergence study for the cloud and flow variables is presented in Tables 1 and 2, respectively. The experimental order of convergence (EOC) has been computed in the following way: where v N,∆t is the numerical solution computed on a grid with N × N grid cells and using a fixed time step ∆t. As one can clearly see, the expected second order of accuracy has been achieved. For comparison, we present in Figures 2 and 3 the errors measured in the L 1 -, L 2 -and L ∞ -norms. They all give similar results.

4.2.
Rayleigh-Bénard convection. In this experiment, we study a natural convection that is used to model structure formation. It occurs in a planar flow between two horizontal plates, where the lower one is heated from below and the upper one is cooled from above. Due to the presence of buoyancy, and hence gravity, the fluid develops a regular pattern of convection roles, known as the Bénard cells. In 3-D, these convection roles form additionally hexagonal structures; see, e.g., [1,13,31].
Example 2: 2-D case. In Figures 4-7, we present time snapshots of the numerical solution computed in a domain Ω = [0, 5000] × [0, 1000] m 2 that has been discretized using 320 × 320 mesh cells. The potential temperature, water vapor mixing ratio, cloud mass and rain mass concentration are presented at two distinct times (t = 800 and 1400s) in Figures 4-7, where we can clearly see the formation and evolution of a pattern. At an earlier time t = 800 s, one can observe the formation of small convective cells, visible as narrow finger-like structures reaching towards the top of the domain. Inside these cells, the potential temperature is enhanced, partly due to the upward transport of higher values from below and partly due to phase changes and thus latent heat release. Also the mass concentrations of water vapor and cloud water follow the small scale structure and show enhanced values inside the fingers. Even at this early stage, rain can be formed at the top of the domain, since there the cloud water concentration is high enough for autoconversion and accretion. Nevertheless, the structure of rain water is very different since after the formation of rain, it is vertically transported due to sedimentation leading to vertically smeared structures. By a later time t = 1400 s, much larger structures, which are similar to classical structures for dry thermal convection, have been formed. In the variables θ, q v and q c , the spatially extended convective cells can be clearly seen. In contrast, rain water is not following the convective structure although some larger features can be seen. In general, smearing due to sedimentation is again a major feature of the rain mass concentration.     ( Figure 8) and cloud water ( Figure 10). As time progresses, these structures aggregate and reorganize to quasi-hexagonal structures, which would be typical for classical dry thermal convection-Rayleigh-Bénard convection. These structures in potential temperature seem to be quite robust as the overall scales and pattern do not change from t = 1800 s to t = 2000 s. As one can see in Figure 10, the structures in the cloud water are very similar since the variables θ and q c are closely connected; the structures are mainly visible in horizontal planes. For the rain distribution, the evolved structures are quite different since rain is formed at regions with high cloud water (that is, at the top layers) and then transported by sedimentation leading to more pronounced pattern in the vertical direction, since sedimentation is the dominant process after the rain has formed.

Stochastic mathematical model
In meteorological applications, it is typical that initial or boundary data as well as parameters are uncertain. In order to take this into account and analyze the influence of such uncertainties on the output of numerical simulations, we apply the stochastic Galerkin method that will be described below. In this paper, we consider the case where the uncertainty arises from the initial data or some coefficients in the microphysical cloud parametrizations. In order to mathematically describe the uncertainty, we introduce a random variable ω. We assume that either the initial data or some well-chosen model parameters depend on ω, that is, Consequently, the solution at later time will also depend on ω, that is, (ρq )(x, t, ω) for ∈ {v, c, r}, and the system (2.2) for cloud variables will read as .
From now on we will stress the dependence on ω, but we will omit the dependence on x and t to simplify the notation. We would like to point out that the solution of the Navier-Stokes equations (2.8) will also depend on ω, because of the source term S θ . In this paper, we will consider a simplified situation by replacing in (2.8) byS θ which only depends on the expected values of the cloud variables This ensures that all of the fluid variables, ρ , ρu and (ρθ) , remain deterministic.

Numerical scheme for the stochastic model
In this section, we describe a generalized polynomial chaos stochastic Galerkin (gPC-SG) method for the system of cloud equations (5.1). 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 gPC-SG method, the solution is sought in the form of a polynomial expansion where Φ k (ω), k = 0, . . . , M , are polynomials of k-th degree that are orthogonal with respect to the probability density function µ(ω). The choice of the orthogonal polynomials {Φ k (ω)} M k=0 depends on the distribution of ω. In our case, we use a uniformly distributed ω ∈ Γ = (−1, 1), which defines the Legendre polynomials, which satisfy where δ kk is the Kronecker symbol and Γ is the sample space.
We use the same expansion for the uncertain coefficients, for the source terms on the RHS of (5.1), as well as for the raindrop fall velocity, Since ρ(x, t) = ρ 0 (x, t), we also obtain Applying the Galerkin projection to (5.1) yields for k = 0, . . . , M , where ·, · is the scalar product in our probability space which is given through We now substitute (6.1), (6.4)-(6.6) into (6.7) and use the orthogonality property (6.2) to obtain the following 3(M + 1) × 3(M + 1) deterministic system for the gPC coefficients: Example 4: 2-D case with stochastic initial data. We begin by considering the following experiment with a 10% perturbation of the initial water vapor concentration: We compute the solution using different meshes until the final time t = 10s.
The experimental convergence study for the cloud and flow variables is presented in Figure 12. Similarly to the deterministic case, one can observe second-order convergence in space and time. In order to test the convergence in the stochastic space, we obtain a reference solution computed by the stochastic Galerkin method with 20 stochastic modes. The convergence study is presented in Figure 13, where we plot the difference between the approximate and reference solutions, both computed using a mesh with 160 × 160 cells and ∆t = 0.01 at time t = 10s. One can clearly see a spectral convergence with the rate e −0.3M .  Example 5: 2-D case with stochastic parameters. In this experiment, we perturb the three selected model parameters (k 1 , k 2 and α) by 50% each: These parameters were proposed in [34] as the most sensitive model parameters. We study the convergence in the stochastic space. To this end, we plot in Figure 14 the difference between the approximate and reference (obtained with 20 stochastic modes) solutions, both computed using a mesh with 160 × 160 cells and ∆t = 0.01 at time t = 10s. As in Example 4, one can observe a spectral convergence with the rate e −0.3M . Example 6: Comparison of the stochastic Galerkin and stochastic collocation methods in the 3-D case. In this comparison test, we consider free convection of a moist smooth warm air bubble with perturbed initial data as in Example 4. Numerical solutions are computed on 50 × 50 × 50 mesh at time t = 1s. We compare the performance of the stochastic Galerkin and stochastic collocation methods. The collocation method (see, e.g., [8]) is an interpolation method in the stochastic space, which uses a deterministic model with the values of the stochastic variable taken at collocation points suitably chosen on the interval (−1, 1); here, we use the Gauss-Legendre points.
In Figure 15, the L 2 -norms of the stochastic coefficients of the stochastic Galerkin ( ( ρq ) m L 2 (Ω) ) and stochastic collocation ( ( ρq ) m L 2 (Ω) ) methods for ∈ {v, c, r} and m = 0, . . . , 19 are shown. One can observe an exponential decay with respect to m and good agreement between both methods that demonstrates the reliability of the stochastic Galerkin method. We note, however, that the norms of the solutions computed by these two different methods are not equal since the stochastic Galerkin method uses the expected values of the cloud variables in the Navier-Stokes equations. In Figure 16, we compare the CPU times consumed by the stochastic Galerkin and stochastic collocation methods with the same number M of stochastic modes/collocation points. Since the stochastic Galerkin method solves the Navier-Stokes equations just once instead of M times, as needed by the stochastic collocation method does, it is expected to outperform the stochastic collocation method. This has been confirmed by our simulations. 7.2. Rayleigh-Bénard convection. In this section, we present results of uncertainty study for the Rayleigh-Bénard convection in both 2-D and 3-D. We investigate uncertainty propagation, which is triggered either by the initial data or cloud parameters.
Example 7: 2-D case with stochastic initial data. In this experiment, we choose the same initial data for the flow variables as in Section 4.2 and perturb the initial data in q v by 5%, 10%, 20% and 50%: where ν = 0.05, 0.1, 0.2 and 0.5, respectively. It should be observed that uniform perturbations in the initial conditions for q v may lead to either reduced or enhanced water vapor concentrations as compared to the deterministic simulations. Since the temperature gradient is quite small, a change in q v translates to an (almost) linear change in saturation ratio, which directly controls cloud formation. Thus, in the case of positive perturbations, a higher water vapor concentration leads to earlier cloud formation and, in addition, a higher potential temperature change since more water is available in the system. On the other hand, lower values of q v lead to a time delay in the formation of clouds, even if small convective cells are driven by the dry unstable situation. In a feedback cycle, a reduced or even delayed formation of cloud water propagates further to a weaker rain formation. Finally, the evaporation of rain water leads to a strong cooling effect of the lower layers of the domain, which also crucially depends on the amount of sedimenting rain water. These effects have to be taken into account for the evaluation of the different perturbation scenarios.
Numerical solutions for both the potential temperature θ ( Figure 17) and the expected value of the cloud drops concentration q c (Figure 18) are computed at time t = 1400s using 320 × 320 mesh cells and presented for ν = 0, 0.05, 0.1, 0.2 and 0.5. For a better comparison, we have used the same range of values for different perturbations in all of the plots.
From the results shown in Figure 17 for the potential temperature θ, one can observe two major features that increase in strength with increasing strength of perturbations. First, larger perturbations lead to more concise fine structures, that is, for the deterministic simulation (ν = 0), the variable is quite smooth, whereas for large perturbations, large gradients on a very small scale appear. This feature can also be recognized for the cloud water, that is, for the expected values E[q c ] of cloud water concentrations in Figure 18. This is probably due to the fact that even small variations in water vapor have a strong impact on cloud formation, since the activation of cloud droplet is basically a threshold process. For values closer to saturation, even small variations in vertical upward motions can trigger cloud formation. Thus, the small scale variations are more prominent in the massively perturbed scenarios. The second feature is even more striking. Increasing the perturbation in the initial water vapor distribution leads to a stronger vertical gradient in potential temperature, that is, at low levels the temperature is much smaller than in the deterministic case. This feature is due to cooling of sedimenting rain water. For simulations with a high water vapor loading, more cloud and thus more rain is formed, which is subsequently falling down into lower levels and cools the environment by evaporation. Since this process is very effective, the temperature can be reduced drastically.
Note that this feature is well-known for the real atmosphere: Falling rain can cool lower levels efficiently, so that a transition from melting rain droplets to snow can be possible for winter seasons. The efficient formation of rain also leads to a strong reduction in cloud water, since accretion can eat up cloud droplets in the lower level; see Figure 18. For reduced values of initial water vapor, processes of cloud and precipitation formation is strongly reduced. However, on average the positive perturbations dominate the statistical picture.
In Figure 19, we plot the standard deviation of the cloud drops concentration. Here, we compute the standard deviation for a function f ( which follows from the orthogonal property (6.2).
In order to investigate the influence of perturbations, the time evolution of the mean expected value per m 2 as well as the mean standard deviation per m 2 for the cloud variables are presented in Figure 20. In d-space dimensions these quantities can be computed in the following way: where N is the number of mesh cells and ∈ {v, c, r}.
These averaged quantities show the qualitative difference in the different perturbations scenarios, as already described above. For low perturbations, the difference of the expectation value is quite small and also the standard deviation remains small. For increasing perturbations, the spread is increased. As noted before, the averaged quantities are dominated by the positive perturbations, leading to (i) earlier cloud formation, (ii) thicker clouds due to more available water vapor, and (iii) to enhanced rain formation. These three features can be seen very nicely in the strongest perturbations (ν = 0.5), with a large drop in water vapor concentration accompanied by a strong increase in cloud water and earlier onset of precipitation. We would also like to note that the spread is only given by the standard deviation, whereas the actual minima (for instance, almost no cloud formation) cannot be seen directly, although these scenarios are possible.
Example 8: 2-D case with stochastic parameters. In the following experiment we study uncertainty propagation due to incomplete information about the model parameters which is a very typical problem arising in atmospheric science. We chose the same initial data for the flow and cloud variables as in Section 4.2. More precisely, we take the following initial cloud variables: Consequently, in order to investigate uncertainty propagation in the numerical solution we choose 10%, 20% and 50% perturbation of these coefficients, namely, we take where ν = 0.1, 0.2 and 0.5, respectively. The numerical solution is computed at time t = 1400s on a 320×320 mesh. Figures 21 and 22 present the potential temperature θ and the expected values of the cloud drops concentration q c , respectively, for ν = 0, 0.1, 0.2 and 0.5.
In these scenarios with perturbed cloud model parameters, the overall structures are more stable and changes are less pronounced than in Example 7. Nevertheless, the main features of the variations are obviously driven by precipitation processes since the only perturbed parameters are those determining rain processes. Again, one key feature of the perturbed scenarios is the signature of evaporating rain in lower levels of the 2-D domain. For positive perturbations, that is, larger parameters k 1 , k 2 and α, rain formation is enhanced (more rain is formed from cloud water, due to larger k 1 and k 2 ) and sedimentation is enhanced (more rain is falling downwards due to larger α). Thus, more rain water is transported downwards into subsaturated regions, which is then evaporated inducing cooling due to latent heat consumption. These positive variations again dominate the potential temperature field; see Figure 21. For the cloud water field (Figure 22), one can observe higher expected values for stronger perturbations. This is probably due to the fact that more rain is evaporated in lower levels, thus more water vapor is then available for cloud formation in upward motions of the convective cells. This redistribution of water vapor as well as the reduction of rain for negative perturbation lead to larger variations of cloud water in the 2-D domain, as can be seen in the standard deviation for the cloud drops concentration q c , as depicted in Figure 23.
The expected values as well as the standard deviation for the cloud variables vary considerably with respect to the perturbation of model parameters; see Figure 24. Our numerical experiments indicate that the standard deviation increases in time and also depends on the size of the parameter perturbation. Indeed, the larger the parameter perturbation, the higher is the standard deviation of the cloud variables. The size of the corresponding standard deviations is depicted in the right column of Figure 24. As expected, rain formation sets on earlier for large perturbations, leading also to a strong decrease in the overall water vapor in comparison to less perturbed scenarios. However, the variations in all water variables increase a lot from the onset of precipitation to later times, and even the spread in cloud water increases in contrast to the time evolution in Example 7. Generally, the spread in the mean values is smaller than in Example 7; changes in initial data produce a larger variation, that means the model physics is quite stable with respect to perturbations in rain process formulations.
Example 9: 3-D case with stochastic initial data. Similarly to Example 7, we now investigate the uncertainty quantification in the 3-D Rayleigh-Bénard convection for stochastically perturbed initial data of the cloud variables given by (7.1). The numerical solution is computed in a domain Ω = [0, 5000] × [0, 5000] × [0, 1000] m 3 , which is discretized using 50 × 50 × 50 mesh cells. In Figure 25, we show the influence of 0%, 10%, 20% and 50% perturbation on the potential temperature θ as well as the expected value of the cloud drop concentration q c . For the perturbation scenarios in the case of the 3-D moist Rayleigh-Bénard convection, one can see overall a qualitatively similar picture as for the corresponding 2-D simulations. However, due to an additional spatial direction, the 3-D structures can further change their patterns. For positive perturbations, clouds can be formed even at low vertical upward motions. The latent heat release increases the vertical motions in the convective cells, which leads to additional feedback, such as stronger cloud formation, which in turn leads to formation of larger amount of rain water. As in the 2-D case, the potential temperature distribution changes considerably at lower levels of the 3-D domain since evaporative cooling of precipitation is a dominant process. Similarly, the mean cloud water distribution is crucially changed. For positive perturbations, which are dominant on average in the beginning, more cloud water is formed and is later removed. Thus, less cloud water is available in the domain at later times. This can also be seen in the time evolution of spatially averaged variables, as shown in Figure 26, where the time evolution of the mean expected value per m 3 and the mean standard deviation per m 3 of the cloud variables are plotted. The time evolution of these variables is very similar to the 2-D case, but the 3-D scenarios show a very interesting feature, which is not available in 2-D. For larger perturbations, the pattern of the convective cells is changed. While convective cells have a quasi hexagonal shape for the deterministic simulation, they are ordered in a different way for larger perturbations. At upper levels, a more roll-like structure or even rectangular shape is formed. Thus, there is a transition of structures due to perturbations in the initial conditions. This feature has not been documented until now. A more detailed analysis of these structures is left for future studies.

Conclusion
In the present paper, we have studied uncertainty propagation in an atmospheric model that combines the Navier-Stokes equations for weakly compressible fluids (2.1) with the cloud equations (2.2). The latter has been recently proposed in [34] and is based on the so-called single moment approach considering the evolution equations for the mass concentrations of the water vapor, cloud drops and rain. Our numerical strategy is based on the stochastic Galerkin method that combines a finite-volume method for space-time discretization with a spectral approximation in the stochastic space. We point out that atmospheric flows are weakly compressible which leads to the low Mach number problem. One therefore needs to use a finite-volume method, which is accurate and efficient in the low Mach number regime; see [5,6]. To this end, we have chosen a suitable linear-nonlinear splitting between the fast and slow flow variables and the second-order IMEX discretization in time (the ARS (2,2,2) scheme) as described in Section 3. Coupling between the cloud model (2.2) and the Navier-Stokes system (2.8) is realized numerically by the second-order Strang splitting. The cloud equations are approximated in space by the finite-volume method and in time using the explicit third-order Runge-Kutta method with an enlarged stability region as explained in Section 3. Note that microscopic cloud dynamics requires a smaller time step than the flow dynamics and thus several microscopic cloud subiterations are realized within one macroscopic splitting time step, whose size is dictated by the flow dynamics. To the best of our knowledge, this is the first contribution that combines an accurate and efficient method for the weakly compressible Navier-Stokes equations with the stochastic Galerkin method for the uncertainty quantification of time evolution of the mass densities of water vapor, cloud drops and rain.
We have conducted extensive numerical benchmarking for both the deterministic and stochastic models and present the obtained numerical results in Sections 4 and 7. In the latter, we took into account the uncertainties in both initial data and cloud model parameters. Our numerical study clearly demonstrates applicability of the stochastic Galerkin method for the uncertainty quantification in complex atmospheric models. We have obtained interesting results illustrating the behavior of clouds in different perturbed scenarios and demonstrated that perturbations in the initial conditions can crucially change the time evolution of the moist Rayleigh-Bénard convection. In particular, it has been shown that for larger perturbation, on average the positive perturbations dominate the expectation values, although the standard deviation can be quite large. The main feature is a strong evaporative cooling in lower levels of the (2-D/3-D) domain due to enhance rain formation and sedimentation into low humidity levels. In the 3-D case, a change in the formed pattern can be seen, changing from hexagons to rolls/rectangles, which is quite surprising. For perturbations in parameters for rain processes, the results are also dominated by the positive part of the perturbations. Since rain processes are affected, the spread in the time evolution is increasing since the changes depend crucially on the interaction of formed rain drops with other variables. Overall, it seems that for cloud physics, the expectation values are dominated by the positive perturbations leading to a change in the distributions. This is an interesting topic for detailed studies in future. Our further goal is to extend the developed numerical method to the fully random Navier-Stokes-cloud system by considering random weakly compressible Navier-Stokes equations. We are also interested in considering different random effects, such as initial data, boundary data and model parameters simultaneously, which would require a multivariate stochastic Galerkin method. This will allow one to quantify more precisely the propagation of small scale stochastic errors initiated at cloud scales to macroscopic scales of flow dynamics.
Lukáčová-Medvid'ová and B. Wiebe would like to thank L. Yelash (University of Mainz) for fruitful discussions. P. Spichtinger would like to thank P. Reutter (University of Mainz) for fruitful discussions on the results of the 2-D/3-D simulations.
9. Appendix A: Closure for single moment schemes The number concentration of rain drops can be approximated by a function of the respective mass concentration n r = f (q r , c r ). Since we implicitly assume that rain drops are distributed according to a size distribution, this approach should be used for mimicking the shape of the distribution in a proper way. If we use a constant mean mass of rain drops, the function will be a simple linear relation n r = 1 mr q r . We extend this approach and propose the following nonlinear relation: n r = c r · q γ r , 0 < γ ≤ 1.
Using this approach, one can replace the quantity n r in the processes related to rain drop number concentration. For the simple case of a constant mean mass m r , we can determine the constants as c r = m −1 r and γ = 1. This approach would be meaningful for the case of a symmetric size distribution of rain droplets centered around the mean mass. However, it is well-known that size distributions of rain are usually skew to larger sizes and thus a linear relation is inappropriate. For sizes of rain drops, an exponential distribution is often assumed (see [25]), namely: Using these relations, one can derive the following function for the number concentration n r : n r = B We stress that c r is, in fact, a function of the air density ρ, that is, c r = c r0 · ρ − 3 4 .

Appendix B: Explicit formulation of the cloud equations
We present the equations of microphysical processes in an explicit way as they are used in our numerical experiments. In Tables 3 and 4, we present physical constants and model parameters with their values used in our numerical simulations.  Table 4. Model parameters, [34].