University of Warsaw Lagrangian Cloud Model (UWLCM) 1.0: a modern large-eddy simulation tool for warm cloud modeling with Lagrangian microphysics

A new anelastic large-eddy simulation (LES) model with an Eulerian dynamical core and Lagrangian particle-based microphysics is presented. The dynamical core uses the multidimensional positive-definite advection transport algorithm (MPDATA) advection scheme and the generalized conjugate residual pressure solver, whereas the microphysics scheme is based on the super-droplet method. Algorithms for coupling of Lagrangian microphysics with Eulerian dynamics are presented, including spatial and temporal discretizations and a condensation substepping algorithm. The model is free of numerical diffusion in the droplet size spectrum. Activation of droplets is modeled explicitly, making the model less sensitive to local supersaturation maxima than models in which activation is parameterized. Simulations of a drizzling marine stratocumulus give results in agreement with other LES models. It is shown that in the super-droplet method a relatively low number of computational particles is sufficient to obtain correct averaged properties of a cloud, but condensation and collision–coalescence have to be modeled with a time step of the order of 0.1 s. Such short time steps are achieved by substepping, as the model time step is typically around 1 s. Simulations with and without an explicit subgrid-scale turbulence model are compared. Effects of modeling subgrid-scale motion of super-droplets are investigated. The model achieves high computational performance by using graphics processing unit (GPU) accelerators.

1 Introduction

Over the last decade, Lagrangian particle-based cloud microphysics schemes have been drawing increasing attention. They are similar to Eulerian bin schemes in that they explicitly model the size spectrum of droplets and explicitly resolve microphysical processes, but have a number of advantages over them (Grabowski et al., 2018b). One of the advantages is that Lagrangian schemes have no numerical diffusion in the spectrum of droplet sizes. Several Lagrangian schemes for warm cloud microphysics have been developed thus far (Andrejczuk et al., 2008; Shima et al., 2009; Riechelmann et al., 2012). Arguably, the most important difference between these schemes is in the way collision–coalescence is modeled. The coalescence algorithm used in the super-droplet method (SDM) of Shima et al. (2009) seems to be the most promising, as it was found to be the most accurate of the coalescence algorithms used in various Lagrangian microphysics schemes (Unterstrasser et al., 2017, where it is called the “all-or-nothing” algorithm). A numerical implementation of the SDM is a major part of the libcloudph++ library (Arabas et al., 2015) developed by the cloud modeling group at the University of Warsaw.

In this paper, we document the development of a new large-eddy simulation (LES) model called the University of Warsaw Lagrangian Cloud Model (UWLCM). It is an anelastic model with a finite-difference Eulerian dynamical core and a Lagrangian microphysics scheme. The Lipps–Hemler anelastic approximation (Lipps and Hemler, 1982) is used, which is applicable to a wide range of atmospheric flows (Klein et al., 2010; Smolarkiewicz, 2011). The dynamical core is implemented using the libmpdata++ software library (Jaruga et al., 2015) also developed by the cloud modeling group at the University of Warsaw. libmpdata++ is a collection of solvers for the generalized transport equation. In libmpdata++, advection is modeled using the multidimensional positive-definite advection transport algorithm (MPDATA) – see Smolarkiewicz (2006) for a recent review. Liquid water is modeled with the SDM implemented in libcloudph++. We do not assume any artificial categorization of liquid water particles. Consequently, all particles, i.e., humidified aerosols, cloud droplets and rain drops, evolve according to the same set of basic equations.

One of the key reasons for developing a new model is to use a modern software development approach. The model code is written in the C++ programming language and makes use of many mature libraries available in that language (e.g., Blitz++, Boost and Thrust). The code is open-source and under a version-control system. A set of automated tests greatly helps in ensuring the correctness of the model. The automated tests include a 2-D moist thermal simulation, a 2-D kinematic stratocumulus simulation and a test of different combinations of model options. Moreover, modeling of physical processes, e.g., condensation, advection, coalescence, sedimentation, is tested separately by the libmpdata++ and libcloudph++ test suites. UWLCM makes efficient use of modern computers that have both central processing units (CPUs) and graphics processing units (GPUs). The Eulerian computations of the dynamical core are carried out on CPUs and, simultaneously, Lagrangian microphysical computations are undertaken on GPUs. However, it is also possible to run the Lagrangian microphysics on CPUs.

Some results obtained using earlier versions of UWLCM have already been published. In Grabowski et al. (2018a), UWLCM was used to model a 2-D moist thermal and in Grabowski et al. (2018b), an idealized 3-D cumulus cloud was modeled. Here, we present simulations of a drizzling marine stratocumulus using the second Dynamics and Chemistry of Marine Stratocumulus (DYCOMS-II) second research flight (RF02) setup. UWLCM results are compared with 11 LES models that took part in the Ackerman et al. (2009) intercomparison. Sensitivity of UWLCM to the parameters of the microphysics scheme and to the description of the subgrid-scale (SGS) turbulence is studied. It is of particular interest how much drizzle a LES model with Lagrangian microphysics produces, compared to models with bin or bulk microphysics that took part in the intercomparison. To our knowledge, LES simulations with warm cloud Lagrangian microphysics were used to study drizzling stratocumulus only by Andrejczuk et al. (2008, 2010). This type of model has more often been employed to study cumulus clouds (Riechelmann et al., 2012; Naumann and Seifert, 2015; Arabas and Shima, 2013; Hoffmann et al., 2015, 2017).

Section 2 presents the governing equations of the model, and  Sect. 3 describes the numerical algorithms for solving these equations. The stratocumulus simulation results are discussed in Sect. 4. Section 5 contains a summary and planned developments of the model. A list of symbols used and their definitions are given in Appendix A, Appendix B compares two substepping algorithms for condensation and Appendix C contains a brief description of the software implementation of the model.

2 Governing equations

2.1 Eulerian variables

Eulerian prognostic variables of the model are the potential temperature θ, the water vapor mixing ratio q_v and the air velocity u. Equations governing the time evolution of these variables are obtained via the Lipps–Hemler approximation, which relies on the assumption that the atmosphere does not depart far from some stationary state, called the “reference state” (Lipps and Hemler, 1982). The reference state is assumed here to be a dry, hydrostatically balanced state with constant stability S^r. S^r is equal to the average stability of the sounding used to initialize a simulation. Surface density and pressure of the reference state are equal to those of the initial sounding. Vertical profiles of potential temperature and of density of dry air in the reference state are (Clark and Farley, 1984) as follows:

$$\begin{array}{}\text{(1)}& {\displaystyle }{\mathit{\theta }}^r\left(z\right)={\mathit{\theta }}_{\mathrm{v}}^{\mathrm{0}}\exp \left(S^{r}z\right),\text{(2)}& {\displaystyle }\begin{array}{rl}{\mathit{\rho }}_{\mathrm{d}}^r\left(z\right)=& {\mathit{\rho }}^{\mathrm{0}}\exp (-S^{r}z)\\ & {\left[\mathrm{1}-{\displaystyle \frac{g}{c_{\mathrm{pd}}{S}^{\mathrm{r}}{\mathit{\theta }}_{\mathrm{v}}^{\mathrm{0}}}}\left(\mathrm{1}-\exp (-S^{r}z)\right)\right]}^{(c_{\mathrm{pd}}/R_{\mathrm{d}})-\mathrm{1}},\end{array}\end{array}$$

where ${\mathit{\theta }}_{\mathrm{v}}^{\mathrm{0}}$ and ρ⁰ are values of the virtual potential temperature and of the air density taken from the initial sounding at the ground level, respectively. An auxiliary “environmental state” is introduced to increase accuracy of numerical calculations (Smolarkiewicz et al., 2014, 2019). It is a hydrostatically balanced moist state with stationary profiles θ^e(z), p^e(z), T^e(z), $q_{\mathrm{v}}^{\mathrm{e}}\left(z\right)$ and $q_{\mathrm{l}}^{\mathrm{e}}\left(z\right)$ calculated from the initial sounding. If the initial sounding is supersaturated, all supersaturation is assumed to be condensed in the environmental state.

The set of anelastic Lipps–Hemler equations (Lipps and Hemler, 1982; Grabowski and Smolarkiewicz, 1996; Clark and Farley, 1984) that govern the time evolution of the Eulerian prognostic variables is as follows:

$$\begin{array}{}\text{(3)}& {\displaystyle }{D}_t\mathit{u}& {\displaystyle }=-\mathrm{\nabla }\mathit{\pi }+\mathit{k}B+{\mathit{F}}_u+{\mathcal{D}}_u,\text{(4)}& {\displaystyle }{D}_t\mathit{\theta }& {\displaystyle }={\displaystyle \frac{{\mathit{\theta }}^{\mathrm{e}}}{T^{\mathrm{e}}}}\left({\displaystyle \frac{l_{\mathrm{v}}}{c_{\mathrm{pd}}}}C\right)+F_{\mathit{\theta }}+{\mathcal{D}}_{\mathit{\theta }},\text{(5)}& {\displaystyle }{D}_t{q}_{\mathrm{v}}& {\displaystyle }=-C+F_{q_{\mathrm{v}}}+{\mathcal{D}}_{q_{\mathrm{v}}}.\end{array}$$

Here D_t denotes the material derivative, $D_t={\partial }_t+\mathit{u}\cdot \mathrm{\nabla }$, and π is the normalized pressure perturbation. Following Grabowski and Smolarkiewicz (1996), buoyancy is defined as

$$\begin{array}{}\text{(6)}& B=g\left[{\displaystyle \frac{\mathit{\theta }-{\mathit{\theta }}^{\mathrm{e}}}{{\mathit{\theta }}^{\mathrm{r}}}}+\mathit{ϵ}\left(q_{\mathrm{v}}-q_{\mathrm{v}}^{\mathrm{e}}\right)-\left(q_{\mathrm{l}}-q_{\mathrm{l}}^{\mathrm{e}}\right)\right].\end{array}$$

The condensation rate C in Eqs. (4) and (5) and the liquid water mixing ratio q_l in Eq. (6) come from the Lagrangian microphysics scheme. The terms F_* represent a total forcing due to surface fluxes, radiative heating/cooling, large-scale subsidence and absorbers, whereas the terms 𝒟_* represent contributions from a SGS turbulence model. The dry-air density is assumed to be equal to the reference state density profile ${\mathit{\rho }}_{\mathrm{d}}^{\mathrm{r}}$ and, characteristically for the anelastic approximation, the dry-air density at given position does not change with time: ${\partial }_t{\mathit{\rho }}_{\mathrm{d}}^{\mathrm{r}}=\mathrm{0}$. By putting ${\partial }_t{\mathit{\rho }}_{\mathrm{d}}^{\mathrm{r}}=\mathrm{0}$ into the continuity equation, the following constraint on the velocity field is obtained:

$$\begin{array}{}\text{(7)}& \mathrm{\nabla }\cdot \left({\mathit{\rho }}_{\mathrm{d}}^{\mathrm{r}}\mathit{u}\right)=\mathrm{0},\end{array}$$

which will be referred to as the “anelastic constraint”. Throughout the model, the pressure is assumed to be equal to the environmental pressure profile p^e(z). The only exception is the pressure gradient term appearing in Eq. (3), in which the pressure is adjusted so that u satisfies the anelastic constraint (Eq. 7) (Lipps and Hemler, 1982; Grabowski and Smolarkiewicz, 1996).

UWLCM offers two methods for modeling diffusion of Eulerian variables due to the SGS turbulence. The first is an implicit LES (ILES) approach, in which there is no explicit parametrization of SGS mixing, i.e., ${\mathcal{D}}_{*}\equiv \mathrm{0}$. Instead, numerical diffusion of the advection scheme is used to mimic the SGS turbulence (Grinstein et al., 2007). The MPDATA algorithm is argued to be well-suited for ILES simulations (Margolin and Rider, 2002; Margolin et al., 2006). The other method is a Smagorinsky-type model  (Smagorinsky, 1963; Lilly, 1962) with the SGS effects parameterized as follows:

$$\begin{array}{}\text{(8)}& {\displaystyle }{\mathcal{D}}_u& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{d}}^{\mathrm{r}}}}\mathrm{\nabla }\cdot \left({\mathit{\rho }}_{\mathrm{d}}^{\mathrm{r}}{K}_m\mathbf{E}\right),\text{(9)}& {\displaystyle }{\mathcal{D}}_{\mathit{\theta }}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{d}}^{\mathrm{r}}}}\mathrm{\nabla }\cdot \left({\mathit{\rho }}_{\mathrm{d}}^{\mathrm{r}}{K}_h\mathrm{\nabla }\mathit{\theta }\right),\text{(10)}& {\displaystyle }{\mathcal{D}}_{q_{\mathrm{v}}}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{d}}^{\mathrm{r}}}}\mathrm{\nabla }\cdot \left({\mathit{\rho }}_{\mathrm{d}}^{\mathrm{r}}{K}_q\mathrm{\nabla }{q}_{\mathrm{v}}\right).\end{array}$$

Here K_m is the eddy viscosity, K_h and K_q are the eddy diffusivities, and $\mathbf{E}=\mathrm{\nabla }\mathit{u}+(\mathrm{\nabla }\mathit{u}{)}^T-\frac{\mathrm{2}}{\mathrm{3}}(\mathrm{\nabla }\cdot \mathit{u})\mathbf{I}$ is the deformation tensor. The eddy viscosity is given by

$$\begin{array}{}\text{(11)}& K_m=\left\{\begin{array}{ll}\left(c_{\mathrm{s}}\mathit{\lambda }{)}^{\mathrm{2}}\right|\mathbf{E}|{\left(\mathrm{1}-\frac{K_h}{K_m}Ri\right)}^{\mathrm{1}/\mathrm{2}},& \text{if }\frac{K_h}{K_m}Ri<\mathrm{1}\\ \mathrm{0}& \text{otherwise},\end{array}\right.\end{array}$$

where c_s is the Smagorinsky constant, λ is the mixing length, and Ri is the Richardson number. The eddy diffusivities are

$$\begin{array}{}\text{(12)}& K_h=K_q=K_m/Pr,\end{array}$$

where Pr is the Prandtl number. Following Schmidt and Schumann (1989) the mixing length is set to $\mathit{\lambda }=min(\mathrm{\Delta },c_{\mathrm{L}}z)$. Given the highly anisotropic grid cells used in stratocumulus simulations we set Δ=Δz, as in the Colorado State University System for Atmospheric Modeling (Khairoutdinov and Randall, 2003). The values of the numerical constants are taken following Schmidt and Schumann (1989) as c_s = 0.165, Pr=0.42 and c_L=0.845.

2.2 Lagrangian particles

Liquid water is modeled with a Lagrangian, particle-based microphysics scheme from the libcloudph++ library (Arabas et al., 2015). It is an implementation of the super-droplet method (SDM) (Shima et al., 2009). The key idea is to represent all liquid particles using a small number of computational particles, called super-droplets (SDs). Each SD represents a large number of real particles. The number of real particles represented by a given SD is called the multiplicity (also known as the weighting factor), and is denoted by ξ. Other attributes of SDs are the dry radius r_d, the wet radius r, the hygroscopicity parameter κ and the position x in the model domain.

The condensational growth rate of a SD is equal to that of a single real particle. We calculate it using the Maxwell–Mason approximation (see Arabas et al., 2015):

$$\begin{array}{}\text{(13)}& r{\displaystyle \frac{\mathrm{d}r}{\mathrm{d}t}}={\displaystyle \frac{D_{\text{eff}}^{\prime }}{{\mathit{\rho }}_{\mathrm{w}}}}\left(q_{\mathrm{v}}-q_{\mathrm{vs}}{a}_{\mathrm{w}}\left(r,r_{\mathrm{d}},\mathit{\kappa }\right)\exp (A/r)\right),\end{array}$$

where

$$\begin{array}{}\text{(14)}& {\displaystyle \frac{\mathrm{1}}{D_{\text{eff}}^{\prime }}}={\left(D{\mathit{\rho }}_{\mathrm{d}}\right)}^{-\mathrm{1}}+K^{-\mathrm{1}}{q}_{\mathrm{vs}}{\displaystyle \frac{l_{\mathrm{v}}}{T}}\left({\displaystyle \frac{l_{\mathrm{v}}}{R_{\mathrm{v}}T}}-\mathrm{1}\right)\end{array}$$

and water activity is calculated using the κ-Köhler parametrization (Petters and Kreidenweis, 2007):

$$\begin{array}{}\text{(15)}& a_{\mathrm{w}}\left(r,r_{\mathrm{d}},\mathit{\kappa }\right)={\displaystyle \frac{r^{\mathrm{3}}-r_{\mathrm{d}}^{\mathrm{3}}}{r^{\mathrm{3}}-r_{\mathrm{d}}^{\mathrm{3}}\left(\mathrm{1}-\mathit{\kappa }\right)}}.\end{array}$$

Following Lipps and Hemler (1982), the relative humidity is defined as $\mathit{\varphi }=q_{\mathrm{v}}/q_{\mathrm{vs}}$ and the saturation water vapor mixing ratio is calculated using the formula $q_{\mathrm{vs}}=\left(R_{\mathrm{d}}/R_{\mathrm{v}}\right)e_{\mathrm{s}}/\left(p^{\mathrm{e}}-e_{\mathrm{s}}\right)$. Formulas for A and l_v can be found in Arabas and Pawlowska (2011). The vapor and heat diffusion coefficients D and K include gas kinetic and ventilation effects and are evaluated as in Arabas et al. (2015).

Collision–coalescence of droplets is treated as a stochastic process (Gillespie, 1972). Collisions are possible only between droplets that are located within the same spatial cell, called the coalescence cell. It is assumed that coalescence cells are well-mixed, i.e., that droplets are randomly and uniformly distributed within a coalescence cell. Then, the probability that any two real droplets j and k that are located in the same coalescence cell coalesce during the time interval Δt_c is given by the equation (Shima et al., 2009)

$$\begin{array}{}\text{(16)}& P_{j,k}=K_{j,k}{\displaystyle \frac{\mathrm{\Delta }{t}_{\mathrm{c}}}{\mathrm{\Delta }V}},\end{array}$$

where K_j,k is the collision–coalescence kernel for these two droplets and ΔV is the volume of the coalescence cell. The probability of coalescence of SDs needs to be increased to account for the fact that each SD represents a large number of real droplets. The probability that any two SDs j and k that are in the same coalescence cell coalesce during the time interval Δt_c is related to the probability of coalescence of real droplets in the following manner (Shima et al., 2009):

$$\begin{array}{}\text{(17)}& P_{j,k}^{\text{SD}}={\mathit{\xi }}_k{P}_{j,k},\end{array}$$

where SDs are labeled so that ξ_j≤ξ_k and this convention is assumed throughout the rest of this paragraph. Coalescence of the two SDs is interpreted as a coalescence of ξ_j pairs of real droplets. Each pair consists of one real droplet represented by the jth SD and one real droplets represented by the kth SD. The remaining ξ_k−ξ_j real droplets represented by the kth SD are not affected by the coalescence of these two SDs.

Such treatment of coalescence, sometimes referred to as the “all-or-nothing algorithm”, assures that coalescence does not increase the number of SDs. This algorithm was found to give the best results in a recent comparison of various coalescence algorithms used in Lagrangian schemes for microphysics (Unterstrasser et al., 2017). Dziekan and Pawlowska (2017) showed that the all-or-nothing algorithm produces correct realizations of the stochastic coalescence process described in Gillespie (1972), but only for ξ=1. For ξ>1, an average over realizations of the all-or-nothing algorithm is in good agreement with the expected value of the stochastic process, but the variability between realizations is much higher. This is because the number of SDs is much smaller than the number of real droplets. Consequently, the statistical sample for ξ>1 is much smaller than in the more realistic case of ξ=1. The collision–coalescence algorithm is not the only cause of the high variability for ξ>1. Motion of SDs is also expected to give a high variability, because when a SD moves from one spatial cell to another, a large number of real particles is abruptly moved between these cells. It is not certain if the high variability in the SDM associated with collision–coalescence of SDs and motion of SDs has any impact on averaged properties of a modeled cloud. To determine if it does have an effect, we conduct simulations for various numbers of SDs (see Sect. 4.2).

Super-droplets are treated as non-inertial particles that always sediment with their terminal velocity. There is an option to model the diffusion of liquid water due to the SGS turbulence by adding a random velocity component ${{\mathit{u}}^{\mathbf{\prime }}}_{\mathrm{SD}}$ that is specific to each SD. Each component of this velocity perturbation evolves according to Eq. (10) from Grabowski and Abade (2017). It is important to note that this SGS velocity can only be added when the Smagorinsky scheme is used for Eulerian variables. Altogether, the velocity of a SD is equal to ${\mathit{u}}_{\mathrm{SD}}=\mathit{u}+{{\mathit{u}}^{\mathbf{\prime }}}_{\mathrm{SD}}+\left(\mathrm{0},\mathrm{0},w_{\mathrm{t}}\right)+\left(\mathrm{0},\mathrm{0},w_{\text{LS}}\right)$. This formula represents the combined effects of transport by the resolved air flow, SGS turbulence, sedimentation and large-scale subsidence.

3 Numerical algorithms

3.1 Numerical integration of Eulerian equations

Numerical integration of the governing Eulerian equations is carried out using the MPDATA algorithm implemented in libmpdata++ (Jaruga et al., 2015). MPDATA is an algorithm for solving the generalized transport equation (Smolarkiewicz, 2006)

$$\begin{array}{}\text{(18)}& {\partial }_t\left(G\mathit{\psi }\right)+\mathrm{\nabla }\cdot \left(G\mathit{u}\mathit{\psi }\right)=GR,\end{array}$$

where ψ is a scalar field advected by the velocity field u, R is the source/sink right-hand side (RHS), and G can represent the fluid density, the Jacobian of coordinate transformation or their product. The equivalent of Eq. (18) in Lagrangian description is

$$\begin{array}{}\text{(19)}& D_t\mathit{\psi }=R.\end{array}$$

Eq. (3) for components of vector u and Eqs. (4) and (5) have the same form as Eq. (19). Equation (19) introduces notation that is convenient for presenting the numerical integration procedure of UWLCM. All RHS terms, except buoyancy and pressure gradient terms in Eq. (3), are integrated with the forward Euler method. These terms are denoted by R_E. The buoyancy and pressure gradient terms, denoted by R_T, are applied using the trapezoidal rule. The integration algorithm is as follows:

$$\begin{array}{}\text{(20)}& \begin{array}{rl}{\displaystyle }{\mathit{\psi }}^{\left[n+\mathrm{1}\right]}=& {\displaystyle }\mathrm{ADV}\left({\mathit{\psi }}^{\left[n\right]}+\mathrm{\Delta }t{R}_{\text{E}}^{\left[n\right]}+\mathrm{0.5}\mathrm{\Delta }t{R}_{\text{T}}^{\left[n\right]},{\mathit{u}}^{\left[n+\mathrm{1}/\mathrm{2}\right]}\right)\\ {\displaystyle }& {\displaystyle }+\mathrm{0.5}\mathrm{\Delta }t{R}_{\text{T}}^{\left[n+\mathrm{1}\right]},\end{array}\end{array}$$

where ADV(ψ,u) is an operator representing MPDATA advection of a scalar field ψ by the velocity field u. Superscripts denote the time level. The mid-time-level velocity field ${\mathit{u}}^{\left[n+\mathrm{1}/\mathrm{2}\right]}$ is obtained by linear extrapolation from u^[n−1] and u^[n].

The pressure perturbation π is adjusted so that the velocity field satisfies Eq. (7). By applying Eq. (7) to the equation for u^[n+1] discretized in the form of Eq. (20), the following elliptic equation for π^[n+1] is obtained:

$$\begin{array}{}\text{(21)}& \mathrm{\nabla }\cdot \left[{\mathit{\rho }}_{\mathrm{d}}^{\mathrm{r}}\left(\widehat{\mathit{u}}+\mathrm{0.5}\mathrm{\Delta }t\mathit{k}{B}^{\left[n+\mathrm{1}\right]}-\mathrm{0.5}\mathrm{\Delta }t\mathrm{\nabla }{\mathit{\pi }}^{[n+\mathrm{1}]}\right)\right]=\mathrm{0},\end{array}$$

where

$$\begin{array}{}\text{(22)}& \begin{array}{rl}{\displaystyle }\widehat{\mathit{u}}=& {\displaystyle }\mathrm{ADV}\left[{\mathit{u}}^{\left[n\right]}+\mathrm{\Delta }t\left({\mathit{F}}_u^{\left[n\right]}+{\mathcal{D}}_u^{\left[n\right]}\right)\right.\\ {\displaystyle }& {\displaystyle }\left.+\mathrm{0.5}\mathrm{\Delta }t\left(-\mathrm{\nabla }{\mathit{\pi }}^{\left[n\right]}+\mathit{k}{B}^{\left[n\right]}\right),{\mathit{u}}^{\left[n+\mathrm{1}/\mathrm{2}\right]}\right]\end{array}\end{array}$$

and the thermodynamic fields required in B^[n+1] are already available when the equation has to be solved. The pressure problem stated in Eq. (21) is solved with the generalized conjugate residual solver (Smolarkiewicz and Margolin, 2000; Smolarkiewicz and Szmelter, 2011).

Figure 1 UML sequence diagram showing the order of operations within a single time step. Calls in boldface start microphysical calculations that are carried out on GPUs simultaneously with solver operations undertaken on CPUs. The RHS is divided into condensational and non-condensational parts, $R=R_{\mathrm{n}}+R_{\mathrm{c}}$.

3.2 Numerical algorithms for super-droplets

For numerical reasons, condensational growth of SDs is solved in terms of the squared wet radius (Chen, 1992; Shima et al., 2009). Integration of Eq. (13) is carried out with a scheme that is implicit with respect to the wet radius and explicit with respect to q_v and θ:

$$\begin{array}{}\text{(23)}& r^{\mathrm{2}[n+\mathrm{1}]}=r^{\mathrm{2}\left[n\right]}+\mathrm{\Delta }t{\left.{\displaystyle \frac{\mathrm{d}{r}^{\mathrm{2}}}{\mathrm{d}t}}\right|}_{r^{\mathrm{2}\left[n+\mathrm{1}\right]},q_{\mathrm{v}}^{\left[n\right]},{\mathit{\theta }}^{\left[n\right]}}.\end{array}$$

Solution for Eq. (23) is found with a predictor–corrector procedure. We refer the reader to Arabas et al. (2015) for details on this procedure. Condensation can rapidly change the radii of small droplets. Therefore to correctly model condensation, in particular during the crucial moment of droplet activation, it is necessary to model condensation with a relatively short time step. Tests performed in a kinematic 2-D model of stratocumulus clouds have shown that the number of activated droplets converges for a condensation time step of around 0.1 s. A typical time step Δt of a LES model is around 1 s. Therefore it is necessary to undertake several condensation time steps in a single LES time step, a procedure we call substepping. To explain the idea of the substepping algorithm, we introduce the following notation: S_c for the number of substeps, $\mathit{\psi }=\left(\mathit{\theta },q_{\mathrm{v}}\right)$ for a vector of Eulerian variables, ψ_old for values of Eulerian variables after the substepping algorithm finished in the previous time step and ψ_new for values of Eulerian variables before the start of the substepping algorithm in the current time step. In the first substep, Eulerian variables are set to ${\mathit{\psi }}_{\text{old}}+\frac{{\mathit{\psi }}_{\text{new}}-{\mathit{\psi }}_{\text{old}}}{S_{\mathrm{c}}}$ and then condensation is calculated using the procedure defined in Eq. (23). Please note that this condensation procedure changes Eulerian variables. In each subsequent time step, $\frac{{\mathit{\psi }}_{\text{new}}-{\mathit{\psi }}_{\text{old}}}{S_{\mathrm{c}}}$ is added to Eulerian variables and then the condensation procedure is run again. Two types of the substepping algorithm are considered that differ only in the spatial cell from which the value of ψ_old is diagnosed. In the “per-particle algorithm”, ψ_old is diagnosed from the cell in which the given SD was in the previous time step. In the “per-cell algorithm”, ψ_old is diagnosed from the cell in which given SD is in the current time step. The per-cell algorithm is computationally less demanding, because ψ_old is the same for all SDs in a given spatial cell. In the per-particle algorithm ψ_old can be different for different SDs in the same cell, so each SD needs to remember its own value of ψ_old. Moreover, in the per-particle algorithm, values of pressure and density also need to vary between substeps. This is not necessary in the per-cell algorithm, because pressure and density in a given cell are constant in time. A more detailed description of the substepping algorithms and a comparison of the results they produce is given in Appendix B. The conclusion from that comparison is that the per-cell algorithm is correct for stationary clouds, but gives significant errors if cloud edge moves, whereas the per-particle algorithm is correct in both cases. All presented results of modeling stratocumulus clouds were obtained using the per-cell algorithm.

The stochastic collision–coalescence process described in Sect. 2.2 is modeled with a Monte Carlo algorithm developed by Shima et al. (2009). The key feature of this algorithm is that each SD can only collide with one other SD during a time step. Hence, the computational cost of the algorithm scales linearly, and not quadratically, with the number of SDs. Dziekan and Pawlowska (2017) showed that this “linear sampling” technique does not affect the mean (or the standard deviation) of the results. Note that in the coalescence algorithm of Shima et al. (2009), the same pair of SDs can collide multiple times during one time step. This feature was not implemented in libcloudph++ when the Arabas et al. (2015) study was published. libcloudph++ has since been modified and multiple collisions are now allowed. It is possible to run the coalescence algorithm with a shorter time step than the model time step; thus, coalescence is calculated more than once in each model time step, a procedure we call coalescence substepping.

The procedure for the initialization of SD sizes is described in detail in Dziekan and Pawlowska (2017), where it is called the “constant SD” initialization. In short, the range of initial values of r_d is divided into N_SD bins, which have the same size in log (r_d). In each bin, a single value or dry radius is randomly selected and assigned to a single SD. Multiplicity of the SD is readily calculated from the initial aerosol size spectrum. Next, the wet radius is initialized to be in equilibrium with the initial relative humidity. If the initial relative humidity is higher than 0.95, the wet radii are initialized as if it was equal to 0.95. This procedure is performed for each spatial cell. This initialization algorithm gives a good representation of the initial size spectrum even for small values of N_SD.

Advection of SDs is modeled with a predictor–corrector algorithm described in Grabowski et al. (2018a). Simpler, first-order algorithms for advection were found to cause inhomogeneous spatial distributions of SDs, with less SDs in regions of high vorticity.

3.3 Order of operations

The sequence of operations undertaken in a single time step is presented on a Unified Modeling Language (UML) sequence diagram in Fig. 1. The diagram is a convenient way of showing how coupling between Eulerian dynamics and Lagrangian microphysics is carried out. The diagram also shows operations that are undertaken simultaneously on CPUs and GPUs. Please note how the liquid water mixing ratio q_l is treated. Liquid water is resolved by the SDM and q_l could be diagnosed from the super-droplet size spectrum each time it is needed in the buoyancy term in Eq. (3) or radiative term in Eq. (4). Buoyancy is integrated with a trapezoidal scheme, which requires q_l after advection to be known. In a straightforward implementation, in which q_l is diagnosed from SDs after the advection of SDs, pressure solver calculations can only be started after advection of SDs has been calculated. Hence, there is little parallelism of calculations on GPUs and CPUs. To achieve more parallelism, we introduce an auxiliary Eulerian field for q_l. The value of q_l is diagnosed from SDs once per time step, after the condensation calculation. Then, q_l advection is carried out using a first-order accurate upwind scheme. Using the auxiliary q_l field, it is possible to calculate coalescence and motion of SDs simultaneously with calculations of advection of Eulerian fields and of the pressure problem.

3.4 Spatial discretization

Eulerian dependent variables of the model are co-located. Their positions form the nodes of the primary grid. However, the libmpdata++ advection algorithms are formulated using a dual, staggered Arakawa-C grid (Arakawa and Lamb, 1977). The cell centers of the dual grid are the nodes of the primary mesh. A schematic of a 2-D computational domain with the Arakawa-C grid is shown in Fig. 2. Throughout this paper, by “'grid cells”, “Eulerian cells” or simply “cells”, we refer to the cells of the dual grid. To form the Arakawa-C arrangement, components of the vector u are linearly interpolated to the edges of the dual grid (see Jaruga et al., 2015 for details). Super-droplets are restricted to the physical space, which is the shaded region in Fig. 2. Coupling of Eulerian variables with SDs is done using the dual grid. All SDs that are located in the same cell of the dual grid are subjected to the same conditions that are equal to the values of scalars residing at the center of the cell. Similarly, condensation of a given SD affects scalars in the center of the dual grid cell, in which this SD is located. To calculate the velocity of air that advects a given SD, velocities, which reside at the edges of the dual grid, are interpolated to the position of the SD. The interpolation is done linearly, separately in each dimension, as advocated by Grabowski et al. (2018a). Spatial discretization is also necessary in the algorithm for modeling collision–coalescence (cf. Sect. 2.2). We also use the dual grid cells as coalescence cells, with the exception of the cells at the domain edges. There, only the physical (shaded) part of dual grid cells is used as coalescence cells.

Figure 2 Schematic of a 2-D computational domain. Bullets mark the data points for the dependent variable ψ in Eq. (18), solid lines depict the edges of the primary grid and dashed lines mark the edges of the dual grid. Reproduced from Jaruga et al. (2015).

4 Comparison with other models – marine stratocumulus simulations

UWLCM simulations of a marine stratocumulus cloud are presented in this section. The main goal is to validate UWLCM by comparing it with other LES models. We study sensitivity of the model to the way the SGS turbulence is modeled and to values of the microphysics scheme parameters. Results of this sensitivity study may provide some guidance to other users of Lagrangian microphysical schemes. The simulation setup is based on observations made during the second Dynamics and Chemistry of Marine Stratocumulus (DYCOMS-II) field study (Stevens et al., 2003). The setup, described in detail in Ackerman et al. (2009), is an idealization of conditions observed during the second research flight (RF02) of this campaign. Both heavily drizzling open cells and lightly drizzling closed cells were sampled by RF02. The initial thermodynamic conditions are an average of both types of cells and the microphysical conditions are an average over heavily drizzling cells only. Comparison of simulation results from 11 different LES models is presented in Ackerman et al. (2009). There is a large variability in the amount of drizzle predicted by different models, which illustrates how difficult it is for LES models to reproduce precipitation formation. One of the reasons why we chose to test UWLCM using this setup is to test how well Lagrangian microphysics performs in modeling drizzle. The models that took part in the intercomparison use either bin microphysics (one model with single-moment bin and one model with double-moment bin) or bulk microphysics (two models with single-moment bulk and seven models with double-moment bulk).

4.1 Simulation setup

The simulation setup follows Ackerman et al. (2009). The domain size is 6.4 km × 6.4 km × 1.5 km with a regular grid of cells of 50 m × 50 m × 5 m size. Rigid and periodic boundary conditions are used at the vertical and horizontal edges of the domain, respectively. The simulations are run for 6 h. The initial profiles of q_v and θ give high values of supersaturation in the layer in which a cloud was observed. However, the simulation is initialized without any cloud water, because it is not known analytically what the initial wet radius distribution should be. The first part of the simulation, called the “spin-up” period, is dedicated to obtaining a stationary distribution of wet radii. During the spin-up, the collision–coalescence process is turned off and the supersaturation in the condensational growth equation is limited to 1 %. Please note that in Ackerman et al. (2009) this supersaturation limit is applied only to the activation and not to the condensational growth. This approach can not be used in UWLCM, because in UWLCM activation is not modeled as a separate process. The spin-up period is 1 h long, which was found to be long enough to reach a stationary concentration of cloud droplets – indicating a stationary spectrum of the wet radius. Aerosol is assumed to consist of ammonium sulfate with the initial size distribution as defined in Appendix A of Ackerman et al. (2009). Following Petters and Kreidenweis (2007), the hygroscopicity parameter for ammonium sulfate is κ=0.61. Collision efficiencies are taken from Hall (1980) for large droplets and from Davis (1972) for small droplets. Coalescence efficiency is set to one. Terminal velocities are calculated using a formula from Khvorostyanov and Curry (2002). libmpdata++ allows the user to choose from a number of MPDATA options. In the presented simulations, we use the “infinite-gauge” option iga for handling variable-signed fields combined with the non-oscillatory option fct.

4.2 The 2-D simulations: sensitivity study of SDM

The 2-D simulations are used to investigate sensitivity of the model to parameters of the SDM: the coalescence time step length Δt_coal and the initial number of SDs N_SD. Results are compared with 3-D simulations from Ackerman et al. (2009) in order to assert if 2-D simulations, which are computationally cheap, provide a reasonable representation of some of the features of 3-D simulations. However, it has to be kept in mind that the turbulence behavior in 2-D is fundamentally different from that in 3-D. Simulations are run for two model time step lengths: Δt=0.1 s and Δt=1 s. No substepping is carried out for Δt=0.1 s. Results of these simulations provide a reference for simulations with longer time steps. In simulations with Δt=1 s, 10 substeps for condensation are carried out; hence, the condensation time step is Δt_cond=0.1 s. Using a longer condensation time step results in the activation of too many aerosols (result not shown in the following figures for clarity). Two values of the coalescence time step are tested – Δt_coal=1 s (no coalescence substepping) and Δt_coal=0.1 s (10 coalescence substeps) – in combination with two values of the initial number of SDs – N_SD=40 (which is used in 3-D simulations) and N_SD=1000. As the goal of 2-D simulations is to study the microphysical model, SGS turbulence is modeled using the ILES approach. In 2-D, we observe significant variability in results of simulation runs undertaken for the same parameter values. The variability comes from two sources: one is that the initial thermodynamic conditions include a small random perturbation; the other is that initialization of SD radii and collision–coalescence of SDs are modeled with Monte Carlo algorithms. To compensate for this inherent variability, all UWLCM results of 2-D simulations shown are averages from ensembles of 10 simulations.

Figure 3 The 2-D UWLCM results. Time series of the domain averaged liquid water path, entrainment rate (equal to $\mathrm{d}{z}_{\mathrm{i}}/\mathrm{d}t+w_{\text{LS}}{z}_{\mathrm{i}}$), vertical velocity variance VAR(w) maximum, surface precipitation, concentration of cloud droplets in cloudy cells and cloud base height. UWLCM simulations were carried out for different model and coalescence time steps and different initial numbers of SDs per cell. Each colored line represents an average from 10 UWLCM simulations of a given type. Results of 3-D simulations from an ensemble of 11 models are shown for reference (Ackerman et al., 2009). The mean, the middle two quartiles, and the range of that reference ensemble are plotted using the black solid line, the dark shaded region and the light shaded region, respectively.

Time series of selected domain averaged variables are shown in Fig. 3. The only significant relationship between model parameters and results is that the amount of surface precipitation is ca. 2 times higher for Δt_coal=1 s than for Δt_coal=0.1 s. Stronger precipitation induces differences in LWP. Before the onset of precipitation, LWP is the same for all combinations of parameters. The amount of surface precipitation does not depend on the model time step or on the number of SDs. Interestingly, 2-D simulations show an abrupt increase in the entrainment rate and in the maximum of variance of w at a simulation time of around 3 h. This increase is preceded by the moment when the first precipitation reaches the surface. This suggests that the increase in the maximum of the variance of w and in the entrainment rate is caused by rain evaporation. The need for the spin-up period for microphysics is best seen on the N_c time series. Initially, due to the large initial supersaturation, cloud droplets form on all aerosol particles. Afterwards, N_c quickly decreases and after 1 h reaches a value of ca. 60 cm⁻³, which is in agreement with the 3-D reference simulations.

Vertical profiles from the 2-D simulations are shown in Fig. 4. As already observed in time series plots, precipitation flux strongly depends on Δt_coal. The precipitation flux profile reveals that precipitation flux also weakly depends on N_SD: it is slightly lower for N_SD=1000 than for N_SD=40. A similar observation was made in Dziekan and Pawlowska (2017), where the autoconversion efficiency was shown to increase with N_SD. The most striking differences between the 2-D UWLCM and the 3-D reference simulations are seen on the profiles of moments of the vertical velocity distribution. This is associated with the decreased dimensionality of our simulations. Interestingly, profiles of VAR(w) and of the third moment of w are in better agreement with observations (see Fig. 3 in Ackerman et al., 2009) in the 2-D UWLCM than in the 3-D reference simulations.

A conclusion for SDM modeling is that coalescence needs to be resolved with a time step of the order of 0.1 s, although more rigorous convergence tests should be carried out in the future. This conclusion is surprising, because coalescence tests of SDM in box models give correct results for time steps larger than 1 s (result not shown), and Shima et al. (2009) estimated that the coalescence algorithm should work well for Δt_coal of the order of 1 s. Moreover, one might expect that using large Δt_coal should give too little precipitation, as large Δt_coal can cause the mean number of collisions to be lower than expected. This is because SDM handles large Δt_coal by allowing multiple collisions between SDs and sometimes, when one of the SDs has low multiplicity, not all of these multiple collisions can be realized. However, we see that surface precipitation increases with Δt_coal. A possible explanation is that for large Δt_coal some of the SDs become extremely lucky and grow much faster than expected due to multiple collisions. Hence, even if the mean number of collisions is lower than it should be, some SDs become very large and cause the observed high surface precipitation. The second conclusion for SDM modeling that can be drawn from the sensitivity test is that N_SD of the order of 40 is sufficient to obtain correct domain averaged results. Certainly, this does not mean that this relatively low number of SDs is sufficient in all cases. For example, a larger number of SDs would probably be needed in simulations in which SDs have more attributes, e.g., when modeling aqueous chemistry. Also, we expect that observables other than domain averages, e.g., related to the spatial structure of a cloud, are more sensitive to the number of SDs. Schwenkel et al. (2018) present how cloud structure depends on the number of SDs in more detail. In general, 2-D UWLCM results do not deviate very much from the 3-D simulations from Ackerman et al. (2009). The biggest difference is in the cloud base height – cloud layer is significantly deeper in 2-D. This shows that cheap 2-D simulations can be used to coarsely study microphysical effects in stratocumulus clouds.

Figure 4 The 2-D UWLCM results. As in Fig. 3, but showing horizontally averaged profiles of liquid water potential temperature (defined in Ackerman et al., 2009), total water mixing ratio, liquid water mixing ratio, cloud fraction (defined in Appendix A), precipitation flux (defined in Appendix A), variance of vertical velocity, third moment of vertical velocity, supersaturation and concentration of droplets in cloudy cells. The vertical axis is altitude normalized by inversion height. The profiles are averaged over the 2 to 6 h period.

Figure 5 As in Fig. 3, but for 3-D UWLCM simulations. No averaging over ensembles is carried out, i.e., each line comes from a single UWLCM run.

4.3 The 3-D simulations: model validation and SGS effects

Based on conclusions of the 2-D sensitivity test, 3-D simulations are undertaken for N_SD=40, Δt=1 s, Δt_cond=0.1 s and Δt_coal=0.1 s. Three different models of SGS turbulence are tested: implicit LES, the Smagorinsky scheme and the Smagorinsky scheme with turbulent SGS motion of SDs. Contrary to the 2-D simulations, the 3-D simulations show very little variability between realizations, owing to the larger simulation domain. Therefore, averaging over an ensemble of simulations, which was necessary in the 2-D case, is not needed here and the results shown come from single simulation runs. Time series of the results are shown in Fig. 5. The biggest difference between different descriptions of the SGS turbulence is in the liquid water content. In ILES, in which there is no diffusion of liquid water, LWP increases over time and is much higher than in the reference simulations. Using the Smagorinsky scheme alone, i.e., without diffusion of liquid water, gives less liquid water than ILES. This indicates that diffusion of Eulerian variables in simulations with the Smagorinsky scheme is higher than in ILES. Still, LWP in this case is close the maximum from the reference models. Using the Smagorinsky scheme with turbulent motion of SDs, which models SGS diffusion of liquid water, further decreases LWP and results in better agreement with the reference models. Using the Smagorinsky scheme also gives the best agreement with the reference models for other variables: entrainment rate, maximum of VAR(w) and cloud base height. Diffusion of liquid water has a visible impact on LWP and therefore needs to be included in Lagrangian microphysics models. Unfortunately, this can not be done in ILES with SDM, because in this case a measure of the SGS energy dissipation is not readily available. Aside from decreasing LWP, SGS diffusion of liquid water is seen to decrease N_c, giving better agreement with the reference models right after the spin-up period. Afterwards, SGS diffusion of liquid water causes N_c to slowly decrease with time. A possible explanation is that in regions with little vertical motion, cloud droplets that diffuse out of supersaturated cells will evaporate, but aerosol particles that diffuse into supersaturated cells will not necessarily be activated, because condensational growth of larger cloud droplets already present in this region may consume all available supersaturation. Surface precipitation is very low in all 3-D UWLCM simulations. The 2-D UWLCM simulations give larger surface precipitation than the 3-D UWLCM simulations, which is attributed to a deeper cloud layer in 2-D. There is a very large spread in surface precipitation in the reference models, with some of them producing as little as 3-D UWLCM. The subject of discrepancy in surface precipitation is discussed in more detail in Sect. 4.4.

Profiles from 3-D simulations are presented in Fig. 6. Using the Smagorinsky scheme with SGS diffusion of liquid water gives the best agreement with the reference models, with the exception of the cloud fraction that is smaller than reference. However, the cloud fraction profile strongly depends on the definition of cloudy cells. Following Ackerman et al. (2009), we define cloudy cells as those with a concentration of cloud droplets greater than 20 cm⁻³. Conversely, most of the reference models use parameterized microphysics. Therefore, these models define cloudy cells as saturated cells. Using this definition, all UWLCM runs give maximum cloud fraction of ca. 95 %, which is in agreement with the reference models. Also, cloud cover, defined as the fraction of columns with LWP >20 g m⁻², is close to 100 % in all 3-D UWLCM simulations. The choice of the SGS diffusion model also affects the structure of the velocity field. Increasing SGS diffusion strength decreases the variance of w and increases the skewness of w, which shifts from negative for ILES to positive for the Smagorinsky scheme with the diffusion of liquid water.

Figure 6 As in Fig. 4, but for the 3-D UWLCM simulations. No averaging over ensembles is carried out, i.e., each line comes from a single UWLCM run.

4.4 Precipitation results

The purpose of this section is to study discrepancy between the amount of surface precipitation observed during the DYCOMS-II campaign (from ca. 0.25 to ca. 0.45 mm d⁻¹, Ackerman et al., 2009) and modeled by 3-D UWLCM (almost none). The precipitation flux in UWLCM is ca. 2 times lower than the average of reference simulations (cf. Fig. 6). To better understand this issue we make a comparison with the only models with bin microphysics that took part in the reference intercomparison: the Distributed Hydrodynamic Aerosol and Radiative Modeling Application (DHARMA, Stevens et al., 2002) and the Regional Atmospheric Modeling System (RAMS, RAMS Technical Description, 2016). DHARMA uses single-moment bin microphysics and RAMS uses double-moment bin microphysics. More details about DHARMA and RAMS simulations of the DYCOMS RF02 case can be found in Appendix B of Ackerman et al. (2009). We only compare our results with these two models, because, contrary to the bulk schemes, bin schemes explicitly resolve the size spectrum of droplets and do not rely on parameterizations of the collision–coalescence process, i.e., they are at a similar level of precision as the SDM. Bin microphysics are troubled by artificial broadening of the size spectrum of droplets due to numerical diffusion associated with advection in the physical space (Morrison et al., 2018). Such artificial broadening increases the rate of collision–coalescence, meaning that models with bin microphysics might produce too much precipitation. Lagrangian, particle-based schemes such as SDM have no numerical diffusion in the size spectrum.

Figure 7 Selected time series of domain averages (a–c) and vertical profiles (d–f) from the 3-D UWLCM and the two models with bin microphysics that took part in the Ackerman et al. (2009) intercomparison: DHARMA and RAMS. From the two DHARMA runs carried out in the intercomparison, the DHARMA_BO run is shown, because it uses a coalescence efficiency closer to unity, which is the value used in UWLCM and RAMS. Profiles are averaged and scaled as in Fig. 4. Reference results from all models discussed in Ackerman et al. (2009) are depicted as in Fig. 4.

Time series and profiles showing the amount of liquid water, surface precipitation and concentration of cloud droplets from UWLCM, DHARMA and RAMS are plotted in Fig. 7. The precipitation flux and surface precipitation are similar in UWLCM and RAMS. Both models produce almost no surface precipitation, aside from a short period at the start of the RAMS simulation, when the simulation has not yet reached a stationary state. The DHARMA model stands out in that the amount of surface precipitation it produces is higher, which is in agreement with observations. Cloud depth, LWP and N_c are similar in DHARMA and UWLCM. So why does DHARMA give a much higher precipitation flux? One possible explanation is that it is a result of the artificial broadening of the size spectra caused by numerical diffusion. Conversely, why does UWLCM give less precipitation than observed? Possibly, precipitation is affected by some physical processes that are not currently modeled by UWLCM. The list of such unresolved processes that could affect precipitation includes, but is not limited to, the following: SGS turbulence affecting condensation and coalescence of droplets, the lucky droplets effect and giant CCN initiating rain formation.

In the bin microphysics of RAMS and DHARMA, water droplets are artificially divided into haze particles and cloud droplets, and droplet activation is modeled as an instantaneous process (Stevens et al., 1996; Ackerman et al., 1995). Therefore even a short-lived maximum of supersaturation results in the activation of new droplets (Hoffmann, 2016). In contrast, in the Lagrangian microphysics of UWLCM all water droplets grow according to the same equations, which include curvature and solute terms, and droplet activation is not modeled as a separate process. Therefore, activation is not instantaneous, but happens over some time. Differences between the treatment of activation in bin and Lagrangian microphysics are apparent in the profile of N_c in Fig. 7. DHARMA and RAMS predict local maxima of N_c near the cloud base, where supersaturation is highest. In UWLCM, owing to the explicit treatment of activation, the timescale of activation is resolved and the local maximum of supersaturation near cloud base does not cause activation of new droplets. Therefore N_c in UWLCM monotonously increases near the cloud base.

5 Summary

We presented the University of Warsaw Lagrangian Cloud Model (UWLCM), a new large-eddy simulations model with Lagrangian particle-based cloud microphysics. The model is built by combining two open-source libraries, one for handling Eulerian dynamics and the other for implementing the Lagrangian microphysics scheme. Methods for coupling Lagrangian microphysics with Eulerian dynamics were presented, including spatial discretization, substepping algorithms and an algorithm for simultaneous computations of Eulerian and Lagrangian components. Simulations of a marine stratocumulus show that the model produces results that are in agreement with reference results from 11 other LES models, which proves the capability of Lagrangian microphysics to model realistic clouds. Both 2-D and 3-D simulations of the stratocumulus were performed. The 2-D simulations with UWLCM give reasonable results regarding microphysical phenomena at a fraction of the computational cost of the 3-D simulations, and were used to study sensitivity of the Lagrangian microphysics scheme. It was found that the condensational and collisional growth of droplets has to be modeled with a 0.1 s time step and that the number of computational particles does not affect domain averages, apart from small changes in the precipitation flux. The 0.1 s time step for condensation and coalescence is realized by carrying out multiple time steps for these processes at each model time step. Different approaches to modeling SGS turbulence were compared in 3-D simulations. The best agreement with other models is obtained using the Smagorinsky scheme and an algorithm for the SGS turbulent motion of computational particles. The implicit LES approach is troubled by the lack of diffusion of liquid water represented by Lagrangian computational particles. Surface precipitation modeled by UWLCM is lower than observed. This suggests that some physical phenomena not modeled by UWLCM, such as SGS turbulence affecting condensation and coalescence of droplets or giant CCN, are important for precipitation formation. In UWLCM, all particles, including humidified aerosols, evolve according to the same set of equations. Therefore, it is not necessary to include droplet activation as an additional process. Advantages of such an approach are most apparent near the cloud base, where bin schemes produce local maxima of cloud droplet concentration, whereas in UWLCM cloud droplet concentration increases monotonously.

Appendix A: List of symbols

Arabas et al., 2015

Table A1 List of symbols. As in Ackerman et al. (2009), cloudy cells are those with a concentration of cloud droplets greater than 20 cm⁻³. Cloud droplets are liquid particles with a radius in the range of 0.5 µm $<r<\mathrm{25}$ µm. The cloud fraction is the ratio of cloudy cells to the total number of cells. Precipitation flux in a cell is calculated as $\left(\sum \mathit{\xi }\frac{\mathrm{4}}{\mathrm{3}}\mathit{\pi }{r}^{\mathrm{3}}{w}_{\mathrm{t}}\right){\mathit{\rho }}_{\mathrm{w}}{l}_{\mathrm{v}}/V$, where V is volume of the grid cell and the sum is calculated over all SDs in the cell.

Appendix B: Condensation substepping algorithm

Consider condensation of SDs within cell i at time step n. The vector of thermodynamic conditions in that cell at the moment right before condensation is calculated is denoted by ${\mathit{\psi }}_i^{\left[n\right]}={\left({\mathit{\theta }}^{\left[n\right]},q_{\mathrm{v}}^{\left[n\right]}\right)}_i$. The number of time steps is denoted by S_c and substeps are indexed by ν, starting at ν=1. Super-droplets within cell i are numbered by μ. The vector of thermodynamic conditions that a given SD experiences at substep ν is denoted by ${\stackrel{\mathrm{˘}}{\mathit{\psi }}}_{\mathit{\mu }}^{\left[\mathit{\nu }\right]}$. Using this notation, the substepping algorithm is

$$\begin{array}{}\text{(B1)}& {\displaystyle }{\stackrel{\mathrm{˘}}{\mathit{\psi }}}_{\mathit{\mu }}^{\left[\mathit{\nu }+\mathrm{1}/\mathrm{2}\right]}={\stackrel{\mathrm{˘}}{\mathit{\psi }}}_{\mathit{\mu }}^{\left[\mathit{\nu }\right]}+{\displaystyle \frac{{\mathit{\psi }}_i^{\left[n\right]}-{\stackrel{\mathrm{˘}}{\mathit{\psi }}}_{\mathit{\mu }}^{\left[\mathit{\nu }=\mathrm{1}\right]}}{S_{\mathrm{c}}}},\text{(B2)}& {\displaystyle }{r}_{\mathit{\mu }}^{\mathrm{2}[\mathit{\nu }+\mathrm{1}]}=r_{\mathit{\mu }}^{\mathrm{2}\left[\mathit{\nu }\right]}+{\displaystyle \frac{\mathrm{\Delta }t}{S_{\mathrm{c}}}}{\left.{\displaystyle \frac{\mathrm{d}{r}^{\mathrm{2}}}{\mathrm{d}t}}\right|}_{r_{\mathit{\mu }}^{\mathrm{2}\left[\mathit{\nu }+\mathrm{1}\right]},{\stackrel{\mathrm{˘}}{\mathit{\psi }}}_{\mathit{\mu }}^{\left[\mathit{\nu }+\mathrm{1}/\mathrm{2}\right]}},\text{(B3)}& {\displaystyle }\begin{array}{rl}& {\stackrel{\mathrm{˘}}{\mathit{\psi }}}_{\mathit{\mu }}^{\left[\mathit{\nu }+\mathrm{1}\right]}={\stackrel{\mathrm{˘}}{\mathit{\psi }}}_{\mathit{\mu }}^{\left[\mathit{\nu }+\mathrm{1}/\mathrm{2}\right]}+\mathit{A}{\displaystyle \frac{\mathrm{4}}{\mathrm{3}}}{\displaystyle \frac{\mathit{\pi }{\mathit{\rho }}_{\mathrm{w}}V}{{\mathit{\rho }}_{\mathrm{d}}^{\mathrm{r}}}}\sum _{\mathit{\mu }=\mathrm{1}}^{\mathit{\mu }=N_i^{\left[n\right]}}\\ & {\mathit{\xi }}_{\mathit{\mu }}\left[{\left(r_{\mathit{\mu }}^{\mathrm{2}[\mathit{\nu }+\mathrm{1}]}\right)}^{\mathrm{3}/\mathrm{2}}-{\left(r_{\mathit{\mu }}^{\mathrm{2}\left[\mathit{\nu }\right]}\right)}^{\mathrm{3}/\mathrm{2}}\right],\end{array}\end{array}$$

where $r_{\mathit{\mu }}^{\mathrm{2}}$ is the square of the wet radius of the μth SD, $N_i^{\left[n\right]}$ is the number of SDs in cell i at time step n and $\mathit{A}=\left({\mathit{\theta }}^{\mathrm{e}}{l}_{\mathrm{v}}/\left(c_{\mathrm{pd}}{T}^{\mathrm{e}}\right),-\mathrm{1}\right)$. The sum in Eq. (B3) is calculated over all SDs in cell i at time step n. For details on the predictor–corrector algorithm for the calculation of the change of radius in Eq. (B2), see Eqs. (17)–(19) in Arabas et al. (2015). After the last substep, the value of ${\stackrel{\mathrm{˘}}{\mathit{\psi }}}_{\mathit{\mu }}^{\left[\mathit{\nu }=S_{\mathrm{c}}\right]}$ is the same for all SDs in the cell and the condensational RHS returned from the condensation algorithm is

$$\begin{array}{}\text{(B4)}& R_{\mathrm{c}}^{\left[n\right]}={\displaystyle \frac{{\stackrel{\mathrm{˘}}{\mathit{\psi }}}_{\mathit{\mu }=\mathrm{1}}^{\left[\mathit{\nu }=S_{\mathrm{c}}\right]}-{\mathit{\psi }}_i^{\left[n\right]}}{\mathrm{\Delta }t}}.\end{array}$$

The initial value ${\stackrel{\mathrm{˘}}{\mathit{\psi }}}_{\mathit{\mu }}^{\left[\mathit{\nu }=\mathrm{1}\right]}$ is equal to the thermodynamic conditions after condensation finished in the previous time step. Two ways of defining ${\stackrel{\mathrm{˘}}{\mathit{\psi }}}_{\mathit{\mu }}^{\left[\mathit{\nu }=\mathrm{1}\right]}$ are considered that differ regarding the spatial cell from which this initial condition is diagnosed:

$$\begin{array}{}\text{(B5)}& {\stackrel{\mathrm{˘}}{\mathit{\psi }}}_{\mathit{\mu }}^{\left[\mathit{\nu }=\mathrm{1}\right]}={\left({\mathit{\psi }}^{\left[n-\mathrm{1}\right]}+{R_{\mathrm{c}}}^{\left[n-\mathrm{1}\right]}\right)}_{i\left(n-\mathrm{1}\right)},\end{array}$$

referred to as “per-particle” substepping, and

$$\begin{array}{}\text{(B6)}& {\stackrel{\mathrm{˘}}{\mathit{\psi }}}_{\mathit{\mu }}^{\left[\mathit{\nu }=\mathrm{1}\right]}={\left({\mathit{\psi }}^{\left[n-\mathrm{1}\right]}+{R_{\mathrm{c}}}^{\left[n-\mathrm{1}\right]}\right)}_{i\left(n\right)},\end{array}$$

a procedure we call “per-cell” substepping. The notation i(n) stands for the index of the cell in which the μth SD was at time step n. The per-cell substepping is less accurate, but requires less computational time and uses less memory. The reason for this is that in the per-cell method, all SDs in a given cell have the same values of ${\stackrel{\mathrm{˘}}{\mathit{\psi }}}_{\mathit{\mu }}^{\left[\mathit{\nu }\right]}$. Moreover, values of pressure and density do not need to be substepped in the per-cell method, as they are constant in time in each cell.

We expect the difference between the two substepping algorithms to be particularly large near a moving cloud edge. We test this hypothesis by simulating cloud edge advection in an idealized 1-D setup, considering two spatial cells: one at saturation (q_l=0.029 g kg⁻¹, N_c≈52 cm⁻³) and the other subsaturated (relative humidity of 94 %). The boundaries are periodic and the only processes are diffusional growth and advection. A single time step (Δt=2 s) is performed, in which Eulerian fields and SDs are advected with the Courant number equal to one. The expected result is that the two cells exchange their contents without any condensation/evaporation taking place. The simulations are conducted for different substepping algorithms and different numbers of substeps. Deviations from the expected result are presented in Table B1.

Table B1 Errors caused by substepping in the 1-D simulation of cloud edge advection. The error is defined as $\mathit{ϵ}=\left(q_{\mathrm{l}}^{\mathrm{sim}}-q_{\mathrm{l}}^{\exp }\right)/q_{\mathrm{l}}^{\exp }$, where q_l is the liquid water mixing ratio diagnosed at the end of the simulation from the cell to which cloud droplets were advected, i.e., the cell that was initially subsaturated. The superscripts “sim” and “exp” denote the simulated and expected values, respectively.

The per-cell algorithm causes artificial evaporation of droplets, and this error increases with the number of substeps, whereas the per-particle algorithm produces correct results independently of the number of substeps. Therefore, we conclude that in simulations in which there is a lot of cloud edge advection it is necessary to use the per-particle algorithm.

However, in simulations with little cloud edge advection, the per-cell algorithm might be sufficient. We check this by using different substepping algorithms in a 2-D marine stratocumulus simulation (the DYCOMS RF02 case, cf. Sect. 4). To limit differences in dynamics between runs, we use the piggybacking approach (Grabowski, 2014): velocity fields from a dynamical “driver” simulation with Δt=0.1 s are used in two other “piggybacking” simulations with Δt=1 s. No substepping is carried out in the driver simulation and the piggybacking simulations have 10 substeps for condensation, one with the per-cell and the other with the per-particle algorithm. Time step length in the “driver” run is shorter than in the piggybacking runs, because we want to properly model condensation in the “driver” simulation without substepping, in order to obtain reference results. Obviously, due to the difference in time step length, the velocity field used in the “piggybacking” simulations is not exactly the same as in the “driver” simulations. However, averaged vertical profiles of moments of the vertical velocity are similar in the “driver” and the “piggybacker”, so we expect the concentration of cloud droplets in the “piggybacker” to be similar to that of the “driver”. We are interested in the cloud droplet concentration, because droplet activation is the process that requires short time steps for condensation. To limit variability between runs caused by the Monte Carlo scheme used to initialize SD radii, a large number of SDs is used in these simulations (N_SD=1000). Vertical profiles of concentration of cloud droplets from the “driver” and “piggybacker” runs are shown in Fig. B1. Interestingly, the per-cell algorithm is in better agreement with the “driver” simulation than the per-particle algorithm. The latter slightly underestimates the concentration of cloud droplets, by ca. 5 %. Nevertheless, we conclude that both substepping algorithms work well in simulations in which cloud edge does not move significantly.

Figure B1 Vertical profiles of the concentration of cloud droplets from the 2-D piggybacking simulations with different substepping algorithms, averaged over the 2–4 h period. The black line and shaded regions show results from Ackerman et al. (2009) (cf. Fig. 3).

Appendix C: Software implementation

UWLCM is coded in the C++ language. It relies heavily on two C++ libraries developed by the cloud modeling group at the University of Warsaw: libmpdata++ (Jaruga et al., 2015) for the Eulerian component and libcloudph++ (Arabas et al., 2015) for the Lagrangian component of the model. The structure of the libmpdata++ and libcloudph++ codes, and how they are used, in UWLCM is schematically depicted in Fig. C1.

libmpdata++ is a set of solvers for the generalized transport equations that use the MPDATA advection scheme. The solvers are organized in a hierarchy, ordered from solvers for simple flows to solvers for more complex flows. Each more complex solver inherits from the simpler solver in the hierarchy. Such a design simplifies code development, maintenance and reusability. Apart from the hierarchy of solvers, libmpdata++ contains three other independent modules: boundary conditions, concurrency handlers and output handlers, all of which are used in UWLCM.

libcloudph++ is an implementation of three microphysical models: SDM, a single-moment bulk model and a double-moment bulk model. The SDM is implemented using the Thrust library and the CUDA programming language. Therefore, the SDM can be run on multi-threaded CPUs as well as on multiple GPUs.

The UWLCM code is built on top of the libmpdata++ solvers. Separate parts of the UWLCM code handle different types of simulations. The “piggybacking” code makes it possible to run kinematic simulations, i.e., simulations with a prescribed velocity field. The “2-D/3-D” part of the code handles the dimensionality of the problem. The “forcings” code specifies external forcings, so it is the part of the code that depends on the simulation setup. The “microphysics” module is responsible for handling the choice of microphysics (only Lagrangian microphysics is available in the current UWLCM release). Owing to the code structure, different types of simulations, e.g., 2-D and 3-D simulations, different simulation setups or kinematic simulations, mostly use the same source code. The highest performance is achieved when UWLCM is run on a system with GPUs. In this case, the Eulerian component is calculated on CPUs and the Lagrangian component on GPUs. A large part of these computations is carried out simultaneously (cf. Fig. 1). The UWLCM code is open-source, under a version-control system and available from a git repository. Model output is produced in the HDF5 format, ready for plotting in ParaView. UWLCM also includes simple software for plotting time series and vertical profiles. A number of test programs was developed for UWLCM, libcloudph++ and libmpdata++. At this time UWLCM can only run in parallel on shared-memory systems. An implementation for distributed-memory systems is currently under development. The UWLCM code is inspired by the icicle kinematic model developed by Sylwester Arabas and Anna Jaruga (https://github.com/igfuw/libcloudphxx/models/kinematic_2D, last access: 26 June 2019).

Figure C1 Schematic depiction of the structure of the code of UWLCM, libmpdata++ and libcloudph++. Black arrows denote inheritance between classes.