Stochastic Galerkin method for cloud simulation. Part II: a fully random Navier-Stokes-cloud model

This paper is a continuation of the work presented in [Chertock et al., Math. Cli. Weather Forecast. 5, 1 (2019), 65--106]. We study uncertainty propagation in warm cloud dynamics of weakly compressible fluids. The mathematical model is governed by a multiscale system of PDEs in which the macroscopic fluid dynamics is described by a weakly compressible Navier-Stokes system and the microscopic cloud dynamics is modeled by a convection-diffusion-reaction system. In order to quantify uncertainties present in the system, we derive and implement a generalized polynomial chaos stochastic Galerkin method. Unlike the first part of this work, where we restricted our consideration to the partially stochastic case in which the uncertainties were only present in the cloud physics equations, we now study a fully random Navier-Stokes-cloud model in which we include randomness in the macroscopic fluid dynamics as well. We conduct a series of numerical experiments illustrating the accuracy and efficiency of the developed approach.


Introduction
In this paper we continue the study of generalized polynomial chaos (gPC) stochastic Galerkin methods for multiscale cloud dynamics with uncertainty. The gPC expansion based methods, such as stochastic Galerkin [3,5,9,11,12,15] and stochastic collocation [6,13,14], are well-known techniques for uncertainty quantification in physical and engineering applications. In these methods the solution is approximated by a truncated generalized Fourier series in terms of orthogonal polynomials corresponding to the underlying probability density function. The stochastic collocation method belongs to the class of non-intrusive methods in which the original deterministic system is solved at certain collocation points in the stochastic space and then the gPC coefficients are obtained by using a polynomial interpolation and a suitable numerical quadrature. The stochastic Galerkin method on the other hand is an intrusive method in which the Galerkin projection is employed after substituting the gPC expansions into the original system. This leads to a larger coupled system of deterministic PDEs for the gPC coefficients for which an accurate and efficient numerical method is to be developed. model. The numerical method for the latter is presented in §3, which is followed by numerical experiments in §4. We demonstrate the experimental convergence of the developed method as well as its applicability for uncertainty quantification in atmospheric flows through well-known meteorological benchmark tests of a rising warm and moist bubble and Rayleigh-Bénard convection.

(2.2)
Here, t is the time variable, x = (x 1 , . . . , x d ) ∈ R d is the space vector, ρ is the density, u = (u 1 , . . . , u d ) 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. Furthermore, e d = e 3 = (0, 0, 1) and e d = e 2 = (0, 1) in the three-dimensional (3-D) and two-dimensional (2-D) cases, respectively. The cloud variables representing the mass concentration of water vapor, cloud droplets and rain drops, q v , q c and q r , respectively, are given by q = mass of the respective phase mass of dry air for ∈ {v, c, r}.
The terms C and E represent phase changes between vapor and cloud water (droplets), and A 1 and A 2 represent collision processes, which lead to the formation of large droplets and thus precipitation.
Note that the systems (2.1) and (2.2) are coupled through the source term S θ , which represents the impact of phase changes and is given by For a detailed description of S θ and the terms E, C, A 1 and A 2 see [3]. The temperature T can be obtained from the moist adiabatic ideal gas equation where p 0 = 10 5 Pa is the reference pressure at the 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, that is when q v = q c = q r = 0, R m reduces to R, S θ = 0 and the moist ideal gas equation as well as the moist equation of state become their dry analogue.
Solving the Navier-Stokes equations (2.1) in a weakly compressible regime is known to trigger 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.3) Let p , ρ , u , θ and (ρθ) stand for the corresponding perturbations of the equilibrium state, then p =p + p , ρ =ρ + ρ , θ =θ + θ , u = u , ρθ =ρθ +ρθ + ρ θ + ρ θ = ρθ + (ρθ) .
The pressure perturbation p is derived from (2.3) and (2.4) using the following Taylor expansion which results in Thus, a numerically preferable perturbation formulation of the Navier-Stokes equations (2.1) reads (2.5) Meteorological applications typically inherit several sources of uncertainty, such as model parameters, initial and boundary conditions. Consequently, stochastic models need to be designed to analyze the influence of uncertainties on the fluid and cloud dynamics. In general, there are different ways to represent and take into account model uncertainty. In this paper, we choose a widely used approach in which the uncertainty is described by random fields. To this end we define an abstract probability space (Γ, Σ, P ) and denote by ω an event ω ∈ Γ. We assume that the initial data depend on ω, that is, for the fluid variables and for the cloud variables. Consequently, the solution at later time will also depend on ω, that is, we will have ρ (x, t, ω), (ρu)(x, t, ω), (ρθ) (x, t, ω) and (ρq )(x, t, ω) for ∈ {v, c, r}.
Remark 2.1. The system parameters and boundary conditions could also depend on the random variable. Some of that cases were considered in our previous work [3], but in this paper we restrict our attention to the situation in which randomness arises in the initial data only.

Numerical scheme
In this section, we describe a gPC-SG method for the coupled system (2.5), (2.2). First, in §3.1 we derive a system for the gPC coefficients and then in §3.2 we describe a method used to numerically solve the resulting system.

3.1.
System of gPC coefficients. In the gPC-SG setup, the solution is sought in the form of polynomial expansions where Φ k (ω), k = 0, . . . , M , are polynomials of degree k that are orthogonal with respect to the probability density function µ(ω). Assuming that Γ is a compact metric event space, the corresponding Riemann integrals can be defined (see [10]), and then the orthogonality property implies where δ kk is the Kronecker symbol and c k are constants depending on the probability density function µ.
In this work, we will focus on two distributions that are important for meteorological applications: 1], which corresponds to the Legendre polynomials if n is even, , if n is odd, which are orthogonal with respect to the probability density function µ(ω) = 1 /2 and the constants in where µ H and σ H are the mean value and the standard deviation, respectively. One can show that the Hermite polynomials Φ k are orthogonal with respect to µ(ω) = and the constants in (3.3) are c k = k! (3.5) We use gPC expansions for the source term S θ on the right-hand side (RHS) of (2.5), for the terms on the RHS of (2.2) and for the raindrop fall velocity, Substituting (3.1), (3.6) into (2.5) and (3.2), (3.7), (3.8) into (2.2), applying the Galerkin projection and using the orthogonality property (3.3) yield the following deterministic system consisting of (d + 2)(M + 1) equations for the gPC coefficients of the fluid variables and 3(M + 1) equations for the gPC coefficients of the cloud variables for k = 0, . . . , M . It should be pointed that for the pressure term in the momentum equation in (3.9) the expected values for q v , q c and q r are used for computing R m , which gives In this way we keep the linear formulation of the pressure as in the deterministic case. In what follows, we explain how the expansion coefficients in (3.9) and (3.10) are obtained.

• Discrete transform (DT)
The discrete transform starts with the expansion of a function f in the stochastic space (3.11) and by using the orthogonality property (3.3) ends up with the expansion coefficients We approximate the above integral using an appropriate Gaussian quadrature rule. We distinguish between the two cases considered in this paper-Legendre and Hermite polynomials.
(1) For Legendre polynomials the coefficients c k are determined in (3.4), Γ = [−1, 1], and the expansion coefficients in (3.12) are given by Approximating the above integrals using the Gauss-Legendre quadrature leads to where ω n and β n are the Gauss-Legendre nodes and weights, respectively.
(2) For Hermite polynomials the coefficients c k are determined in (3.5), Γ = (−∞, ∞), and the expansion coefficients in (3.12) are given by Approximating the above integral using the Gauss-Hermite quadrature leads to where ω n and β n are the Gauss-Hermite nodes and weights, respectively.
• Inverse discrete transform (IDT) The inverse discrete transform maps the expansion coefficients To this end we simply compute the point values of f using the gPC expansion (3.11): (3.15) Remark 3.1. The number of quadrature points N can be chosen equal to the number of expansion coefficients M , or even higher for a more accurate approximation.
Remark 3.2. We stress that the quadrature weights β n and the values Φ k (ω n ), 0 ≤ k ≤ M , 0 ≤ n ≤ N , which are used in (3.13)-(3.15), can be pre-computed for the code efficiency.
The coefficients {( q ) k (x, t)} M k=0 are computed in the following way: where the DT and IDT operators are defined in (3.13)-(3.15). Additionally, the coefficients (3.10) are obtained through the DT and IDT in a similar manner as in (3.16). Here, the nonlinear Navier-Stokes advection coefficients and the nonlinear diffusion coefficients Finally, the nonlinear cloud dynamics advection coefficients and the nonlinear diffusion coefficients 3.2. Discretization of the gPC system. In this section we describe the numerical methods used to solve the resulting gPC system (3.9), (3.10). To this end we denote by w := ρ , ρu, (ρθ) and w q := ρq v , ρq c , ρq r the solution vectors of (3.9) and (3.10), respectively. Here, the underline (·) denotes the vector of the respective coefficients. For instance, for the solution coefficients we have Then, the coupled system can be written as where F and F q are convective fluxes and D, R and D q , R q denote the diffusion and reaction operators of the respective systems. They are given by where the respective components of the above vectors were defined in (3.17)-(3.20). It should be observed that (3.21)-(3.24) is a deterministic system for the expansion coefficients. This system has the same structure as the deterministic system studied in [3]. Therefore, one can directly apply the finite volume method from [3] for the spatial discretization of the system (3.21)-(3.24) taking into account that the number of equations has increased by a factor of M + 1 and that the DT and IDT should be applied for each evaluation of the nonlinear terms on the RHS of (3.21) and (3.22).
Concerning the time discretization, we also follow the approach presented in [3] and implement a Strang splitting approach between the Navier-Stokes and cloud dynamics systems. We then approximate the cloud dynamics (3.22), (3.24) with a third-order large stability domain explicit Runge-Kutta method (DUMKA3 see [7,8]) and the fluid dynamics (3.21), (3.23) with the IMEX AP ARS(2,2,2) method from [1]. For the latter, we split the Navier-Stokes operators in (3.23) into two parts and then define the stiff linear operator L := −∇ · F L ( w) + R L ( w) and the nonstiff nonlinear operator , which are treated implicitly and explicitely, respectively.
Finally, we note that as in [3], the time steps for both the Navier-Stokes and cloud dynamics splitting substeps are chosen adaptively at every time level.

Numerical experiments
In this section, we present experimental results for the fully random Navier-Stokes-cloud model (2.5), (2.2). In Examples 4.1 and 4.2, we investigate the experimental convergence of our numerical scheme using the wellknown meteorological benchmark describing the free convection of a smooth warm and moist air bubble; see, e.g., [2,4]. In Example 4.1, we demonstrate the spatio-temporal convergence as well as the convergence in the stochastic space for the case in which the initial vapor concentration q v is perturbed by 10% which is realized with a uniform distribution of the randomness. In Example 4.2, we investigate the convergence for the same setup as in Example 4.1 but with normally distributed randomness. Since we use the same numerical method for the space and time discretization as in Example 4.1, in Example 4.2 we just investigate the convergence in the stochastic space. In Examples 4.3 and 4.4, we present the results of the uncertainty study for the Rayleigh-Bénard convection in both 2-D and 3-D. We also compare the results obtained in this work for the fully random model with the deterministic one, in which both the Navier-Stokes euqations (2.5) and the cloud equations (2.2) are deterministic, and with the semi-random one, in which the deterministic Navier-Stokes equations are coupled with the random cloud dynamics (this semi-random model was studied in the first part of this work in [3]). We note that the parameters used in the numerical experiments presented below are slightly different from the one used in [3]; however, the main qualitative features remain the same. In both Rayleigh-Bénard experiments (Examples 4.3 and 4.4), we investigate uncertainty propagation, which is triggered by the initial data of the water vapor concentration q v which we perturbed uniformly. A comparison with a normally distributed initial perturbation or even perturbations of certain parameters is beyond the scope of this work and is left for future study.
Example 4.1 (Stochastic initial data with uniformly distributed perturbation).
In this experiment, we simulate free convection of a smooth warm and moist air bubble in 2-D. Due to the shear friction with the surrounding air at the warm/cold air interface, the warm air bubble rises and deforms axisymmetrically and gradually forms a mushroom-like shape. The bubble is placed at (2500 m, 2000 m) in a domain Ω = [0, 5000]×[0, 5000] m 2 . We consider a 10% perturbation of the initial water vapor concentration. This is realized through the following initial conditions in the case of a uniformly distributed randomness for the cloud variables: and for the Navier-Stokes variables: Additionally, we setθ = 285 K, p 0 =p = 10 5 Pa and with c p = 1005 J /(kg·K), c v = 718 J /(kg·K) and γ = c p /c v . We start here with nonzero values for the cloud drops concentration q c and the rain concentration q r to avoid values close to the machine precision since the main purpose of the test is the convergence study. Furthermore, we apply the no-slip boundary conditions for the velocities and homogeneous Neumann boundary conditions for the remaining variables, that is, ∇ρ · n = 0, ∇(ρθ) · n = 0 and ∇(ρq ) · n = 0, ∈ {v, c, r}.
In Figure 1, we depict the expected values of the potential temperature θ and the cloud variables q v , q c and q r , computed using a 160 × 160 uniform mesh at time 200s with M = L = 3. For comparison purposes, in Figure 2 we show the potential temperature θ and the water vapor concentration q v computed using the deterministic Navier-Stokes-cloud model and the potential temperature θ and the expected values of the water vapor concentration q v computed using the semi-random Navier-Stokes-cloud model. Note that for a better comparison, we have used the same vertical scales for presenting the results in Figures 1 and 2. It can be observed that the fully random results are more smeared compared to the deterministic ones and that in the fully random experiment no additional vortices beneath the bubble have been developed and the results are slight variations of the deterministic ones, which is to be expected. In order to investigate the appearance of the vortices in the semi-random case (see the second row Figure 2), we depict in Figure 3 the semi-random results obtained for the smaller initial water vapor perturbation taken as 1%, 5% and 7%. One can clearly see, the vortices develop gradually with higher perturbation and that the 1% perturbation results are very close to the deterministic ones. Thus, the vortex features of the solutions obtained with the semi-random model seem to result from the missing feedback to the dynamics of the fluid and are not a defect of the numerical method. We also note that in the semi-random model, the energy conservation is (slightly) violated. This might lead to the differences in the simulations, since the dynamics is then driven by the averaged latent heat release and not by the one in the realization. In Figure 4, we present the spatio-temporal convergence study for the expected values of the cloud and flow variables at the time t = 10s. We compute the solutions on different N ×N uniform meshes with M = L = 3. As in the deterministic case presented in [3], one can clearly see a second-order convergence for the studied fully random model.
The stochastic convergence studies are presented in Figures 5 and 6 for the cloud and Navier-Stokes variables, respectively, at time t = 10s using a 160 × 160 uniform mesh. We compute the difference between the approximate solutions with different numbers of modes M and L = M and the reference solution obtained with 20 stochastic modes and L = 19. One can observe a spectral convergence with an approximate rate of e −0.3M . One can also see that the error of the rain drops in Figure 6 (right) basically stays constant at some point because in this case it approaches the machine precision. In this experiment, we demonstrate that the convergence of the stochastic Galerkin method for the fully stochastic model does not depend on the choice of the distribution of the randomness. For this purpose, we choose the same initial conditions as in Example 4.1, but this time with a normally N (0, 1) distributed perturbation. In Figure 7, we compare the solutions (the potential temperature θ and the water vapor concentration q v ) computed using the deterministic, semi-random and fully random Navier-Stokes-cloud models. For a better comparison, we have used the same vertical scales for presenting the results. As in the previous example, one can observe that in the fully random experiment no additional vortices beneath the bubble have been developed and the results are slight variations of the deterministic ones. Thus, the vortex features of the semi-random results are independent of the distributions of the initial perturbation and caused by the missing feedback to the dynamics of the fluid.
Next, we investigate the influence of the choice of distribution for the initial perturbation. In Figure 8, we depict the he cloud drops concentration q c computed using the fully random model with the initial normally and uniformly distributed perturbations; the latter one was computed in Example 4.1. For a better comparison, we have used the same vertical scales for presenting the results. Since the initial perturbation is rather small, the results look very alike. In general, both simulations smear the boundaries of the bubble. However, the smearing with the normal distribution is not as strong as with the uniform distribution. This effect is due to the concentrated shape of the normal distribution around the expected value; thus, the different realizations are closer to the averaged potential temperature as a feedback to the energy equation.
The convergence studies in the stochastic space are presented in Figures 9 and 10 for the cloud and Navier-Stokes variables, respectively, at time t = 10s using a 160 × 160 uniform mesh. We computed the difference between the approximate solutions with different numbers of modes M and L = M and the reference solution obtained with 12 stochastic modes and L = 11. As in the case with a uniform distribution studied in Example 4.1, one can observe a spectral convergence with an approximate rate of e −0.3M . This demonstrates that the experimental convergence rate is independent of the chosen distribution.      We consider a 2-D stochastic Rayleigh-Bénard convection simulated on a domain Ω = [0, 5000] × [0, 1000] m 2 that has been discretized using a uniform 160 × 160 mesh. The initial data for the cloud variables are ( q r ) 0 (x, 0) = 10 −6 (θ (x, 0)) + , ( q r ) k (x, 0) = 0 for 1 ≤ k ≤ M, and for the Navier-Stokes variables are .   We implement the following Dirichlet boundary conditions for the potential temperature: θ(x 2 = 0) = 284 K and θ(x 2 = 1000) = 283 K, as well as the periodic boundary conditions for all of the variables in the horizontal direction, no-slip boundary conditions for the velocities at the vertical boundaries, and zero Neumann conditions for the remaining variables in the vertical direction, that is, ∇ρ · n = 0. Also, these boundary conditions have to be projected onto the stochastic space. The projections of the periodic, no-slip and Neumann boundary conditions are straightforward and lead to the same conditions as in the deterministic case for all of the expansion coefficients of the respective variable. Here, we briefly explain how the projection of the Dirichlet boundary conditions works. We implement the Dirichlet boundary conditions for the potential temperature using ρθ(x, t, ω) = (ρθ) (x, t, ω) + ρθ(x). Rearranging and inserting the expansion for ρ (x, t, ω) and (ρθ) (x, t, ω) At θ is constant at the boundary, applying the projection leads to and analogously for x 2 = 1000.
In Figures 11 and 12, we present snapshots of the expected values and standard deviations of the potential temperature and the cloud variables at times t = 1000 and 6000s, respectively. Additionally, in Figure 13, we plot the the differences between the expected value of the water vapor concentration and the saturation mixing ratio E[q v ] − q * at the same times. As one can observe, the potential temperature exhibits a strong vertical gradient at time t = 1000s. Similarly to the deterministic and semi-random cases, at a later time t = 6000s, supersaturated regions are formed in the rolls where the convection takes place leading to the overall roll-like cloud flow structure. Figure 11. Example 4.3: Expected value and standard deviation of the potential temperature θ at t = 1000s and 6000s.
In Figures 14 and 15, we present the time evolution of the mean expected value per m 2 as well as the mean standard deviation per m 2 for the potential temperature and the cloud variables. In d space dimensions these quantities can be computed for uniformly distributed perturbations in the following way: where N d is the number of mesh cells and ∈ {v, c, r}. We compare the solutions using 0% (purely deterministic model) and 10% of perturbation of the initial data in q v , where for 10% of perturbation the solutions are added in both fully-and semi-random Navier-Stokes-cloud models. The time evolution of the averaged quantities clearly shows the differences between the semi-random and fully random models. In all shown cases (including the purely deterministic one), the time evolution starts with cloud formation and thus increase of cloud water on the expense of water vapor and also latent heat release (increase of θ). However, for the semi-random model the rain formation starts earlier than in the deterministic and fully random simulations. Since rain is falling into subsaturated regions which induces evaporation, this leads to a different time evolution in all variables. Generally, we observe a much stronger cooling effect of the system due to evaporation of rain in the semi-random model than in the other simulations. This is probably due to the use of the expected values of the terms for phase changes in the energy equation. Although the general qualitative behavior in the time evolution of the expected values is quite similar, the absolute values differ quite substantially. The same is true for the standard deviations of the cloud variables q v q c and q r as shown in the right column of Figure 15; the variation in the water variables (and also in θ) is much larger for the fully random model than for the semi-random one. This is also reasonable, since the fully random model can capture the correct feedback of the latent heat release in the phase changes for the "different realizations", whereas the semi-random model only feeds back the averaged potential temperature, leading to a smaller variability. Some examples of the expected values and the related standard deviations for θ, q c and the super/subsaturation (in terms of E[q v ] − q * ) are shown in Figures 11-13. for the potential temperature θ per m 2 using 0% (purely deterministic case) and 10% perturbation of the initial data in qv, where the latter was simulated using both fully-and semi-random models. In the final example, we consider a 3-D stochastic Rayleigh-Bénard convection. The initial data for the cloud variables are ( q r ) 0 (x, 0) = 10 −6 (θ (x, 0)) + , ( q r ) k (x, 0) = 0 for 1 ≤ k ≤ M, with ν = 0, 0.1, 0.2 or 0.5, which correspond to 0% (pure deterministic case), 10%, 20% or 50% perturbation of the initial water vapor concentration. For the Navier-Stokes variables we take purely deterministic initial and standard deviation (right column) using 0% (purely deterministic case) and 10% perturbation of the initial data in qv, where the latter was simulated using both fully-and semirandom models.
data, which in terms of their expansion coefficients read as In Figures 16-18, we present the influence of the 10%, 20% and 50% initial water vapor perturbation on the expected values of the potential temperature, cloud droplets and rain drops concentration at times t = 1000s and 6000s. The influence on the supersaturated and subsaturated regions is highlighted as a 2-D slice along x 1 = 3000 in Figure 19, where we depict the difference between the expected water vapor concentration and the saturation mixing ratio. For a better comparison, we have used the same vertical scales in all of the plots. Here, one can clearly observe a different behavior compared with the semi-random case. The vertical gradient of the potential temperature increases as the size of the initial perturbations increases (see Figure  16), while the pattern of the developed convection cells is similar for different perturbations. The latent heat release increases the vertical motions in the convective cells, which leads to additional feedback, such as stronger and more cloud formation (see Figure 17), which in turn leads to the formation of a much larger amount of rain water, especially at a later time t = 6000s (see Figure 18). At the time t = 1000s one can see that the roll-like structure of the clouds in the deterministic case (that is, with 0% perturbation) again end up in a more cell-like structure in the initially perturbed cases.
In Figures 20 and 21, we show the time evolution of the mean expected value per m 3 as well as the mean standard deviation per m 3 for the potential temperature and cloud variables in the cases with 0% (purely deterministic), 10%, 20% and 50% perturbation of the initial water vapor concentration. For increasing perturbations, the spread is increased, mostly for the water vapor concentration q v and the rain concentration q r . 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) enhanced rain formation. These three features can be clearly seen in the case of the largest initial perturbation (50%), where a large spread in water vapor concentration is accompanied by a strong increase in cloud water and an earlier onset of strong precipitation. Due to the strong rain formation the cloud concentration decreases when the perturbation size increases and also the amount of supersaturated regions is much smaller as can be observed in Figure  19 which leads to less new formation of clouds. 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. Overall, one can see that the time evolution for the deterministic simulation as well as for perturbations with 10% and 20% behave quite similarly and the averaged quantities follow closely the same evolution, although the standard deviations increase quite substantially. However, for larger perturbations (50%), the time evolution of the expected values of q c and q r is strongly disturbed and shows large deviations from the other simulations. This can also be seen in the 3-D panels at the later time t = 6000s.

Conclusion
This paper is a continuation of our previous study on uncertainty propagation in atmospheric flows containing phase changes. In particular, we consider warm cloud dynamics of weakly compressible fluids. This model consists of a multiscale system of PDEs in which the macroscopic dynamics of the fluid is described by a weakly compressible Navier-Stokes system and the microscopic cloud dynamics is described by a system of convection-diffusion-reaction equations. We have extended the gPC-SG method from [3], where we considered a semi-random model with the deterministic macroscopic dynamics of the fluid coupled with the random microscopic cloud dynamics, to the case of a fully random multiscale system. To this end, we have first derived a system for the gPC coefficients and then presented a method we have used to numerically solve the resulting system. The latter is an extension of the numerical method developed in [3] and approximates the gPC coefficients for the dynamics of the fluid by an IMEX AP finite volume method and the gPC coefficients for the cloud dynamics by an explicit finite volume method with an enlarged stability region.
The aim of this work is to demonstrate the applicability, accuracy and efficiency of the gPC-SG method for atmospheric flows. Comprehensive studies of uncertainty propagation in these models considering different perturbation scenarios are left for a future work. Additionally, we will investigate the accuracy and performance of different uncertainty quantification methods, for instance, stochastic Galerkin, stochastic collocation and Monte Carlo method, in a review paper. Here, we have focused on numerical convergence and benchmark experiments as well as comparison with the results of the previous semi-random model presented in [3]. We have demonstrated that the gPC-SG method for the fully random model preserves the secondorder spatial experimental convergence rate when the time increments are chosen according to the time step restriction and additionally exhibited an experimental exponential convergence rate in the stochastic space. This experimental convergence rate has been observed for both perturbation scenarios: uniform and normal distribution of the initial data perturbation. Additionally, we have studied the numerical solutions of the fully random cloud model for both the 2-D and 3-D Rayleigh-Bénard convection. By illustrating the behavior of clouds in different perturbed scenarios, we have demonstrated that perturbations of the initial conditions of cloud variables can crucially change the time evolution. The results have also exhibited a clear difference of the solutions of the semi-and fully random models in both the 2-D and 3-D Rayleigh-Bénard convection, which indicates that initial perturbations of cloud variables propagate to the Navier-Stokes equations and have a significant effect on the fluid variables. Our numerical study clearly demonstrates the applicability of the stochastic Galerkin method for the uncertainty quantification in complex atmospheric models and paves the path for more extensive practically relevant numerical studies.