On numerical broadening of particlesize spectra: a condensational growth study using PyMPDATA 1.0
 ^{1}Faculty of Physics, Astronomy and Applied Computer Science, Jagiellonian University, Kraków, Poland
 ^{2}Faculty of Mathematics and Computer Science, Jagiellonian University, Kraków, Poland
 ^{3}Zentrum für Datenverarbeitung, Johannes Gutenberg University Mainz, Mainz, Germany
 ^{4}Institute for Atmospheric Physics, Johannes Gutenberg University Mainz, Mainz, Germany
 ^{5}Institute of Atmospheric Physics, German Aerospace Center (DLR), Oberpfaffenhofen, Germany
 ^{6}Department of Atmospheric Sciences, University of Illinois at UrbanaChampaign, Urbana, IL, USA
 ^{1}Faculty of Physics, Astronomy and Applied Computer Science, Jagiellonian University, Kraków, Poland
 ^{2}Faculty of Mathematics and Computer Science, Jagiellonian University, Kraków, Poland
 ^{3}Zentrum für Datenverarbeitung, Johannes Gutenberg University Mainz, Mainz, Germany
 ^{4}Institute for Atmospheric Physics, Johannes Gutenberg University Mainz, Mainz, Germany
 ^{5}Institute of Atmospheric Physics, German Aerospace Center (DLR), Oberpfaffenhofen, Germany
 ^{6}Department of Atmospheric Sciences, University of Illinois at UrbanaChampaign, Urbana, IL, USA
Correspondence: Michael A. Olesik (michael.olesik@doctoral.uj.edu.pl) and Sylwester Arabas (sylwester.arabas@uj.edu.pl)
Hide author detailsCorrespondence: Michael A. Olesik (michael.olesik@doctoral.uj.edu.pl) and Sylwester Arabas (sylwester.arabas@uj.edu.pl)
This work discusses the numerical aspects of representing the condensational growth of particles in models of aerosol systems such as atmospheric clouds. It focuses on the Eulerian modelling approach, in which fixedbin discretisation is used for the probability density function describing the particlesize spectrum. Numerical diffusion is inherent to the employment of the fixedbin discretisation for solving the arising transport problem (advection equation describing size spectrum evolution). The focus of this work is on a technique for reducing the numerical diffusion in solutions based on the upwind scheme: the multidimensional positive definite advection transport algorithm (MPDATA). Several MPDATA variants are explored including infinitegauge, nonoscillatory, thirdorder terms and recursive antidiffusive correction (doublepass donor cell, DPDC) options. Methodologies for handling coordinate transformations associated with both particlesize spectrum coordinate choice and with numerical grid layout choice are expounded. Analysis of the performance of the scheme for different discretisation parameters and different settings of the algorithm is performed using (i) an analytically solvable boxmodel test case and (ii) the singlecolumn kinematic driver (“KiD”) test case in which the sizespectral advection due to condensation is solved simultaneously with the advection in the vertical spatial coordinate, and in which the supersaturation evolution is coupled with the droplet growth through water mass budget. The boxmodel problem covers sizespectral dynamics only; no spatial dimension is considered. The singlecolumn test case involves a numerical solution of a twodimensional advection problem (spectral and spatial dimensions). The discussion presented in the paper covers sizespectral, spatial and temporal convergence as well as computational cost, conservativeness and quantification of the numerical broadening of the particlesize spectrum. The boxmodel simulations demonstrate that, compared with upwind solutions, even a 10fold decrease in the spurious numerical spectral broadening can be obtained by an apt choice of the MPDATA variant (maintaining the same spatial and temporal resolution), yet at an increased computational cost. Analyses using the singlecolumn test case reveal that the width of the droplet size spectrum is affected by numerical diffusion pertinent to both spatial and spectral advection. Application of even a single corrective iteration of MPDATA robustly decreases the relative dispersion of the droplet spectrum, roughly by a factor of 2 at the levels of maximal liquid water content. Presented simulations are carried out using PyMPDATA – a new opensource Python implementation of MPDATA based on the Numba justintime compilation infrastructure.
1.1 Motivation and outline
The focus of this paper is on the problem of particlesize evolution for a population of droplets undergoing condensational growth. Representing the particlesize spectrum using a number density function, the problem can be stated as a populationbalance equation expressing conservation of the number of particles. Herein, the numerical solution of the problem using the MPDATA family of finite difference schemes originating in Smolarkiewicz (1983, 1984) is discussed. MPDATA stands for multidimensional positive definite advection transport algorithm and is a higherorder iterative extension of the forwardintime upwind scheme. Iterative application of upwind scheme, first using the physical velocity and subsequently with socalled antidiffusive velocities, results in MPDATA being characterised by small amounts of implicit diffusion but retaining the salient features of the upwind scheme: conservativeness, small phase error and sign preservation.
MPDATA features a variety of options allowing an algorithm variant that is appropriate to the problem at hand to be picked. This work highlights the importance of the MPDATA algorithm variant choice for the resultant spectral broadening of the particlesize spectrum. The term spectral broadening refers to the increase in width of the droplet spectrum during the lifetime of a cloud. The broadening may be associated with both physical mechanisms (including turbulent mixing, particle composition diversity, radiative heat transfer effects; see for example Feingold and Chuang, 2002) as well as with spurious artefacts stemming from the employed numerical solution technique, the latter being the focus of this work.
Cloud simulations with a detailed treatment of droplet microphysics face a twofold challenge in resolving the droplet spectrum width. First, it is challenging to model and numerically represent the subtleties of condensational growth which link the physicochemical properties of single particles with ambient thermodynamics through latent heat release and multiparticle competition for available vapour (e.g. Arabas and Shima, 2017; Yang et al., 2018), even more so when considering the interplay between particle growth and supersaturation fluctuations (e.g. Jeffery et al., 2007; Abade et al., 2018). Second, the discretisation strategies employed in representing the particlesize spectrum and its evolution are characterised by inherent limitations which constrain the fidelity of spectral width predictions (e.g. Arabas and Pawlowska, 2011; Morrison et al., 2018). Finally, corroboration of spectral width estimates from both theory and modelling against experimental data faces the problems of instrumental broadening inherent to the measurement techniques (e.g. Devenish et al., 2012, Sect. 3.2) and the problem of sampling volume choice (e.g. Kostinski and Jameson, 2000).
The width of the spectrum plays a key role in determining both the droplet collision probabilities (Grabowski and Wang, 2013) and the characteristics relevant for radiative transfer (Chandrakar et al., 2018). These in turn are reflected in parameterisations of cloud processes in largescale models. Taking climatetimescale simulation as an example, the representation of clouds remains the largest source of uncertainty (Schneider et al., 2017). Parameterisations used in climate models, for example of such cloud microphysical processes as cloud condensation nuclei (CCN) activation, are developed based on smallerscale simulations resolving particlesize spectrum evolution. Consequently, it is of importance to quantify the extent to which the dropletsize spectrum width is a consequence of (a) the physics of particle growth embodied in the governing equations and (b) the discretisation and the associated numerical diffusion.
The following introductory subsections start with a chronologically presented literature review of applications of MPDATA to the problem of condensational growth of particles. Section 2 focuses on a simple analytically solvable boxmodel test case with no spatial dimension considered and serves as a tutorial on MPDATA variants. It is presented to gather the information that is scattered across works focusing on more complex computational fluid dynamics applications of MPDATA. Presented simulations pertaining to the evolution of cloud droplet size spectrum in a cumulus cloud depict how enabling the discussed algorithm variants affects simulated droplet spectrum width. An analysis of the computational cost of different algorithm variants is carried out and corroborated with previously published data. While comprehensive from the point of view of the considered problem of diffusional growth (and hence limited to onedimensional homogeneous advection of positivesign fields), the presented material merely hints at the versatility of the algorithm. For a proper review of the MPDATA family of algorithms highlighting the multidimensional aspects and its multifaceted applications, we refer to Smolarkiewicz and Margolin (1998), Smolarkiewicz (2006), and Kühnlein and Smolarkiewicz (2017).
Section 3 covers the application of MPDATA for coupled sizespectral and spatial advection in a singlecolumn kinematic setup from Shipway and Hill (2012). First, the methodology to handle the spectralspatial liquid water advection problem taking into account the coupling with the vapour field is detailed. Second, the results obtained using different MPDATA variants are discussed focusing on the measures of spectral broadening.
Section 4 concludes the work with a summary of findings. Appendix A contains a convergence analysis based on results of multiple simulations using the boxmodel test case run with different temporal and spatial (sizespectral) resolutions and compared with the analytical solution.
All presented simulations are performed with PyMPDATA (Bartman et al., 2022a) – an opensource Python implementation of MPDATA built on top of the Numba justintime compilation framework. In terms of numerics, PyMPDATA follows libmpdata++ (Jaruga et al., 2015).
1.2 Background
There exist two contrasting approaches for modelling the evolution of dropletsize spectrum (see Grabowski, 2020, for a review): the Eulerian (fixedbin) and the Lagrangian (movingbin, movingsectional or particlebased). Overall, while the Lagrangian methods are the focus of active research and development (Morrison et al., 2020), the Eulerian schemes have been predominantly used in largescale modelling (Khain et al., 2015).
Following Liu et al. (1997) and Morrison et al. (2018), the earliest documented study employing the Eulerian numerics for condensational growth of a population of particles is that of Kovetz and Olund (1969). Several earlier works, starting with the seminal study of Howell (1949), utilised the Lagrangian approach. The numerical scheme proposed in Kovetz and Olund (1969, Eq. 10) resembles an upwind algorithm, being explicit in time and orienting the finitedifference stencil differently for condensation and evaporation.
An early discussion of numerical broadening of the cloud droplet spectrum can be found in Brown (1980) where the numerical scheme from Kovetz and Olund (1969) was improved in several ways, including the sampling of the drop growth rate at the bin boundaries (as is done herein). Brown (1980) also covers quantification of the error of the method by comparisons to analytical solutions.
In Tsang and Brock (1982), an Eulerian–Lagrangian scheme is considered that combines a methodofcharacteristics solution with spline interpolation onto a fixed grid. Based on comparison with the EulerianLagrangian scheme, the authors conclude that the upwind differencing is not suitable for aerosol growth calculations due to its unacceptable numerical diffusion. It is worth noting that the study includes considerations of the Kelvin effect of surface tension on the drop growth (not considered herein; see discussion of Eq. 5 below).
The first mention of an application of the MPDATA scheme for the problem of condensational growth can be found in Smolarkiewicz (1984). The problem is given as an example where the divergentflow option of the algorithm may be applicable (see Sect. 2.8 below).
In Tsang and Korgaonkar (1987), which is focused on evaporation, MPDATA is used as a predictor step followed by a corrective step using a Galerkin finiteelement solver. In two subsequent studies from the same group (Tsang and Rao, 1988, 1990), MPDATA is compared to other algorithms in terms of conservativeness and computational cost. In Tsang and Rao (1988), the basic threeiteration MPDATA was used. Intriguingly, it is noted there that “If the antidiffusion velocities are increased by some factor between 1.04 and 1.08, use of [corrective iteration] only once can reduce 50 % of the computing time […] without much sacrifice of accuracy”. In conclusion, the authors praised MPDATA for providing narrow size spectra. Tsang and Rao (1988) pointed out that MPDATA performs worse than the upwind scheme in terms of the prediction of the mean radius.
The “Aerosol Science: Theory and Practice” book of Williams and Loyalka (1991) contains a section (5.19) on MPDATA (termed “Smolarkiewicz method”) within a chapter focusing on the methods of solving the dynamic equation describing aerosol spectrum evolution. The basic variant of MPDATA (Smolarkiewicz, 1983) is presented with an outline of its derivation.
In Kostoglou and Karabelas (1995) and Dhaniyala and Wexler (1996), the authors mention that MPDATA has the potential to reduce the numerical diffusion as compared to the upwind scheme in the context of particlesize evolution calculations. However, the high computational cost of the method is given as a reason to discard the scheme from further analysis.
In Morrison et al. (2018), a comparison of different numerical schemes for the condensational growth problem is performed. Fixedbin and movingbin approaches are compared. MPDATA (referred to as MPDG therein) with the nonoscillatory option enabled is reported to produce most significant numerical diffusion and spectral broadening among employed fixedbin methods. Intriguingly, as can be seen in Fig. 7 therein, the broad spectrum appears as early as after 20 time steps, at an altitude of 20 m (out of 520 m of the simulated parcel displacement). Thus, it is questionable if the broadening attributed to the diffusivity of MPDATA is not an artefact of the tophat initial condition.
In Wei et al. (2020), MPDATA is employed for integrating droplet spectrum evolution for comparison with a Lagrangian scheme. The work concludes that the spurious broadening of the spectrum cannot be alleviated even with a spectral discretisation involving 2000 size bins.
The discussion presented in Morrison et al. (2018) prompted further analyses presented in Hernandéz Pardo et al. (2020) and Lee et al. (2021). These studies highlighted that, in principle, the problem is a fourdimensional transport problem (three spatial dimensions and the spectral dimension) and that the interplay of spectral and spatial advection further nuances the issue of spectral broadening.
It is worth noting that none of the works mentioned above discussed coordinate transformations to nonlinear grid layouts with MPDATA (a discussion of handling nonuniform mesh with the upwind scheme can be found in Li et al., 2017, their Appendix A). Wei et al. (2020) and Morrison et al. (2018) are the only works mentioning any options other than the basic flavour of the scheme, yet only the nonoscillatory option was considered. Herein, the applicability for solving the condensational growth problem of multiple variants of MPDATA and their combinations is expounded.
1.3 Governing equations
To describe the conservation of particle number N under the evolution of the particlesize spectrum ${n}_{\mathrm{p}}\left(p\right)=\frac{\mathrm{d}N}{\mathrm{d}p}$ (n denoting number density as a function of particlesize parameter p such as radius or volume), one may take the onedimensional continuity equation (i.e. Liouville equation expressing the conservation of probability; for discussion see Hulburt and Katz, 1964), in a generalised coordinate system:
where G≡G(x) represents the coordinate transformation from p to x with x being an equidistant mesh coordinate used in the numerical solution, n_{p}≡n_{p}(p(x)) being number density function and u≡u(x) denoting the pace of particle growth in the chosen coordinate x. Equation (1) in this context is also referred to as a population balance equation (see for example Ramkrishna, 2000, chap. 2.7).
The coordinate transformation term G may play a 2fold role in this context. First, one has the choice of the particlesize parameter used as the coordinate (i.e. the argument p of the density function n(p)). For the chosen coordinates $p\in [r,s\sim {r}^{\mathrm{2}},v\sim {r}^{\mathrm{3}}]$, the appropriate distributions will be n_{r}(r), n_{s}(s) and n_{v}(v), where s=4πr^{2} and $v=\mathrm{4}/\mathrm{3}\phantom{\rule{0.125em}{0ex}}\mathit{\pi}{r}^{\mathrm{3}}$ denote particle surface and volume, respectively. The size spectrum n_{p}(p) in a given coordinate is related with n_{r}(r) via the following relation of measures: n_{p}(p)dp=n_{r}(r)dr such that the total number $N=\int {n}_{\mathrm{r}}\mathrm{d}r$ is conserved. Second, one has the choice of the grid layout p(r(x)), that is, how the parameters r, s or v are discretised to form the equidistant grid in x. This can be used, for instance, to define a massdoubling grid layout (x=ln _{2}(r^{3})) as used in Morrison et al. (2018) and herein.
Combination of the two transformations yields the following definition of G:
which defines the transformation from the coordinate p of the density function to the numerical mesh coordinate x. For further discussion of the coordinate transformation approaches in the embraced framework (including multidimensional setting), see Smolarkiewicz and Clark (1986) and Smolarkiewicz and Margolin (1993).
2.1 Upwind discretisation
The numerical solution of Eq. (1) is obtained on a grid defined by $x=i\cdot \mathrm{\Delta}x$ and at discrete time steps defined by $t=n\cdot \mathrm{\Delta}t$. Henceforth, ${\mathit{\psi}}_{i}^{n}$ and G_{i} denote the discretised number density n_{p} and the discretised coordinate transformation term, respectively. The dimensionless advective field is denoted by $GC=\frac{\mathrm{d}p}{\mathrm{d}x}u\mathrm{\Delta}t/\mathrm{\Delta}x$, where C stands for the Courant number, i.e. the velocity in terms of temporal and spatial grid increments. Note that the values of the Courant number itself are not used, only the product GC (of the coordinate transformation term G and the Courant number C) is used. A staggered grid is employed and indicated with fractional indices for vector fields, e.g. $G{C}_{i+\mathrm{1}/\mathrm{2}}\equiv {\left(\right)}_{\left(GC\right)}i+\mathrm{1}/\mathrm{2}$ in the case of the discretisation of GC. A finitedifference form of the differential operators is introduced embracing the socalled upwind approach (dating back at least to Courant et al., 1952, Eq. 16 therein):
with
where the introduced flux function F defines the flux of ψ across a gridcell boundary. Hereinafter a shorthand notation ${F}_{i+\frac{\mathrm{1}}{\mathrm{2}}}\left(\mathit{\psi}\right)\equiv F({\mathit{\psi}}_{i},{\mathit{\psi}}_{i+\mathrm{1}},G{C}_{i+\frac{\mathrm{1}}{\mathrm{2}}})$ is used.
2.2 Boxmodel test case and upwind solution
The test case is based on Fig. 3 from East (1957) – one of the early papers on the topic of cloud droplet spectral broadening. The case considers the growth of a population of cloud droplets through condensation in the equilibrium supersaturation limit, where
with $\mathit{\xi}={\mathit{\xi}}_{\mathrm{0}}(S\mathrm{1})$ being an approximately constant factor proportional to the supersaturation (S−1) where the saturation S denotes the relative humidity of ambient air. The parameter ξ_{0} is set to 100 µm^{2} s^{−1} to match the results from East (1957). The approximation (Eq. 5) neglects the dependence of particle growth rate on the surface tension (Kelvin term). Taking it into consideration requires replacing (S−1) with $(S{e}^{A/r})$, where A depends on temperature only; for discussion see for example Tsang and Brock (1982).
For the initial number density function, an idealised fairweather cumulus droplet size spectrum is modelled with a lognormal distribution:
with r_{0}=7 µm and κ=22 (East and Marshall, 1954), and n_{0}=465 cm^{−3} to match the liquid water content of 1 g kg^{−1} as indicated in East (1957).
For the boundary conditions (implemented using halo grid cells), linear extrapolation is applied to G, while both ψ and GC are set to zero within the halo.
Analytical solution to Eq. (1) is readily obtainable for $\dot{r}=\mathit{\xi}/r$ and for any initial size spectrum. Noting that introducing x=r^{2} coordinates, the transport Eq. (1) becomes a constantcoefficient advection equation; the problem reduces to translation in x by 2ξt. Cast in the r coordinate, the solution can be expressed as follows (Kovetz, 1969):
where $\stackrel{\mathrm{\u0303}}{r}=\stackrel{\mathrm{\u0303}}{r}(r,t)=\sqrt{{r}^{\mathrm{2}}\mathrm{2}\mathit{\xi}t}$.
The upper panels in Figs. 1 and 2 depict the droplet size spectrum evolution through condensational growth from an initial liquid water mixing ratio of M_{0}=1 g kg^{−1} under supersaturation of $S\mathrm{1}=\mathrm{0.075}\phantom{\rule{0.125em}{0ex}}\mathit{\%}$.
Two grid layout (x) and size parameter (p) choices are depicted. Figure 1 presents a simulation carried out with a p=r^{2} coordinate and discretised on a massdoubling grid (x=ln _{2}(r^{3})). Figure 2 presents simulation results obtained with x=r and p=r. In all cases, the time step is set to $\mathrm{\Delta}t=\frac{\mathrm{1}}{\mathrm{3}}$ s. The domain span is 1–26 µm. The grid is composed of 75 grid cells. Such settings correspond to GC≈0.26 for p=r^{2}, and a variable Courant number approximately in the range of 0.03 to 0.07 for p=r.
The times for which results are depicted in the plots are selected by finding t for which the integrated liquid water mixing ratio of the analytical solution is equal to 1, 4 and 10 g kg^{−1} (assuming air density of 1 kg m^{−3}). In both Figs. 1 and 2, the upper panels show the number density and the bottom panels show the normalised mass density. The bottom panels thus depict the same quantities as Fig. 3 in East (1957). A similar comparison of upwind and analytical solutions is also presented in Fig. 1a in Li et al. (2017).
The normalised mass density of bin i is evaluated as $\mathrm{4}/\mathrm{3}\mathit{\pi}{\mathit{\rho}}_{L}{m}_{i}^{(l=\mathrm{3})}/M$, where ρ_{L}=1000 kg m^{−1}, by calculating the third statistical moment of the number distribution n_{r}(p) with the formula:
where r_{1}, r_{2} are the boundaries of ith bin, and ψ_{i} is the value of n_{p} associated with the bin (i.e. n_{p} is assumed to be binwise constant; note that the physical unit associated with n_{p} depends on the choice of p). The normalisation factor M is the water mixing ratio (e.g. $M={M}_{\mathrm{0}}=\mathrm{1}$ g kg^{−1} for t=0).
Looking at the mass density plots in Figs. 1 and 2, it is evident that casting the results in the form of mass density shifts positions of the extrema in comparison with the analytical solution. This is one of the consequences of solving the number conservation equation (for discussion see Sect. 2.12). As can be seen in both the number and massdensity plots in Figs. 1 and 2, solutions obtained with the upwind scheme are characterised by a significant drop in the peak value and visible broadening of the spectrum with respect to the analytical solution – both manifesting the numerical diffusion. The broadening and the drop in the peak value are less pronounced in Fig. 2 where employment of a linear grid causes an increase in the resolution in the largeparticle region of the spectrum as compared to massdoubling grid case of Fig. 1.
2.3 Truncation error analysis of the upwind scheme
One of the methods used to quantify the numerical diffusion of a numerical scheme is the modified equation analysis of Hirt (1968) (see Margolin and Shashkov, 2006, for discussion in the context of MPDATA). In a nutshell, the analysis involves (i) Taylor expansion of each term of the numerical scheme; (ii) elimination of higherthanfirstorder time derivatives using the timedifferentiated original advection equation; and (iii) derivation of a partial differential equation, referred to as the modified equation, that a given scheme actually approximates in lieu of the advection equation.
To depict an application of the modified equation analysis in the present context (upwind scheme), a simplified setting where G=1 and C is constant is outlined herein. In the analysis, the Taylor expansion of ψ up to the second order is taken at ${\mathit{\psi}}_{i}^{n+\mathrm{1}}$, ${\mathit{\psi}}_{i+\mathrm{1}}^{n}$ and ${\mathit{\psi}}_{i\mathrm{1}}^{n}$ and substituted into the numerical upwind scheme, in which the flux function (Eq. 4) is expressed using moduli (e.g. Crowley, 1968, Eq. 12):
which results in
which can be further transformed by employing a time derivative of both sides of the original advection equation ${\partial}_{t}\mathit{\psi}=u{\partial}_{x}\mathit{\psi}\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}\u27f6\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}{\partial}_{t}^{\mathrm{2}}\mathit{\psi}=u{\partial}_{x}{\partial}_{t}\mathit{\psi}={u}^{\mathrm{2}}{\partial}_{x}^{\mathrm{2}}\mathit{\psi}$ to substitute the secondorder time derivative with a spatial derivative (Cauchy–Kowalevski procedure; see Toro, 2009) leading to the sought modified equation (Roberts and Weiss, 1966, Eq. 2.9):
The above analysis depicts that the employment of the numerical scheme (Eq. 3) results in a solution of a modified equation (Eq. 11), approximating the original problem up to first order. The leading secondorder error contribution has the form of a diffusive term with a coefficient K (note that the above outline of the modified equation analysis assumes a constant velocity field). The diffusive form of the leading error term explains the smoothing of the spectrum evident in Figs. 1 and 2 and is consistent with the notion of numerical diffusion.
2.4 Antidiffusive velocity and iterative corrections
The problem of numerical diffusion can be addressed by introducing the socalled “antidiffusive velocity” (Smolarkiewicz, 1983). To this end, the Fickian flux can be cast in the form of an advective flux – an approach dubbed pseudovelocity technique in the context of advection–diffusion simulations (Lange, 1973, 1978) or hyperbolic formulation of diffusion (Cristiani, 2015, discussion of Eq. 4 therein) – and is discussed in detail in Smolarkiewicz and Clark (1986, Sect. 3.2):
In Smolarkiewicz (1983, 1984), it was proposed to apply the identity (Eq. 12) to Eq. (11) to suppress the spurious diffusion. The procedure is iterative. The first iteration is the basic upwind pass. Subsequent corrective iterations reverse the effect of numerical diffusion by performing upwind passes with the socalled antidiffusive flux based on Eq. (12) but with K taken with a negative sign and approximated using the upwind stencil (for discussion of the discretisation, see Smolarkiewicz and Margolin, 2001).
Accordingly, the basic antidiffusive field GC^{(k)} is defined as follows (with ϵ>0 being an arbitrarily small constant used to prevent from divisions by zero):
where k is the iteration number, GC^{(1)}≡GC and
where ψ^{*} denotes ψ^{n} in the first iteration, or the values resulting from the application of the upwind scheme with the antidiffusive flux in subsequent iterations. The MPDATA scheme inherits the key properties of the upwind scheme in terms of positivedefiniteness, conservativeness and stability while reducing the effect of numerical diffusion. In all presented formulæ below, it is assumed that ψ is positive (as in the case of particle number density). The references provided include formulation of the algorithm for variablesign fields (e.g. momentum advection).
In Fig. 3, the analytical results obtained with upwind solutions and presented in Fig. 1 are supplemented with results obtained using the MPDATA scheme with two and three iterations. Employment of MPDATA corrects (with respect to analytical solution) both the peak amplitude and the spectrum width, as well as the position of the maximum. It is visible that the effect of the third iteration is less pronounced than that of the second one. Overall, while the MPDATA solutions are superior to the upwind solutions, the drop in amplitude and broadening of the resultant spectrum still visibly differs from the discretised analytical solution.
2.5 Infinite gauge variant
For the possible improvement of the algorithm, one may consider linearising MPDATA about an arbitrarily large constant (i.e. taking ${\mathit{\psi}}^{\prime}=\mathit{\psi}+a\mathit{\chi}$ in the limit a⟶∞ instead of ψ, where χ is a constant scalar background field). Such analysis was considered in Smolarkiewicz and Clark (1986, Eq. 41) and subsequently referred to as the “infinitegauge” (or “iga”) variant of MPDATA (Smolarkiewicz, 2006, their Eq. 34; Margolin and Shashkov, 2006, point (6) on page 1204).
Such gauge transformation changes the corrective iterations of the basic algorithm as follows (replacing Eqs. 14 and 4, which is symbolised with ↝):
Note that the amplitude of the diffusive flux (Eq. 12) is inversely proportional to the amplitude of the transported field. Consequently, such a gauge choice decreases the amplitude of the truncation error (see Smolarkiewicz and Clark, 1986, p. 408, Jaruga et al., 2015, discussion of Fig. 11). However, the infinitegauge variant no longer assures positivedefinite solutions.
Figure 4 depicts how enabling the infinite gauge variant influences results presented in Fig. 3. In each plotted time step, the maximum amplitude of the infinitegauge result is closest to the analytical solution – improving over the basic MPDATA. However, in each case, negative values are also observed (nonphysical in case of the considered problem). Consequently, for the problem at hand, it is effectively essential to combine it with the monotonicitypreserving nonoscillatory option outlined in the next section.
2.6 Nonoscillatory option
In Smolarkiewicz and Grabowski (1990), an extension of the MPDATA algorithm was introduced that makes the solution monotonicitypreserving. In the case of the infinitegauge variant outlined above, it precludes the appearance of negative values in the discussed solution of dropletsize spectrum evolution. The tradeoff is that the order of the algorithm is reduced (see Appendix A).
The nonoscillatory option (later referred to as “nonosc” herein) modifies the algorithm as follows:
where
and
with
Note that in the case of the infinite gauge option being enabled, F function takes the form presented in Eq. (15) (see also Hill, 2011, Sect. 2.5).
Figure 5 juxtaposes infinitegauge solutions with the nonoscillatory option switched on or off. The effectiveness of the latter variant is apparent as spurious negative values no longer occur.
2.7 DPDC
An alternative to the iterative application of the antidiffusive velocities was introduced in Beason and Margolin (1988); Margolin and Smolarkiewicz (1998) and further discussed in Margolin and Shashkov (2006), where the contributions of multiple corrective iterations of MPDATA were analytically summed. This leads to a new twopass scheme dubbed DPDC (doublepass donor cell), featuring the following form of the antidiffusive GC field:
with A_{i+½} defined in Eq. (14). Note that only one corrective iteration is performed with the DPDC variant.
An example simulation combining the doublepass (DPDC), the nonoscillatory and infinitegauge variants is presented in Fig. 6, which depicts how the solution is improved over that in Fig. 5.
2.8 Divergentflow correction
For divergent flows (hereinafter abbreviated dfl), the modified equation analysis yields an additional correction term in the antidiffusive velocity formula (see Smolarkiewicz, 1984, Eq. 38, for uniform coordinates; Margolin and Smolarkiewicz, 1998, Eq. 30, for nonuniform coordinates; and Waruszewski et al., 2018, Sect. 4, for the infinitegauge variant):
As pointed out in Sect. 5.1 in Smolarkiewicz (1984), this option has the potential of improving results for the problem of the evolution of the droplet size spectrum (personal communication with William Hall cited therein). This is due to the drop growth velocity defined by Eq. (5) being dependent on the droplet radius and hence divergent. Yet, applying adequate coordinate transformation (i.e. p=r^{2}), the drop growth velocity in the transformed coordinates becomes constant (see Sect. 2.2 above; see for example Hall, 1980, Sect. 3b). Nevertheless, the antidiffusive velocities employed in corrective iterations of MPDATA are in principle divergent; hence the option has the potential to influence results even with p=r^{2}.
In simulations using the presented setup (also for p≠r^{2}; not shown), only insignificant changes in the results when the divergentflow option was used were observed. However, the problem considered herein does not include, for instance, the surface tension influence on the drop growth rate.
2.9 Thirdorder terms
Another possible improvement to the algorithm comes from the inclusion of the thirdorder terms in the modified equation analysis, which leads to the following form of the antidiffusive velocity (Margolin and Smolarkiewicz, 1998):
Figure 7 depicts how enabling the thirdorder terms improves the solution with respect to the upwind and basic MPDATA solutions.
It is worth noting that discussion of higherorder variants of MPDATA was carried forward in Kuo et al. (1999) and Waruszewski et al. (2018). In the latter case, the focus was placed on accounting for coordinate transformation and variable velocity in the derivation of antidiffusive velocities leading to a fully thirdorder accurate scheme.
2.10 A “best” combination of options
The MPDATA variants presented in the preceding sections can be combined together. In Fig. 8, results obtained with the upwind scheme and with the basic twopass MPDATA are compared to those obtained with three iterations and the thirdorder terms, the infinitegauge and the nonoscillatory options enabled simultaneously. This combination of options is hereinafter referred to as the “best” variant (for the problem at hand).
In the following subsections, the influence of MPDATA algorithm variant choice on the resultant spectrum width and on the computational cost is analysed using the example simulation setup used above (i.e. in all figures except Fig. 2; see Sect. 2.2 for test case definition). Analysis of the scheme solution convergence with changing resolution and Courant number is presented in Appendix A.
2.11 Quantification of numerical broadening
The relative dispersion, defined as the ratio of standard deviation σ to the mean μ of the spectrum, is a parameter commonly used to describe the width of the spectrum (e.g. Chandrakar et al., 2018).
The calculated dispersion ratio over all bins takes the following form:
where m_{i} is defined in Eq. (8) and N is the conserved total number of particles (equal to ${\sum}_{i}{m}_{i}^{(l=\mathrm{0})}$). To quantify the effect of the numerical diffusion on the broadening of the spectrum, the following parameter is introduced based on the numerical and analytical solutions (hereinafter reported in percentages):
Table 1 depicts the gradual narrowing of the spectrum under undisturbed adiabatic growth. The lefthand panel (a) in Fig. 9 provides values of the R_{d} parameter evaluated at six selected times corresponding to $M=\mathrm{1},\mathrm{2},\mathrm{4},\mathrm{6},\mathrm{8},\mathrm{10}$ g kg^{−1}. Although numerical broadening is inherent to all employed schemes and grows in time for all considered variants, the scale of the effect is significantly reduced when using MPDATA. In particular, a 10fold decrease in numerical broadening as quantified using R_{d} is observed comparing upwind and the “best” variant considered herein.
2.12 Notes on conservativeness
Due to the formulation of the problem as number conservation, and discretisation of the evolution equation using fixed bins, even though the numerical scheme is conservative (up to subtle limitations outlined below), evaluation of other statistical moments of the evolved spectrum from the number density introduces an inherent discrepancy from the analytical results (for a discussion on multimoment formulation of the problem, see Liu et al., 1997).
In order to quantify the discrepancy in the total mass between the discretised analytical solution and the numerically integrated spectrum, the following ratio is defined using the moment evaluation formula (Eq. 8):
The righthand panel (b) in Fig. 9 depicts the values of the abovedefined ratio computed for spectra obtained with different variants of MPDATA discussed herein. The departures from analytically derived values are largest for the upwind scheme (up to ca. 5 %) and oscillate around 0 with an amplitude of the order of 1 % for most of the MPDATA solutions.
The consequences of mass conservation inaccuracies in the fixedbin particlesize spectrum representation may not be as severe as in, for example, a dynamical core responsible for the transport of conserved scalar fields. The outlined discrepancies may be dealt with by calculating the change in mass during a time step from condensation, then using it in vapour and latent heat budget calculations so the total mass and energy in the modelled system are conserved.
The embraced algorithm (Eqs. 3–4) is conservative (up to numerical precision) for G=1. However, the formulation of the donor cell scheme ${\mathit{\psi}}^{n+\mathrm{1}}={\mathit{\psi}}^{n}+{G}_{i}^{\mathrm{1}}\left({F}_{i\mathrm{1}/\mathrm{2}}+{F}_{i+\frac{\mathrm{1}}{\mathrm{2}}}\right)$ on the staggered grid with G≠1, for example due to employment of nonidentity coordinate transformations, implies that even though the influx and outflux across boundary of adjacent cells is equal, discretisation of G_{i} at cell centres limits the level of accuracy in number conservation. For further discussion, see Sect. 3 in Smolarkiewicz and Rasch (1991).
The total number of particles in the system may diverge from the analytical expected value even for the initial condition depending on the employed discretisation approach. In the present work, the probability density function is sampled at cell centres effectively assuming piecewiseconstant number density function. An alternative approach is to discretise the initial probabilities by assigning to ψ_{i} the values of $({\mathit{\varphi}}_{i+\frac{\mathrm{1}}{\mathrm{2}}}{\mathit{\varphi}}_{i\mathrm{1}/\mathrm{2}})/({r}_{i+\frac{\mathrm{1}}{\mathrm{2}}}{r}_{i\mathrm{1}/\mathrm{2}})$ where ϕ is the cumulative distribution.
2.13 Computational cost
Table 2 includes an assessment of the relative computational cost of the explored variants of MPDATA. The performance was estimated by repeated measurements of the wall time and selecting the minimal value as representative. Values are reported as multiplicities of the upwind execution time. Simulations were performed using the massdoubling grid.
The table includes, where available, analogous figures reported in earlier studies on MPDATA (see caption for comments on the dimensionality of the employed cases, as it differs and thus does not warrant direct comparison). Among notable traits is the decrease in computational cost when enabling the infinite gauge option that is associated with a reduced number of terms in both the flux function as well as in the antidiffusive velocity formulation (see Sect. 2.5 in Hill, 2011, and Sect. 2.5–2.6 herein). The “best” variant is roughly 10 times more costly than the upwind scheme for the setup studied herein and the employed implementation. Among studies of bin microphysics schemes, analogous measures were reported in Liu et al. (1997) where the variational method presented there was reported to take 3.1 times longer to execute than the firstorder upwind solution, and in Onishi et al. (2010) where the studied semiLagrangian scheme was reported to be characterised by a computational cost over 4 times higher than the upwind scheme (see Table 4 therein). In the latter case, a direct comparison is hindered by significantly different stability constraints on the time step.
Although the discussed problem is onedimensional, a computationally efficient and an accurate solution is essential, as it typically needs to be solved at every time step and grid point of a threedimensional cloud model. While the reported upwindnormalised wall times give a rough estimation of the cost increase associated with a particular MPDATA option, the actual footprint on a complex simulation system will depend on numerous implementation details including parallelisation strategy.
3.1 Problem statement
In multidimensional simulations in which the considered particle number density field is not only a function of time and particle size, but also of spatial coordinates, there are several additional points to consider applying MPDATA to the problem.
First, in the context of atmospheric cloud simulations, owing to the stratification of the atmosphere, the usual practice is to reformulate the conservation problem in terms of specific number concentration being defined as the number of particles n_{p} (cf. Eq. 1) divided by the mass of air (commonly the dry air) effectively resulting in including the (dry) air density in the G factor (see Eqs. 2–3). This is motivated by atmospheric stratification associated with the presence of a vertical air density gradient. In a nondivergent stratified flow, the specific number concentration (summed across all particlesize categories) is not modified by advection along the vertical dimension. On the other hand, particle volume concentration (as opposed to specific number concentration) would vary due to variable density of air (i.e. expansion of air along the vertical coordinate). Note that in Eq. (3) it is further assumed that the G factor does not vary in time.
Second, even with a single spatial dimension (singlecolumn setup), the coupled sizespectral–spatial advection problem is twodimensional. This is where the inherent multidimensionality of MPDATA (the “M” in MPDATA) requires further attention. The onedimensional antidiffusive formulæ discussed in Sect. 2.4–2.9 need to be augmented with additional terms representing crossdimensional contributions to the numerical diffusion. For an introduction, see for example Sect. 2.2 in Smolarkiewicz and Margolin (1998), for original derivation see Smolarkiewicz (1984), and for a recent work discussing the interpretation of all terms in the antidiffusive velocity formulæ, including crossdimensional terms, see Waruszewski et al. (2018).
Third, in any practical application where the drop size evolution is coupled with the water vapour budget (and hence with supersaturation evolution), it is essential to evaluate the total change in mass of liquid water due to condensation which is then to be used to define the source term of the water vapour field (and in latent heat budget representation). It is worth noting that knowing the difference of values at n+1 and at n time steps of the advected specific number concentration field is not sufficient to evaluate the vapour sink–source term. This is because the fluxes across the sizespectral dimension only need to be accounted for (note that the fluxes in all MPDATA iterations need to be summed up).
Several recent papers are highlighting the need for scrutiny when comes to the interplay of sizespectral and spatial advection and the associated numerical broadening (Morrison et al., 2018; Hernandéz Pardo et al., 2020; Lee et al., 2021). In the following subsection, a set of singlecolumn simulations is presented and discussed depicting the performance of MPDATA in a sizespectral–spatial advection problem coupled with vapour advection and supersaturation budget. The simulations are performed using a commonly employed MPDATA setting with only the nonoscillatory option enabled, and the discussion is focused on the sensitivity of the results to spatial, spectral and temporal resolution, as well as to the effect of performing one or two corrective passes of MPDATA (two or three iterations, respectively).
3.2 Test case definition
The test setup is based on the singlecolumn KiD warm case introduced in Shipway and Hill (2012). This prescribedflow framework has been further used, e.g. in Field et al. (2012) (mixedphase scenario), and in Gettelman and Morrison (2015) (both pureice, mixedphase and warmrain scenarios). Here, condensation is the only microphysical process considered.
The simulated 3.2 km high column of air is described by the following:

a constantintime piecewiselinear potential temperature profile (297.9 K from the ground to the level of 740 m, linearly decreasing down to 312.66 K at 3260 m);

constantintime hydrostatic pressure and density profiles computed assuming surface pressure of 1007 hPa;

a piecewise linear initial vapour mixing ratio profile (15 g kg^{−1} at ground, 13.8 g kg^{−1} at 740 m and 2.4 g kg^{−1} at 3260 m); and

a constantinspace but timedependent vertical momentum profile defining the vertical component of the advector field [GC_{r},GC_{z}] as in Eq. (28):
$$\begin{array}{}\text{(28)}& G{C}_{z}(z,t)={\mathit{\rho}}_{\mathrm{d}}{u}_{z}{\displaystyle \frac{\mathrm{\Delta}t}{\mathrm{\Delta}z}}\mathrm{sin}(\mathit{\pi}t/{t}_{\mathrm{1}})(\mathrm{1}H(t{t}_{\mathrm{1}}\left)\right),\end{array}$$where H is the Heaviside step function, w is the vertical velocity, ρ_{d}u_{z}=3 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{\mathrm{1}}\phantom{\rule{0.125em}{0ex}}\mathrm{kg}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{\mathrm{3}}$, ρ_{d}(z) is the hydrostatic dry density profile and t_{1}=600 s.
Note that the vertical velocity thus differs from the original KiD setup where w is held constant across the column (rather than constant momentum density as done herein). This change is motivated by the aim of maintaining the nondivergent flow field condition.
The advection is thus solved for two scalar fields: (i) a onedimensional water vapour mixing ratio field representing the vertical distribution of mass of vapour per mass of dry air and (ii) a twodimensional field representing vertical and spectral variability of liquid particle specific concentration (number of particles per mass of dry air). The spectral coordinate is set to particle radius (p=r) and the bins are laid out uniformly (x=r) over a range of 1 to 20.2 µm. It is worth noting that this results in the sizespectral component of the advection velocity being divergent (while the vertical component is nondivergent).
The initial condition does not feature supersaturation anywhere in the domain. The upward advection of water vapour causes supersaturation to occur and trigger condensation. The sizespectral velocity is defined as in the boxmodel test case (cf. Eq. 5) but with supersaturation being timedependent and derived from the values of vapour mixing ratio, temperature and pressure at a given level. Note that the temperature profile is constant in time and the test case does not feature representation of latent heat release effects, only the ambient air and particle vapour budget is accounted for by subtracting the amount of condensed water from the vapour field in each time step, before performing the subsequent step of advection on the vapour mixing ratio field.
The domain is initially void of liquid water and the only source of it is through the boundary condition in the spectral dimension specified as follows:
with $i=\mathrm{1}$ denoting the halo grid cell at the left edge of the spectral domain on a given vertical level (the summation spans all bins at a given level excluding the halo grid cells). The flux across the domain boundary in the spectral dimension represents cloud droplet activation. Through Eq. (5), the flux is dependent on the supersaturation at a given level, and on the N_{CCN} parameter representing a maximal number of activated droplets (per unit mass of dry air). In the performed simulations, N_{CCN} was set to 500 mg^{−1}. For discussion of other ways to represent activation in bin microphysics models, see for example Grabowski et al. (2011).
The simulations cover a time period of 15 min out of which the first 10 min (until t_{1}=600 s) involve nonzero vertical velocity (cf. Eq. 28). Several temporal, spatial and spectral resolutions are tested with the following settings hereinafter referred to as the base resolution case: 32×32 grid with a vertical grid step Δz=100 m, sizespectral grid step Δr=0.6 µm and time step Δt=0.25 s.
3.3 Discussion of results
Figure 10 depicts qualitatively how MPDATA performs with the singlecolumn simulation (base resolution case) depending on the number of MPDATA iterations employed. The twodimensional liquid water mixing ratio grid is rendered with shaded histogram bars. The vertical axis corresponds to the advected quantity: spatiospectral number density divided by the dry density of air. Histogram bars with values of less than 5 % of the vertical axis range (5 %×2 m^{−1} mg^{−1} µm^{−1}) are not plotted for clarity. Presented plots are aimed at intuitively portraying the model state and the extent to which the introduction of subsequent MPDATA corrective iterations reduces spectrum broadening. Note that besides the depicted liquid water mixing ratio, the model state consists as well of a onedimensional vapour mixing ratio vector (not shown).
In Fig. 11, the base resolution case is depicted with plots constructed following the original methodology from Shipway and Hill (2012) (as in Fig. 1 therein). The greyscale maps depict the evolution in time and vertical dimension of water vapour mixing ratio q_{l}, supersaturation S and the droplet spectrum relative dispersion d. The adjacent profile plots depict the vertical variability of the mapped quantity at four selected times.
Notwithstanding the highly idealised and simplified modelling framework employed herein, one may attempt a comparison with profiles obtained from both in situ aircraft measurements (Arabas et al., 2009, profiles of d in Fig. 1 therein) and detailed threedimensional simulations (Arabas and Shima, 2013, profiles of S and liquid water content in Figs. 2–4 therein) inspired by the same RICO field campaign (Rauber et al., 2007) as the singlecolumn setup of Shipway and Hill (2012). The comparison confirms that the chosen test case covers the parameter space relevant to the studied problem. Resemblance remains, at most, qualitative, as expected given the stark simplicity of the KiD framework.
Interestingly, the parabolic vertical profile of the relative dispersion obtained herein was also reported in Lu and Seinfeld (2006) for binmicrophysics threedimensional simulations of marine stratocumulus. In the discussion of Figs. 2, 3 and 6 therein, it was hypothesised that the parabolic shape is a signature of entrainment as well as updraft–downdraft interactions, none of which are represented in the kinematic framework employed herein.
The liquid water profiles depicted in the top row of Fig. 11 reveal that the cloud structure developed within the first ca. 9 min of the simulation is later maintained, with the profiles at t=9 min and t=12 min being virtually indistinguishable. Middle row plots of supersaturation profiles depict that the considered simulation setup enables the user to capture the characteristic supersaturation maximum just above cloud base. Furthermore, it is evident that the corrective iterations of MPDATA influence the maximal supersaturation values. It is worth noting that this results in different time step (Courant number) constraints depending on the number of iterations used because the spectral velocity is a function of supersaturation.
There is a cloudtop activation feature hinted in all three panels in Fig. 10 as well as indirectly in the supersaturation profiles in Fig. 11. The representation of activation above cloud base is sensitive to both numerical details of vapour and heat transport reflected in the diagnosed supersaturation, as well as to the assumptions behind the activation formulation itself (see for example the discussion of Figs. 2 and 6 in Slawinska et al., 2012, and references therein). Given the simplified treatment of activation defined by Eq. (29), together with the unphysical assumption of constant temperature profile, no physical interpretation of this feature is warranted. Yet, it is worth noting that, consistent with the differences in supersaturation values between upwind and MPDATA solutions, the cloudtop activation is least noticeable in Fig. 10 in the case of the upwind solution.
The bottom row in Fig. 11 depicts the relative dispersion defined and computed as in Sect. 2.11 (discarding levels where the total droplet number mixing ratio summed over all bins on a level is below 5 % of N_{CCN}). Narrowing of the spectrum with height below z=1.5 km revealed by decreasing values of d is a robust feature. Minimum values of d for a given profile vary visibly depending on the number of MPDATA iterations employed.
To provide insight into the sensitivity of the results to temporal, spatial and spectral resolution, Fig. 12 presents the relative dispersion profiles at $t={t}_{\mathrm{1}}=\mathrm{10}$ min for several resolution settings. In the background of the figure, there are three axes plotted pointing in the directions in which the figure panels can be explored to reveal the dependence on the vertical spatial spacing Δz (lefttoright), the spectral spacing Δr (bottomtotop) and the time step (backtoforeground). The base resolution case is plotted at the intersection of the axes. Note that besides the backtoforeground sequence of plots where all but the time step settings is kept equal, the time step also varies with the grid settings to fulfil scheme stability constraint.
The dependence on the temporal resolution, as gauged by comparing the base resolution case with cases in which the time step is halved (Δt=125 ms; background) and is doubled (Δt=500 s; foreground), remains barely observable. This is in general agreement with Morrison et al. (2018) and Hernandéz Pardo et al. (2020) where the dependence on the time step is shown to be much smaller than on the spatial or on the spectral resolution.
The dependence on the spectral resolution is clearly manifested at the lowest spectral resolution where the minimum spectral dispersion d along a profile drops by ca. 0.1 when decreasing Δr=1.2 µm down to Δr=0.3 µm. Little further change can be observed by refining the resolution down to Δr=0.15 µm. Focusing on the minimum values of d for a given profile, in general, the lower the spectral resolution, the more profound the effect of introducing corrective iterations of MPDATA. In most cases, applying even a single corrective step (i.e. two iterations) results in halving of the minimal values of d as compared to the upwind solution (i.e. one iteration).
The spatial resolution setting Δz significantly alters the results, particularly near the cloud base. The values of d at the lower half of the presented profile (i.e. ca. below z=1 km) drop from over 0.3 down to around 0.15 when refining the resolution from Δz=200 m down to Δz=25 m.
The study focuses on the MPDATA family of numerical schemes and its application to the sizespectral as well as spatiospectral transport problem arising in models of condensational growth of cloud droplets. MPDATA iteratively applies the upwind algorithm, first with the physical velocity and subsequently using antidiffusive velocities. As a result, the algorithm is characterised by reduced numerical diffusion compared with upwind solutions, while maintaining conservativeness and positivedefiniteness.
In literature, the derivations of different MPDATA variants are spread across numerous research papers published across almost four decades, and in most cases focused on multidimensional hydrodynamics applications. It is the aim of this study to highlight the developments that followed the original formulation of the algorithm, and to highlight their applicability to the problem of bin microphysics. To this end, it was shown that the combination of such features of MPDATA as the infinitegauge, nonoscillatory and thirdorderterm options, together with the application of multiple corrective iterations offer a robust scheme that grossly outperforms the almost quadragenarian basic MPDATA.
In the case of the singlecolumn test case, discussed simulations feature coupling between droplet growth and supersaturation evolution. The embraced measure of spectrum width, the cloud droplet spectrum relative dispersion, is influenced by numerical diffusion pertinent to both spectral and vertical advection. Focusing on the levels corresponding to the region of maximal liquid water content (ca. between z=1 and 2 km for the case considered), application of even a single corrective iteration of MPDATA robustly reduces (in most cases more than halves) the spectral width. In agreement with conclusions drawn from singlecolumn simulations in Morrison et al. (2018) and Lee et al. (2021), within the range of explored grid settings, the vertical resolution has the most profound effect on the overall characteristics of the spectrum width profile as it significantly influences the justabovecloudbase evolution of the spectral width.
To assess the spatial and temporal convergence of the numerical solutions presented above, a convergence test originating from Smolarkiewicz and Grabowski (1990) is used. For the analysis the following truncationerror L^{2} measure is used (e.g. Smolarkiewicz, 1984):
As a side note, it is worth pointing out that for the chosen coordinates $\left(p={r}^{\mathrm{2}},x={r}^{\mathrm{2}}\right)$, the coordinate transformation term is equal to the identity, so there is no need for including the G factor into the computed error measures. In the general case, convergence will depend on the grid choice and to account for that one may use a modified measure as given in Smolarkiewicz and Rasch (1991, Eq. 24).
To explore the convergence, the error measures are computed for seven different linearly spaced values of C between 0.05 and 0.95, and $nx\in \left\{{\mathrm{2}}^{\mathrm{7}},{\mathrm{2}}^{\mathrm{8}},{\mathrm{2}}^{\mathrm{9}},{\mathrm{2}}^{\mathrm{10}},{\mathrm{2}}^{\mathrm{11}},{\mathrm{2}}^{\mathrm{12}},{\mathrm{2}}^{\mathrm{13}},{\mathrm{2}}^{\mathrm{14}}\right\}$ resulting in 56 simulations for each presented combination of options.
As proposed in Smolarkiewicz and Grabowski (1990), visualisation of the results is carried out on polar plots with radius ρ and angle ϕ coordinates defined as follows:
where ρ was shifted by a constant so that the highest resolution grid corresponds to ρ=1.
Figures A1–A8 depict the convergence rates and are intended for comparison with analogously constructed plots in Smolarkiewicz and Grabowski (1990, Figs. 2–3), Margolin and Smolarkiewicz (1998, Figs. 8.1–8.2) and Jaruga et al. (2015, Figs. 10–11).
The chosen colour increments correspond to the error reduction by a factor of 2, and the warmer the colour, the larger the error. The small grey points behind the isolines represent points for which the error value was calculated. When moving along the lines of constant Courant number, increasing the space and time discretisation, the number of crossed dashed isolines indicate the order of convergence. For the considered problem, it can be seen that the upwind scheme (Fig. A1) has a convergence of the first order (one isoline is crossed when spatial discretisation increases by one order); the MPDATA scheme (Fig. A2) is of the second order, and MPDATA with 3 iterations (Fig. A6) is of the third order.
Moreover, the shape of the dashed isolines tells the dependency of the solution accuracy on the Courant number. When these are isotropic (truncation error being independent of polar angle), the solution is independent of the Courant number.
It is worth noting that in Figs. A3 and A4, a groove of the thirdorder convergence rate is evident around $\mathit{\varphi}=\frac{\mathit{\pi}}{\mathrm{4}}$, normally characteristic for MPDATA with three or more passes. When secondorder truncation error is sufficiently reduced, the thirdorder error, proportional to $(\mathrm{1}\mathrm{3}C+\mathrm{2}{C}^{\mathrm{2}})$ as can be seen in Eq. (23), dominates but vanishes for C=0.5, thus resulting in the existence of the groove.
The convergence test results for the threepass MPDATA with infinite gauge, nonoscillatory and thirdorder terms options enabled (Fig. A8) are consistent with results depicted in Fig. A7, although the order of convergence is reduced due to the employment of nonoscillatory option.
Calculations presented in the paper were performed using Python with a new opensource implementation of MPDATA: PyMPDATA (Bartman et al., 2022a). All of presented figures and tables can be recreated using Jupyter notebooks included in the PyMPDATAexamples package. Both PyMPDATA and PyMPDATAexamples are licensed under the GNU General Public License 3.0 and are available on the PyPI.org Python package repository. Releases are being additionally archived on zenodo.org. The DOI links for versions 1.0 used in this study are https://doi.org/10.5281/zenodo.6329303 (PyMPDATA 1.0, Bartman et al., 2022b) and https://doi.org/10.5281/zenodo.6471494 (PyMPDATAexamples 1.0.1) (Bartman et al., 2022c).
No data sets were used in this article.
The idea of the study originated in discussions between SA, SU and MAO. MAO led the work, and a preliminary version of a significant part of the presented material constituted his MSc thesis prepared under the mentorship of SA. PB architected the key components of the PyMPDATA package. JB contributed the DPDC variant of MPDATA to PyMPDATA. MB participated in composing the paper and devising the result analysis workflow. All authors contributed to the final form of the text.
At least one of the (co)authors is a member of the editorial board of Geoscientific Model Development. The peerreview process was guided by an independent editor, and the authors also have no other competing interests to declare.
Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Comments from Wojciech Grabowski, Adrian Hill, Hugh Morrison, Andrzej Odrzywołek, Piotr Smolarkiewicz and Maciej Waruszewski as well as paper reviews by Josef Schröttle and three anonymous reviewers helped to improve the article.
The project was carried out within the POWROTY/REINTEGRATION programme of the Foundation for Polish Science (https://www.fnp.org.pl/, last access: April 2022) cofinanced by the European Union under the European Regional Development Fund (POIR.04.04.00005E1C/1800). Manuel Baumgartner acknowledges support by the Deutsche Forschungsgemeinschaft (DFG) within the Transregional Collaborative Research Centre TRR165 Waves to Weather (https://www.wavestoweather.de/, last access: April 2022) project Z2. Simon Unterstrasser acknowledges support through the DLRinternal research group H2CONTRAIL. The publication was partially funded through the “Excellence Initiative – Research University” programme at the Jagiellonian University.
This paper was edited by Juan Antonio Añel and reviewed by Josef Schröttle and three anonymous referees.
Abade, G., Grabowski, W. W., and Pawlowska, H.: Broadening of Cloud Droplet Spectra through Eddy Hopping: Turbulent Entraining Parcel Simulations, J. Atmos. Sci., 75, 3365–3379, https://doi.org/10.1175/JASD180078.1, 2018. a
Arabas, S. and Pawlowska, H.: Adaptive method of lines for multicomponent aerosol condensational growth and CCN activation, Geosci. Model Dev., 4, 15–31, https://doi.org/10.5194/gmd4152011, 2011. a
Arabas, S. and Shima, S.I.: LargeEddy Simulations of Trade Wind Cumuli Using ParticleBased Microphysics with Monte Carlo Coalescence, J. Atmos. Sci., 70, 2768–2777, https://doi.org/10.1175/JASD120295.1, 2013. a
Arabas, S. and Shima, S.: On the CCN (de)activation nonlinearities, Nonlin. Processes Geophys., 24, 535–542, https://doi.org/10.5194/npg245352017, 2017. a
Arabas, S., Pawlowska, H., and Grabowski, W.: Effective radius and droplet spectral width from insitu aircraft observations in tradewind cumuli during RICO, Geophys. Res. Lett., 36, L11803, https://doi.org/10.1029/2009GL038257, 2009. a
Bartman, P., Banaśkiewicz, J., Drenda, S., Manna, M., Olesik, M., Rozwoda, P., Sadowski, M., and Arabas, S.: PyMPDATA v1: Numbaaccelerated implementation of MPDATA with examples in Python, Julia and Matlab, J. Open Source Soft., https://joss.theoj.org/papers/10e7361e43785dbb1b3d659c5b01757a (last access: April 2022), in review, 2022a. a, b
Bartman, P., Banaśkiewicz, J., Drenda, S., Manna, M., Olesik, M. A., Rozwoda, P., Sadowski, M., and Arabas, S.: PyMPDATA 1.0, Zenodo [code], https://doi.org/10.5281/zenodo.6329303, 2022b. a
Bartman, P., Banaśkiewicz, J., Drenda, S., Manna, M., Olesik, M. A., Rozwoda, P., Sadowski, M., and Arabas, S.: PyMPDATAexamples 1.0.1, Zenodo [code], https://doi.org/10.5281/zenodo.6471494, 2022c. a
Beason, C. W. and Margolin, L. G.: DPDC (doublepass donor cell): A secondorder monotone scheme for advection, in: Fifth Nuclear Code Developers' Conference, Boulder, CO, 11–14 October 1988, LLNL report UCRL99731, https://www.osti.gov/servlets/purl/7049237 (last access: April 2022), 1988. a
Brown, R.: A numerical study of radiation fog with an explicit formulation of the microphysics, Q. J. Roy. Meteor. Soc., 106, 781–802, https://doi.org/10.1002/qj.49710645010, 1980. a, b
Chandrakar, K. K., Cantrell, W., and Shaw, R. A.: Influence of Turbulent Fluctuations on Cloud Droplet Size Dispersion and Aerosol Indirect Effects, J. Atmos. Sci., 75, 3191–3209, https://doi.org/10.1175/JASD180006.1, 2018. a, b
Courant, R., Isaacson, E., and Rees, M.: On the Solution of Nonlinear Hyperbolic Differential Equations by Finite Differences, Comm. Pure Appl. Math., 5, 243–255, https://doi.org/10.1002/cpa.3160050303, 1952. a
Cristiani, E.: Blending Brownian motion and heat equation, J. Coupled Syst. Multiscale Dyn., 3, 351–356, https://doi.org/10.1166/jcsmd.2015.1089, 2015. a
Crowley, W. P.: Numerical advection experiments, Mon. Weather Rev., 96, 1–11, https://doi.org/10.1175/15200493(1968)096<0001:NAE>2.0.CO;2, 1968. a
Devenish, B. J., Bartello, P., Brenguier, J.L., Collins, L. R., Grabowski, W. W., IJzermans, R. H. A., Malinowski, S. P., Reeks, M. W., Vassilicos, J. C., Wang, L.P., and Warhaft, Z.: Droplet growth in warm turbulent clouds, Q. J. Roy. Meteor. Soc., 138, 1401–1429, https://doi.org/10.1002/qj.1897, 2012. a
Dhaniyala, S. and Wexler, A. S.: Numerical schemes to model condensation and evaporation of aerosols, Atmos. Env., 30, 919–928, https://doi.org/10.1016/13522310(95)00288X, 1996. a
East, T. W. R.: An inherent precipitation mechanism in cumulus clouds, Q. J. Roy. Meteor. Soc., 83, 61–76, https://doi.org/10.1002/qj.49708335506, 1957. a, b, c, d, e
East, T. W. R. and Marshall, J. S.: Turbulence in clouds as a factor in precipitation, Q. J. Roy. Meteor. Soc., 80, 26–47, https://doi.org/10.1002/qj.49708034305, 1954. a
Feingold, G. and Chuang, P. Y.: Analysis of the Influence of FilmForming Compounds on Droplet Growth: Implications for Cloud Microphysical Processes and Climate, J. Atmos. Sci., 59, 2006–2018, https://doi.org/10.1175/15200469(2002)059<2006:AOTIOF>2.0.CO;2, 2002. a
Field, P., Heymsfield, A., Shipway, B., DeMott, P., Pratt, K., Rogers, D., Stith, J., and Prather, K.: Ice in Clouds Experiment–Layer Clouds. Part II: Testing Characteristics of Heterogeneous Ice Formation in Lee Wave Clouds, J. Atmos. Sci., 69, 1066–1079, https://doi.org/10.1175/JASD11026.1, 2012. a
Gettelman, A. and Morrison, H.: Advanced TwoMoment Bulk Microphysics for Global Models. Part I: OffLine Tests and Comparison with Other Schemes, J. Climate, 28, 1268–1287, https://doi.org/10.1175/JCLID1400102.1, 2015. a
Grabowski, W. and Wang, L.P.: Growth of Cloud Droplets in a Turbulent Environment, Annu. Rev. Fluid Mech., 45, 293–324, https://doi.org/10.1146/annurevfluid011212140750, 2013. a
Grabowski, W., Andrejczuk, M., and Wang, L.P.: Droplet growth in a bin warmrain scheme with Twomey CCN activation, Atmos. Res., 99, 290–301, https://doi.org/10.1016/j.atmosres.2010.10.020, 2011. a
Grabowski, W. W.: Comparison of Eulerian Bin and Lagrangian ParticleBased Microphysics in Simulations of Nonprecipitating Cumulus, J. Atmos. Sci., 77, 3951–3970, https://doi.org/10.1175/JASD200100.1, 2020. a
Hall, W. D.: A Detailed Microphysical Model Within a TwoDimensional Dynamic Framework: Model Description and Preliminary Results, J. Atmos. Sci., 37, 2486–2507, https://doi.org/10.1175/15200469(1980)037<2486:ADMMWA>2.0.CO;2, 1980. a
Hernandéz Pardo, L. H., Morrison, H., Machado, L. A. T., Harrington, J. Y., and Lebo, Z. J.: Drop Size Distribution Broadening Mechanisms in a Bin Microphysics Eulerian Model, J. Atmos. Sci., 77, 3249–3273, https://doi.org/10.1175/JASD200099.1, 2020. a, b, c
Hill, A., Shipway, B., and Boutle, I.: How sensitive are aerosolprecipitation interactions to the warm rain representation?, J. Adv. Model. Earth Syst., 7, 987–1004, https://doi.org/10.1002/2014MS000422, 2015.
Hill, R. N.: Numerical modelling of multimaterial interfaces, PhD thesis, Loughborough University, EThOS Persistent ID: uk.bl.ethos.545738, https://hdl.handle.net/2134/8103 (last access: April 2022), 2011. a, b
Hirt, C. W.: Heuristic stability theory for finitedifference equations, J. Comput. Phys., 2, 339–355, https://doi.org/10.1016/00219991(68)900417, 1968. a
Howell, W.: The growth of cloud drops in uniformly cooled air, J. Meteorol., 6, 134–149, https://doi.org/10.1175/15200469(1949)006<0134:TGOCDI>2.0.CO;2, 1949. a
Hulburt, H. and Katz, S.: Some problems in particle technology: A statistical mechanical formulation, Chem. Eng. Sci., 19, 555–574, https://doi.org/10.1016/00092509(64)850478, 1964. a
Jaruga, A., Arabas, S., Jarecka, D., Pawlowska, H., Smolarkiewicz, P. K., and Waruszewski, M.: libmpdata++ 1.0: a library of parallel MPDATA solvers for systems of generalised transport equations, Geosci. Model Dev., 8, 1005–1032, https://doi.org/10.5194/gmd810052015, 2015. a, b, c
Jeffery, C. A., Reisner, J. M., and Andrejczuk, M.: Another Look at Stochastic Condensation for Subgrid Cloud Modeling: Adiabatic Evolution and Effects, J. Atmos. Sci., 64, 3949–3969, https://doi.org/10.1175/2006JAS2147.1, 2007. a
Khain, A. P., Beheng, K. D., Heymsfield, A., Korolev, A., Krichak, S. O., Levin, Z., Pinsky, M., Phillips, V., Prabhakaran, T., Teller, A., van den Heever, S. C., and Yano, J.I.: Representation of microphysical processes in cloudresolving models: Spectral (bin) microphysics versus bulk parameterization, Rev. Geophys., 53, 247–322, https://doi.org/10.1002/2014RG000468, 2015. a
Kostinski, A. B. and Jameson, A. R.: On the Spatial Distribution of Cloud Particles, J. Atmos. Sci., 57, 901–915, https://doi.org/10.1175/15200469(2000)057<0901:OTSDOC>2.0.CO;2, 2000. a
Kostoglou, M. and Karabelas, A. J.: Evaluation of numerical methods for simulating an evolving particle size distribution in growth processes, Chem. Eng. Comm., 136, 177–199, https://doi.org/10.1080/00986449508936360, 1995. a
Kovetz, A.: An Analytical Solution for the Change of Cloud and Fog Droplet Spectra Due to Condensation, J. Atmos. Sci., 26, 302–304, https://doi.org/10.1175/15200469(1969)026<0302:AASFTC>2.0.CO;2, 1969. a
Kovetz, A. and Olund, B.: The Effect of Coalescence and Condensation on Rain Formation in a Cloud of Finite Vertical Extent, J. Atmos. Sci., 26, 1060–1065, https://doi.org/10.1175/15200469(1969)026<1060:TEOCAC>2.0.CO;2, 1969. a, b, c
Kühnlein, C. and Smolarkiewicz, P. K.: An unstructuredmesh finitevolume MPDATA for compressible atmospheric dynamics, J. Comp. Phys., 334, 16–30, https://doi.org/10.1016/j.jcp.2016.12.054, 2017. a
Kuo, H.C., Leou, T.M., and Williams, R. T.: A study on the highorder Smolarkiewicz methods, Comput. Fluids, 28, 779–799, https://doi.org/10.1016/s00457930(98)00036x, 1999. a
Lange, R.: ADPIC: a threedimensional computer code for the study of pollutant dispersal and deposition under complex conditions, LLNL report no. UCRL51462, https://doi.org/10.2172/4308175, 1973. a
Lange, R.: ADPIC–A ThreeDimensional ParticleinCell Model for the Dispersal of Atmospheric Pollutants and its Comparison to Regional Tracer Studies, J. Appl. Meteorol., 17, 320–329, https://doi.org/10.1175/15200450(1978)017<0320:ATDPIC>2.0.CO;2, 1978. a
Lee, H., Fridlind, A., and Ackerman, A.: An Evaluation of SizeResolved Cloud Microphysics Scheme Numerics for Use with Radar Observations. Part II: Condensation and Evaporation, J. Atmos. Sci., 78, 1629–1645, https://doi.org/10.1175/JASD200213.1, 2021. a, b, c
Li, X.Y., Brandenburg, A., Haugen, N. E. L., and Svensson, G.: Eulerian and Lagrangian approaches to multidimensional condensation and collection, J. Adv. Model. Earth Syst., 9, 1116–1137, https://doi.org/10.1002/2017MS000930, 2017. a, b
Liu, Q., Kogan, Y. L., Lilly, D. K., and Khairoutdinov, M. P.: Variational Optimization Method for Calculation of Cloud Drop Growth in an Eulerian DropSize Framework, J. Atmos. Sci., 54, 2493–2504, https://doi.org/10.1175/15200469(1997)054<2493:VOMFCO>2.0.CO;2, 1997. a, b, c
Lu, M.L. and Seinfeld, J.: Effect of aerosol number concentration on cloud droplet dispersion: A largeeddy simulation study and implications for aerosol indirect forcing, J. Geophys. Res., 111, D02207, https://doi.org/10.1029/2005JD006419, 2006. a
Margolin, L. G. and Shashkov, M.: MPDATA: gauge transformations, limiters and monotonicity, Int. J. Numer. Methods Fluids, 50, 1193–1206, https://doi.org/10.1002/fld.1070, 2006. a, b, c
Margolin, L. G. and Smolarkiewicz, P. K.: Antidiffusive Velocities for Multipass Donor Cell Advection, SIAM J. Sci. Comput., 20, 907–929, https://doi.org/10.1137/S106482759324700X, 1998. a, b, c, d
Margolin, L. G., Shashkov, M., and Smolarkiewicz, P. K.: A discrete operator calculus for finite difference approximations, Comput. Methods Appl. Mech. Eng., 187, 365–383, https://doi.org/10.1016/S00457825(00)800018, 2000. a
Morrison, H., Witte, M., Bryan, G. H., Harrington, J. Y., and Lebo, Z. J.: Broadening of Modeled Cloud Droplet Spectra Using Bin Microphysics in an Eulerian Spatial Domain, J. Atmos. Sci., 75, 4005–4030, https://doi.org/10.1175/JASD180055.1, 2018. a, b, c, d, e, f, g, h, i
Morrison, H., van LierWalqui, M., Fridlind, A. M., Grabowski, W. W., Harrington, J. Y., Hoose, C., Korolev, A., Kumjian, M. R., Milbrandt, J. A., Pawlowska, H., Posselt, D. J., Prat, O. P., Reimel, K. J., Shima, S.I., van Diedenhoven, B., and Xue, L.: Confronting the Challenge of Modeling Cloud and Precipitation Microphysics, J. Adv. Model. Earth Syst., 12, e2019MS001689, https://doi.org/10.1029/2019MS001689, 2020. a
Onishi, R., Sugimura, T., and Takahashi, K.: CIPCSLR Scheme for Condensation and Evaporation Calculations of Water Droplets, J. Environ. Eng., 5, 1–14, https://doi.org/10.1299/jee.5.1, 2010. a
Ramkrishna, D.: Population Balances: Theory and Applications to Particulate Systems in Engineering, Academic Press, ISBN 9780125769709, https://doi.org/10.1016/B9780125769709.X50000, 2000. a
Rauber, R., Stevens, B., Ochs III, H., Knight, C., Albrecht, B., Blyth, A., Fairall, C., Jensen, J., LasherTrapp, S., MayolBracero, O., Vali, G., Anderson, J., Baker, B., Bandy, A., Burnet, F., Brenguier, J.L., Brewer, W., Brown, P., Chuang, P., Cotton, W., Di Girolamo, L., Geerts, H., Gerber, H., Göke, S., Gomes, L., Heikes, B., Hudson, J., Kollias, P., Lawson, R., Krueger, S., Lenschow, D., Nuijens, L., O'Sullivan, D., Rilling, R., Rogers, D., Siebesma, A., Snodgrass, E., Stith, J., Thornton, D., Tucker, S., Twohy, C., and Zuidema, P.: Rain in Shallow Cumulus Over the Ocean: The RICO Campaign, B. Am. Meteorol. Soc., 88, 1912–1928, https://doi.org/10.1175/BAMS88121912, 2007. a
Roberts, K. V. and Weiss, N. O.: Convective Difference Schemes, Math. Comput., 20, 272–299, https://doi.org/10.2307/2003507, 1966. a
Schneider, T. Teixeira, J., Bretherton, C. S., Brient, F., Pressel, K. G., Schär, C., and Siebesma, A. P.: Climate goals and computing the future of clouds, Nat. Clim. Change, 7, 3–5, https://doi.org/10.1038/nclimate3190, 2017. a
Shipway, B. and Hill, A.: Diagnosis of systematic differences between multiple parametrizations of warm rain microphysics using a kinematic framework, Q. J. Roy. Meteor. Soc., 138, 2196–2211, https://doi.org/10.1002/qj.1913, 2012. a, b, c, d
Slawinska, J., Grabowski, W., Pawlowska, H., and Morrison, H.: Droplet Activation and Mixing in LargeEddy Simulation of a Shallow Cumulus Field, J. Atmos. Sci., 69, 444–462, https://doi.org/10.1175/JASD11054.1, 2012. a
Smolarkiewicz, P. and Margolin, L.: MPDATA: A FiniteDifference Solver for Geophysical Flows, J. Comp. Phys., 140, 459–480, https://doi.org/10.1006/jcph.1998.5901, 1998. a, b
Smolarkiewicz, P. K.: A Simple Positive Definite Advection Scheme with Small Implicit Diffusion, Mon. Weather Rev., 111, 479–486, https://doi.org/10.1175/15200493(1983)111<0479:ASPDAS>2.0.CO;2, 1983. a, b, c, d, e
Smolarkiewicz, P. K.: A Fully Multidimensional Positive Definite Advection Transport Algorithm with Small Implicit Diffusion, J. Comp. Phys., 54, 325–362, https://doi.org/10.1016/00219991(84)901219, 1984. a, b, c, d, e, f, g
Smolarkiewicz, P. K.: Multidimensional positive definite advection transport algorithm: An overview, Int. J. Numer. Methods Fluids, 50, 1123–1144, https://doi.org/10.1002/fld.1071, 2006. a, b
Smolarkiewicz, P. K. and Clark, T. L.: The Multidimensional Positive Definite Advection Transport Algorithm: Further Development and Applications, J. Comp. Phys., 67, 396–438, https://doi.org/10.1016/00219991(86)902706, 1986. a, b, c, d
Smolarkiewicz, P. K. and Grabowski, W. W.: The Multidimensional Positive Definite Advection Transport Algorithm: Nonoscillatory Option, J. Comp. Phys., 86, 355–375, https://doi.org/10.1016/00219991(90)90105A, 1990. a, b, c, d
Smolarkiewicz, P. K. and Margolin, L. G.: On ForwardinTime Differencing for Fluids: Extension to a Curvilinear Framework, Mon. Weather Rev., 121, 1847–1859, https://doi.org/10.1175/15200493(1993)121<1847:OFITDF>2.0.CO;2, 1993. a
Smolarkiewicz, P. K. and Margolin, L. G.: MPDATA – A Multipass Donor Cell Solver for Geophysical Flows, in: Godunov methods: Theory and applications, edited by: Toro, E., Springer, https://doi.org/10.1007/9781461506638_81, 2001. a
Smolarkiewicz, P. K. and Rasch, P. J.: Monotone Advection on the Sphere: An Eulerian Versus SemiLagrangian Approach, J. Atmos. Sci., 48, 793–810, https://doi.org/10.1175/15200469(1991)048<0793:MAOTSA>2.0.CO;2, 1991. a, b, c
Smolarkiewicz, P. K. and Szmelter, J.: MPDATA: An edgebased unstructuredgrid formulation, J. Comp. Phys., 206, 624–649, https://doi.org/10.1016/j.jcp.2004.12.021, 2005. a
Toro, E.: Riemann Solvers and Numerical Methods for Fluid Dynamics, A Practical Introduction, Springer, 3 edn., https://doi.org/10.1007/b79761, 2009. a
Tsang, T. H. and Brock, J. R.: Simulation of Condensation Aerosol Growth by Condensation and Evaporation, Aerosol Sci. Tech., 2, 311–320, https://doi.org/10.1080/02786828308958637, 1982. a, b
Tsang, T. H. and Korgaonkar, N.: Effect of Evaporation on the Extinction Coefficient of an Aerosol Cloud, Aerosol Sci. Tech., 7, 317–328, https://doi.org/10.1080/02786828708959167, 1987. a
Tsang, T. H. and Rao, A.: Comparison of Different Numerical Schemes for Condensational Growth of Aerosols, Aerosol Sci. Tech., 9, 271–277, https://doi.org/10.1080/02786828808959214, 1988. a, b, c
Tsang, T. H. and Rao, A.: A moving finite element method for the population balance equation, Num. Meth. Fluids, 10, 753–769, https://doi.org/10.1002/fld.1650100704, 1990. a
Waruszewski, M., Kühnlein, C., Pawlowska, H., and Smolarkiewicz, P. K.: MPDATA: Thirdorder accuracy for variable flows, J. Comput. Phys., 359, 361–379, https://doi.org/10.1016/j.jcp.2018.01.005, 2018. a, b, c
Wei, L., Sun, J., Lei, H., Dong, L., and Hu, W.: A Lagrangian Advection Scheme for Solving Cloud Droplet Diffusion Growth, Atmosphere, 11, 632, https://doi.org/10.3390/atmos11060632, 2020. a, b
Williams, M. M. R. and Loyalka, S. K.: Aerosol Science: Theory and Practice: With Special Applications to the Nuclear Industry, Oxford, Pergamon Press, 1991. a
Yang, F., Kollias, P., Shaw, R. A., and Vogelmann, A. M.: Cloud droplet size distribution broadening during diffusional growth: ripening amplified by deactivation and reactivation, Atmos. Chem. Phys., 18, 7313–7328, https://doi.org/10.5194/acp1873132018, 2018. a
 Abstract
 Introduction
 Spectral advection with upwind and MPDATA solutions (boxmodel test case)
 Spectralspatial advection with MPDATA (singlecolumn test case)
 Conclusions
 Appendix A: Convergence analysis
 Code availability
 Data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
 Abstract
 Introduction
 Spectral advection with upwind and MPDATA solutions (boxmodel test case)
 Spectralspatial advection with MPDATA (singlecolumn test case)
 Conclusions
 Appendix A: Convergence analysis
 Code availability
 Data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References