# Eulerian and modified Lagrangian approaches to multi-dimensional condensation and collection

###### Abstract

Turbulence is argued to play a crucial role in cloud droplet growth. The combined problem of turbulence and cloud droplet growth is numerically challenging. Here, an Eulerian scheme based on the Smoluchowski equation is compared with two Lagrangian superparticle (or superdroplet) schemes in the presence of condensation and collection. The growth processes are studied either separately or in combination using either two-dimensional turbulence, a steady flow, or just gravitational acceleration without gas flow. Good agreement between the different schemes for the time evolution of the size spectra is observed in the presence of gravity or turbulence. Higher moments of the size spectra are found to be a useful tool to characterize the growth of the largest drops through collection. Remarkably, the tails of the size spectra are reasonably well described by a gamma distribution in cases with gravity or turbulence. The Lagrangian schemes are generally found to be superior over the Eulerian one in terms of computational performance. However, it is shown that the use of interpolation schemes such as the cloud-in-cell algorithm is detrimental in connection with superparticle or superdroplet approaches. Furthermore, the use of symmetric over asymmetric collection schemes is shown to reduce the amount of scatter in the results.

Journal of Advances in Modeling Earth Systems

Department of Meteorology, and Bolin Centre for Climate Research, Stockholm University, Stockholm, Sweden Nordita, KTH Royal Institute of Technology and Stockholm University, 10691 Stockholm, Sweden Swedish e-Science Research Centre, www.e-science.se, Stockholm, Sweden Laboratory for Atmospheric and Space Physics, University of Colorado, Boulder, CO 80303, USA JILA and Department of Astrophysical and Planetary Sciences University of Colorado, Boulder, CO 80303, USA Department of Astronomy, Stockholm University, SE-10691 Stockholm, Sweden SINTEF Energy Research, 7465 Trondheim, Norway Department of Energy and Process Engineering, NTNU, 7491 Trondheim, Norway Global & Climate Dynamics, National Center for Atmospheric Research, Boulder, CO 80305, USA

Xiang-Yu L, Revision: 1.715

Eulerian Smoluchowski and Lagrangian superdroplet/superparticle approaches to cloud droplet growth through condensation and collection are compared using DNS techniques

Size spectra agree well for both approaches, especially in case of turbulence

The Lagrangian scheme with symmetric collection is found to be optimal and computationally most efficient

## 1 Introduction

In the context of raindrop formation, it is generally accepted that turbulence plays a crucial role in bridging the size gap between efficient condensational growth of small particles (radii below ) and efficient collectional growth due to gravity of larger ones (radii around and above) (Shaw, 2003; Grabowski and Wang, 2013; Khain et al., 2007). Improving the understanding of this important problem in meteorology (Berry and Reinhardt, 1974; Pinsky and Khain, 1997; Falkovich et al., 2002; Naumann and Seifert, 2016a) might also shed light on how to bridge the even more severe size gap in the astrophysical context of planetesimal formation (Johansen et al., 2007, 2012). To address these questions numerically, one has to combine direct numerical simulations (DNS) of turbulent gas motions with those of particles. The particles are cloud droplets in the meteorological context and dust grains in astrophysics. A possible approach to treat collection is to solve the Smoluchowski equation (also known as the stochastic collection equation in the meteorological context) (Ogura and Takahash, 1973; Svensson and Seinfeld, 2002; Bec et al., 2016), which couples the spatio-temporal evolution equations of the particle distribution function for different particle sizes. The particle motion can be treated using a fluid description for each particle size. Thus, not only does one have to solve the Smoluchowski equation at each meshpoint, but, because heavier particles have finite momenta and speeds that are different from those of the gas, one has to solve corresponding momentum equations for each mass species. In the meteorological context, it is also referred to as a binned spectral method, although in that case the momentum equations for the particle bins are normally ignored (Xue et al., 2008). An Eulerian approach is technically more straightforward than a Lagrangian one, but it becomes computationally demanding as the size range of cloud droplets is large.

The Eulerian approach also has conceptual difficulties if the collection probability depends on the mutual velocity difference. This is due to the fact that particles of the same size are described by the same momentum equation and have therefore the same velocity at a given position in space, so the velocity difference vanishes. This means that particles of the same size are not allowed to collide. This is not a problem for freely falling particles of the same size, which would have the same terminal velocity and are not expected to collide. This would however be an unrealistic restriction when particles are subject to acceleration by turbulence. More importantly, as was emphasized in the recent review of Khain et al. (2015), the Smoluchowski equation is a mean-field equation and cannot capture the random properties of the collections if the collision kernel is prescribed a priori. Nevertheless, most numerical cloud microphysical approaches are based on the Smoluchowski equation, which therefore raises questions regarding the accuracy of the basic equations Khain et al. (2015). Thus, new approaches based on inherently different equations are required to model the cloud microphysical processes.

An alternative approach is the Lagrangian one, where one solves for the motion of individual particles and treats collections explicitly. In atmospheric clouds, the number density of micrometer-sized cloud droplets is of the order of , so in a volume of , one has 100 million particles, which is the typical size that can be treated on modern supercomputers. A domain of this size is also about the largest that is possible in direct numerical simulations (DNS) of atmospheric turbulence; the Reynolds number based on the length scale and the corresponding velocity scale is , where is the viscosity of the gas flow. Such a large Reynolds number is just within reach on current supercomputers, but larger domains would remain out of reach for a long time. Several earlier works investigated condensational growth of cloud droplets using Lagrangian tracking in DNS (Paoli and Shariff, 2009; Sardina et al., 2015; de Lozar and Muessle, 2016; Lanotte et al., 2009), but those neglected the collectional growth and only proposed to study the collectional growth in future work. An intermediate approach involves the use of Lagrangian “superparticles” (Johansen et al., 2012; Pruppacher and Klett, 2012; Shima et al., 2009; Zsom and Dullemond, 2008), which represent a “swarm” of particles of certain size and number density. Depending on the values of particle size and number density, there is a certain probability that an encounter between two superparticles leads to collectional growth of some of the particles in each swarm (or superparticle). This superparticle approach has been applied in a recent LES model to represent the cloud microphysical condensation (Andrejczuk et al., 2008) and collection (Andrejczuk et al., 2010; Riechelmann et al., 2012; Naumann and Seifert, 2015) processes.

The purpose of the present paper is to compare the Eulerian approach involving the Smoluchowski equation with the Lagrangian superparticle approach with the aim of identifying a promising DNS scheme for tackling the bottleneck problem of cloud droplets growth. This has been done in the astrophysical context (Ohtsuki et al., 1990; Dra̧żkowska et al., 2014), where the principal problem with the Eulerian approach was emphasized in that it requires high mass bin resolution (MBR) to avoid artificial speedup of the growth rate. Here we also compare with the superdroplet approach of Shima et al. (2009). The original work on this approach was restricted to the case of vanishing particle inertia, but this restriction is not a principal limitation of this scheme, which is in fact well applicable to the case of finite particle inertia.

## 2 Lagrangian and Eulerian approaches

In the following, we refer to the superparticle or superdroplet approaches as the swarm model, where each superparticle represents a swarm of physical particles. By contrast, the Eulerian approach is also referred to as the Smoluchowski model. Here we compare the two approaches in the meteorological context of water droplets using, however, simplifying assumptions such as constant supersaturation and ideal collection efficiency. In this paper, we generally refer to particles and superparticles, which are thus used interchangeably with droplets and superdroplets, respectively. We begin with a discussion of the gas flows that are being used in some of the models.

### 2.1 Evolution equations for the gas flow in both approaches

In all the experiments reported below, where a nonvanishing gas flow is used, we restrict ourselves to two-dimensional (2-D) flows. However, we also perform several experiments with no gas flow (). In those cases the system is spatially uniform and therefore zero-dimensional (0-D). For the swarm model, however, each swarm occupies one grid cell, so it must be treated in at least one dimension (1-D), although the results of higher-dimensional swarm models will also be discussed, and they are computational cheaper because they can take advantage of better parallelization. By comparison, the Eulerian models are strictly 0-D when there is no flow.

#### 2.1.1 Momentum equation of the gas flow

To obtain at each meshpoint, we solve the usual Navier-Stokes equation

(1) |

where is a forcing term, is the gas pressure, is the gas density, which in turn obeys the continuity equation,

(2) |

the viscous force is given by

(3) |

where so that the pressure is proportional to the gas density . Note that gravity has been neglected in equation (1), but this is not a principal restriction and can be relaxed once suitable non-periodic boundary conditions are adopted. For the relatively small domains that can be handled by DNS, gravity will nevertheless have only minor effects on the fluid flow for atmospheric conditions. is the traceless rate-of-strain tensor and commas denote partial differentiation. We assume that the gas is isothermal and has constant sound speed

#### 2.1.2 Straining flow

To obtain a non-vanishing flow, we apply volume forcing via the term . In the case of a time-independent 2-D divergence-free straining flow,

(4) |

we take , where determines the amplitude and the wavenumber of the flow.

#### 2.1.3 Turbulence

In the case of a turbulent flow, is delta-correlated in time and consists of random waves in space (Haugen et al., 2004). The flow is characterized by a typical forcing wavenumber ( for the straining flow or the average wavenumber from a narrow band of wavevectors) and the root-mean-square (rms) velocity . As a relevant timescale characterizing such a flow, we define

(5) |

which is an estimate of the correlation time. This definition is also used for the straining flow, which is a special case in that it is time-independent and therefore would no longer characterize the correlation time of the flow, but it would still be proportional to the turnover time. A simulation without spatial extent can be adopted to investigate the statistical convergence properties of the Eulerian model regarding its computational efficiency.

### 2.2 Condensational growth

The growth of the particle radius by condensation is governed by (Lamb and Verlinde, 2011)

(6) |

where is the supersaturation and is the condensation parameter. Both and are in principle dependent on the flow and the environmental temperature and pressure (see Chapter 8 of Lamb and Verlinde, 2011), but these dependencies are here neglected, because it would complicate the comparison of different numerical schemes even further. Therefore, the condensational growth is driven by constant water vapor flux without latent heat release in the present study. We adopt the value (Lanotte et al., 2009). The assumed constancy of also implies that the total liquid water content is not conserved.

### 2.3 The swarm model

The swarm model is a Monte Carlo type approach that handles particle collections in a swarm of particles in a statistical manner (Zsom and Dullemond, 2008). Each swarm has a particle number density , and occupies a volume , which equals the volume of a fluid grid cell of size in dimensions. All particles in a given swarm have the same mass, radius, and velocity. Following the description of Johansen et al. (2012), the swarm is transported along with its “shepherd particle”, which is also referred to as the corresponding superparticle. The swarm is treated as a Lagrangian point-particle, where one solves for the particle position via

(7) |

and the velocity via

(8) |

in the usual way. Here, is the gravitational acceleration, is the particle inertial response or stopping time of a particle in swarm and is given by

(9) |

where is the radius of particles in swarm , is the particle solid material density, is the density of the gas and the effective viscosity is given by Sullivan et al. (1994)

(10) |

where is the ordinary (microphysical) fluid viscosity, and is the particle Reynolds number, which provides a correction factor to the particle stopping time.

A given swarm may only interact with every other swarm within the same grid cell. The computational cost associated with such collections scales as , where is the number of swarms within a grid cell, but this is computationally not prohibitive as long as is not too large.

We now consider two swarms and residing within the same grid cell. Consider first collections of particles within swarm with a particle of swarm . The inverse mean free path of in is given by

(11) |

where is the collectional cross section with

(12) |

and is the collision efficiency, but in the following we assume
in all cases.^{1}^{1}1
In Shima et al. (2009), the mean free path
is defined by invoking the swarm with the larger
number density of
physical particles; see section 2.3.2 for details.
The particle number density in swarm is
and and represent the radii of the particles in the
two swarms.
From this, one can find the typical rate of collections between a
particle of swarm and particles of swarm as

(13) |

where and are the velocities of swarms and . The probability of a collection between the swarm and any of the particles of swarm within the current time step is then given by

(14) |

This effectively puts a restriction on the time step, since the probability cannot be larger than unity. For each swarm pair in a grid cell, one now picks a random number, , and compares it with . A collection event occurs in the case when .

#### 2.3.1 Collection scheme I

For the swarm model, two different collection schemes have been proposed in the astrophysical and meteorological contexts. We begin discussing the former (scheme I), which is similar to that described by Johansen et al. (2012) in that it maintains a constant mass of the individual swarms. In the context of mathematical probability, this approach is also known as mass flow algorithm (Eibeck and Wagner, 2001; Patterson et al., 2011). Scheme II is discussed in section 2.3.2.

If , one assumes that *all* the particles in swarm
have collided with a particle in swarm . In this collection scheme,
all swarms are treated individually. This means that even though the
particles in swarm have
collided with the particles in swarm , swarm is kept unchanged at this
stage. Instead, swarm is treated individually at a different stage.
Hence, all collections are asymmetric, i.e., .
The new mass of the particles in swarm now becomes

(15) |

where is the mass before the collection and the tilde represents the new value after collection. In order to ensure mass conservation, the total mass of swarm is kept unchanged, i.e.,

(16) |

which implies that the new particle number density, , is given by ; see equation (17) of Patterson et al. (2011) for the corresponding treatment in the mass flow algorithm. By invoking momentum conservation,

(17) |

the new velocity of any particle in swarm is given by .

#### 2.3.2 Collection scheme II

In the meteorological context, the following collection scheme has been proposed (Shima et al., 2009). Assume two swarms and , and consider (without loss of generality) the case . The collection probability of particles in swarm with swarm is, again, given by equation (14). If the two swarms are found to collide, the new masses of the particles in the two swarms are given by

(18) |

but now their new particle number densities are

(19) |

In other words, the number of particles in the smaller swarm remains unchanged (and their masses are increased), while that in the larger one is reduced by the amount of particles that have collided with all the particles of the smaller swarm (and their masses remain unchanged). This implies that in equation (11), the mean free path is defined with respect to the swarm with the larger number density of physical particles, as explained in Shima et al. (2009). Finally, the new momenta of the particles in the two swarms are given by

(20) |

In contrast to scheme I, these collections are symmetric, i.e. . Consequently, both swarms are changed during a collection. However, the asymmetric collection property of scheme I of (Johansen et al., 2012) may not have been previously recognized, nor has its accuracy been compared with other models, which we will further discuss below.

#### 2.3.3 Initial particle distribution

We recall that particles within a swarm may interact with particles of another swarm only if both swarms occupy the same grid cell. The effective volume of each swarm is therefore equal to , where is the spatial dimension introduced in section 2.3. The total number of particles in our computational domain is therefore times the sum of over all swarms. This must also be equal to , where is the total number density represented by the simulation and is the size of the computational domain. Thus, we have

(21) |

Initially (), the particle number densities of all swarms are the same and since is the total number of grid points, we have . Thus, the initial number density of particles within one swarm must be

(22) |

In the following, we choose the initial particle size distribution of total physical particles in the domain to be log-normal, i.e.,

(23) |

where and are the center and width of the size distribution, respectively; is the initial total number density of physical particles. These particles are distributed uniformly over all swarms within the computational domain. This means that particles in each swarm are of the same size, but different from swarm to swarm.

### 2.4 Eulerian approach

To model the combined growth of particles through condensation and collection in a multi-dimensional flow in the Eulerian description, we describe the evolution of particles of different radii (or, equivalently, of different logarithmic particle mass ) at different positions and time . We employ the particle distribution function , or, alternatively in terms of logarithmic particle mass , , such that the total number density of particles is given by

(24) |

or, correspondingly for , we have . Since , we have . Note that obeys the usual continuity equation,

(25) |

where is the mean particle velocity (i.e., an average over all particle sizes) and is a Brownian diffusion term, which is enhanced for numerical stability and will be chosen depending on the mesh resolution. The evolution of the particle distribution function is governed by a similar equation, but with additional coupling terms due to condensation and collection, i.e.

(26) |

where is the derivative with respect to , , as given in equation (6), and describes the change of the number density of particles for smaller and larger radii, as will be defined below. Furthermore, is the particle velocity within the resolved grid cell, which is discussed below. It also determines the mean flow in equation (25).

The modeling of condensation and collection implies coupling of the evolution equations of for different values of . The advantage of using is that it allows us to cover a large range in , because we will use then an exponentially stretched grid in such that is uniformly spaced (Pruppacher and Klett, 2012; Suttner and Yorke, 2001; Johansen, 2004). The total number density within a finite mass interval is then given by . Thus, the total number density of particles of all sizes at position and time is given by

(27) |

where is the variable used in the simulations and is the number of logarithmic mass bins.

Let us first consider the process of condensation, which is described in equation (26) by the term , where is the flux of particle from one size bin to the next. Evidently, the total number density is only conserved if the particle flux vanishes for and , which is the case if the range of is sufficiently large. In particular, , because for . In practice, however, we consider finite lower cutoff values of and therefore expect some degree of mass loss at the smallest mass bins. The same is also true for the largest mass bin once the size distribution has grown to sufficiently large values. In all cases with pure condensation, it is convenient to display solutions in non-dimensional form by measuring time in units of

(28) |

and in units of . We refer to Appendix A for more details on the condensation equation for the Eulerian approach.

Next, we consider collection, which leads to a decrease of , but does not change the mean mass density of liquid water. The evolution of due to collection is governed by the Smoluchowski equation

(29) |

Here, is a kernel, which is proportional to the collision efficiency and a geometric contribution. As mentioned above, we assume and so is given by

(30) |

where and are the radii of the corresponding mass variables, and , while and are their respective velocities, whose governing equation is given below.

In the following, we define the mass and radius bins such that

(31) |

Unfortunately, is in many cases far too coarse, so we take

(32) |

where is a parameter that we chose to be a power of two. For a fixed mass bin range, the number of mass bins increases with increasing . In terms of , equation (29) reads

(33) |

where we have adopted the nomenclature of Johansen (2004), where denotes all values of and for which

(34) |

is fulfilled. The term in equation (33) comes from the fact that collections between cloud droplets from two mass bins may not necessarily result in a cloud droplet mass being exactly in the middle of the nearest mass bin. Johansen (2004) therefore included this factor so that mass is strictly conserved. The discrete kernel is then .

The corresponding momentum equations for the velocities for each logarithmic mass value is

(35) |

Here, is the gas velocity, (for ) is defined by equation (9), and

(36) |

is a viscous force among particles, which should be very small for dilute particle suspensions, but is nevertheless retained in equation (35) for the sake of numerical stability of the code. It is not to be confused with the drag force, between particles and gas. In principle, the expression for should be based on the divergence of the traceless rate-of-strain tensor of , similarly to the corresponding expression for the viscous force of the gas discussed in equation (3). However, since the term is unphysical anyway, we just use the simpler expression proportional to instead.

The linear momentum of all particles is given by , where angle brackets denote volume averages. In order that this quantity is conserved by each collection, the target has to receive a corresponding kick, which leads to the last term in equation (35), but it leaves the velocities of the collection partners unchanged. It is therefore only related to the first term on the right-hand side of equation (33) and not the second, so it is given by (see Appendix B)

(37) |

To our knowledge, this momentum-conserving term has not been included in any of the very few earlier works that include a momentum equation for each particle species (cf. Suttner and Yorke, 2001; Elperin et al., 2015). The reason why this has apparently not previously been discussed in the literature is that in meteorological applications one usually works with the averaged kernel and neglects the evolution of the velocities for the different mass bins (Grabowski and Wang, 2013). This correction term is evidently zero when the momentum of the two collection constituents () is equal to that of the resulting constituent []. Nevertheless, as is shown in Appendix B, the momentum conserving correction changes the time evolution of the droplet spectrum in an unexpected way when the MBR is high, but the results are similar for . Furthermore, for turbulent flows, as is discussed below, these correction terms become insignificant.

As mentioned above, a shortcoming of the Eulerian approach is that no collection is possible from equally sized particles. To assess the consequences of this unphysical limitation, we study the sensitivity of the results to replacing either (i) by or (ii) by , where is an empirical parameter.

### 2.5 Boundary conditions and diagnostics

In the present work, we use periodic boundary conditions for all variables in all directions. Therefore, no particles and no gas are lost through the boundaries of the domain. This approximation is reasonable as long as we are interested in modeling a small domain well within a cloud where also heavier particles can be assumed to enter from above. The use of periodic boundary conditions requires us to neglect gravity in equation (1), which could be relaxed if non-periodic boundary conditions were adopted.

To characterize the size distribution, especially for the larger particles, we consider the evolution of different normalized moments of the size spectra,

(38) |

where is a positive integer. The mean radius is given by . Higher moments represent the tail of the distribution at large radii. In view of raindrop formation, we will be particularly interested in the largest droplets in the distribution. However, very large moments become numerically difficult to compute accurately, but , for example, was still not sufficiently representative of the largest droplets. Therefore we arrived at as a reasonable compromise to characterize the largest droplets in the distribution. Alternatively, the size distribution can be characterized by a gamma distribution, which requires the determination of only three moments in an approach known as the three-moment bulk scheme (Seifert and Beheng, 2001). This will be discussed in more detail in section 3.3.

In the case of collection, the condensation timescale , defined in equation (28), is no longer relevant, but it is instead a collection timescale that can be defined in the Eulerian model as

(39) |

which is, in this definition, a time-dependent quantity. In the Lagrangian model, this quantity can be defined by the collection frequency. Unlike the case of pure condensation, where is the appropriate time unit, can only be used a posteriori as a diagnostic quantity. However, given that the speed of pure collection is proportional to the mean particle density , it is often convenient to perform simulations at increased values of and then rescale time to a fixed reference density and use

(40) |

In the following we use , which is the typical value of in atmospheric clouds. Analogously we also define Finally, the number of particles in the total simulation domain is .

### 2.6 Computational implementation

We use the Pencil Code^{2}^{2}2https://github.com/pencil-code/,
which is a public domain code where
the relevant equations have been implemented (Johansen, 2004; Johansen et al., 2004; Babkovskaia et al., 2015).
We refer to Appendix A for a description of an important modification
applied to the implementation of equation (6).
The implementation of equation (33) has been discussed in detail by
Johansen (2004), and follows an approach described earlier by Suttner and Yorke (2001).
However, momentum conservation during collections
was previously ignored in the Eulerian model.
The current revision number is 73563 when checking out the code via the
svn bridge on the public github repository.

When traditional point particle Lagrangian particle tracking is employed,
it is usually beneficial to employ higher order interpolation between the
neighboring grid cells to find the value of a given fluid variable at
the exact position of the particle.
By default, the cloud-in-cell (CIC) algorithm is used, which involves
first order interpolation for the particle properties on the mesh.
In the swarm approach, however, the particles in each swarm fills
the volume of a grid cell in which the shepherd particle is.
The distribution of the swarm throughout the grid cell is homogeneous and
isotropic, and as such the swarm has no particular position within the grid cell.
It is true that there is a particular position associated with the
swarm, namely the position of the shepherd particle, but this position
has no purpose other than to determine in which grid cell the swarm resides.
Below we shall show that it is *n*ot better to use any kind
of interpolation in determining the value of the fluid variables at the
position of the swarm, but rather to use the values of the grid cell in
which the swarm resides.
This method is technically referred to as nearest grid point mapping (NGP).
Details concerning each experiment are summarized in Table 1.

## 3 Results

### 3.1 Condensation experiments

Run | Scheme | Dim | (m) | IM | Processes | () | Flow | () | () | |||

1A | SwI | 3-D | CIC | Con | – | – | ||||||

2B | Eu | 0-D | – | – | – | Con | – | |||||

3C | Eu | 0-D | – | – | – | Col | grav | |||||

4B | SwI | 3-D | CIC | Col | – | grav | ||||||

6A | SwII | 3-D | CIC | Col | – | grav | ||||||

7A | SwII | 2-D | CIC | Col | – | strain | ||||||

7B | SwII | 2-D | CIC | Col | – | strain | ||||||

7C | SwII | 2-D | CIC | Col | – | strain | ||||||

7D | Eu | 2-D | – | – | Col | strain | ||||||

7E | SwII | 2-D | NGP | Col | – | strain | ||||||

7F | SwII | 2-D | NGP | Col | – | strain | ||||||

9A | SwII | 2-D | CIC | Both | – | strain | ||||||

9B | SwII | 2-D | NGP | Both | – | strain | ||||||

9C | Eu | 2-D | – | – | Both | strain | ||||||

9D | Eu | 2-D | – | – | Both | strain | ||||||

9E | Eu | 2-D | – | – | Both | strain | ||||||

10A | Eu | 2-D | – | – | Col | turb | ||||||

10B | SwII | 2-D | NGP | Col | – | turb | ||||||

“IM” denotes the interpolation method, “Col” refers to collection, “Con” refers to condensation, “Eu” refers to Eulerian model, “SwI” refers to collection scheme I of swarm model, “SwII” refers to collection scheme II of swarm model, “Both” refers to condensation and collection, “grav” refers to gravity (=0), “strain” refers to straining flow, “turb” refers to turbulence, and “Dim” refers to the dimension. |

We compare the Eulerian and Lagrangian models for the pure condensation process without motion, i.e., zero gas velocity. In the case of homogeneous condensation, we can compare the numerical solution with the analytic solution of Seinfeld and Pandis (2006); see their Fig. 13.25. To this end, we make use of the fact that solutions of the condensation equation (6) obey

(41) |

where is a shifted coordinate with . With the log-normal initial distribution given by equation (23), this yields

(42) |

where denotes the position of the peak of the distribution and denotes its width, where is the symbol introduced by Seinfeld and Pandis (2006). What is remarkable here is the fact that vanishes for . This is because in this model, no new particles are created and even particles of zero initial radius will have grown to a radius after time . Furthermore, the small particles with grow faster than any of the larger ones, which leads to a sharp rise in the distribution function at . Thus, has a discontinuity at . This poses a challenge for the Eulerian scheme in which the derivative is discretized; see equation (26). In Figure 1, we compare solutions obtained using both Eulerian and Lagrangian approaches. It is evident that the -dependence obtained from the Eulerian solution is too smooth compared with the analytic one, even though we have used 1281 mass bins with =128 to represent on our logarithmically spaced mesh over the range , which corresponds to ; see equation (32). Better accuracy could be obtained by using a uniformly spaced grid in , but this would not be useful later when the purpose is to consider collection spanning a range of several orders of magnitude in radius. By comparison, the Lagrangian solution shown in the right-hand panel of Figure 1 (here with ) has no difficulty in reproducing the discontinuity in at . Moreover, the Lagrangian solution agrees perfectly with the analytical solution.

In practice, for the Eulerian approach we would use logarithmic spacing on a mesh with or . However, in such cases, the distribution develops a broad tail. This is demonstrated in detail in Appendix C. On the other hand, as we show further below, for turbulent and other velocity fields, the results depend much less on MBR so that computations with can be sufficiently accurate.

### 3.2 Purely gravitational collection experiments

We now consider uniform collection with no spatial variation of the velocity and density fields for both the gas and the particles. For the purely geometrical kernel, no analytic solution exists. However, we can compare the convergence properties of our two quite different numerical approaches and thereby get some sense of their validity in cases when the two agree. We consider pure collection experiments, starting again with a log-normal distribution. The results are presented in terms of normalized time; see equation (40).

#### 3.2.1 Comparison between swarm collection schemes I and II

In Figure 2, we compare schemes I and II of the swarm model together with the Eulerian model. The simulations have been performed with grid points and swarms (the statistics is converged for , as discussed in Appendix D). The left-hand panel of Figure 2 shows that for the results of the swarm simulations with scheme I agree with those of scheme II at early times, but depart at late times. However, for , the agreement is excellent, as shown in the right-hand panel of Figure 2. The evolution of with scheme I shows considerable scatter at late times. We recall that the main difference between schemes I and II is the geometry of collections. The collections simulated with scheme I are asymmetric, while those with scheme II are symmetric. Thus, in scheme II both swarms change either their total mass or their total particle number, while in scheme I the total mass of a swarm is kept constant by adjusting the particle number correspondingly. This property of scheme I may be responsible for creating stronger fluctuations in the mean radius. Therefore, to keep the amount of scatter comparable, scheme II is effectively less demanding. In the following, we will mainly adopt scheme II to save computational time.

#### 3.2.2 Comparison between collection scheme II and the Eulerian model

As we have seen above, the swarm simulations follow the Eulerian results rather well for (see the right-hand panel of Figure 2), but are somewhat different for . At early times, on the other hand, the evolution of obtained with the swarm model with collection scheme I follows more closely that of the Eulerian model. However, at later times, the evolution of obtained with the swarm model departs from the one simulated with the Eulerian model. This is surprising and might hint at a false convergence behavior, especially of the swarm model which has very few particles at small radii. This interpretation is supported by the fact that at larger radii the agreement is better, which also is the physically more relevant case.

We show in Appendix E that, in the case of purely gravity-driven collections, converges only for very large MBR. Thus, the MBR dependency of the numerical solution using the Smoluchowski scheme appears to be a serious obstacle in studying particle growth not only by condensation, but also by collection. This is a strong argument in favor of the Lagrangian scheme. The evolution of , on the other hand, agrees rather well between the swarm and Eulerian models.

We emphasize that is sensitive to subtle changes in the size distribution, but it is at the same time not really relevant to characterizing the collectional growth toward large particles. As is shown in the following sections, the mean particle radius often increases by not much more than a factor of three (see also the left-hand panel of Figure 2), while the size distribution can become rather broad and its tail can reach the size of raindrops within a relatively short time. In addition to the mean radius, we now also consider size spectra to address the collectional growth to larger particles.

The evolution of size spectra simulated with the Eulerian scheme with 3457 mass bins () is shown as blue lines in Figure 3, while the corresponding size spectra obtained with the swarm model (collection scheme II) with 32 particles per grid point are shown as red curves. The agreement between the Eulerian and Lagrangian schemes is good at early times (), but at late times () the size spectra from the Eulerian approach is broader for the largest sizes (). Shima et al. (2009) found that the results of the super-droplet method (collection scheme II) agree fairly well with the numerical solution of a binned spectral method. It is interesting to note that the size spectra simulated with the swarm model (scheme II) converge to those obtained with the Eulerian model with increasing . This can simply be explained by the fact that more swarms contribute as potential collectional partners and thus ensure more reliable statistics, which was also shown in the work of Shima et al. (2009).

### 3.3 Characterizing the size spectra

#### 3.3.1 Estimating the extent of the size distribution

The size spectra obtained in our simulations are rather broad and cover nearly three orders of magnitude in radius (and six in mass). Even at , the mean radius, , has barely reached (see Figure 2) and does not give any indication about the width of the distribution. The values of the higher moments and lie still only in the middle of the size range; see Figure 3, where we have marked the values of , , , , and by arrows for . We see that the value characterizes rather well the maximum radius of the distribution. The actual maximum is at , but this value is rather noisy, because it represents only a single data point in our simulated volume. We have seen that our smallest and largest moments, and , conveniently bracket the extent of the size distribution. However, for the type of size spectra presented here, the higher moments do not contain any new or independent information, because all moments grow exponentially at nearly the same rate; see Figure 4. This is quantified by the instantaneous growth rates, , which are found to be around for large moments . Thus, at least in the range , the value of is always approximately ten times larger than . However, this ratio can be different in different cases. Knowing therefore the values of and gives a fairly reliable indication about the full extent of the size spectra.

Given that the different moments are not independent, it should not be too surprising that useful information can already be extracted from the first three moments, which is at the heart of the so-called three-moment bulk scheme (Naumann and Seifert, 2016b). This will be discussed next.

#### 3.3.2 Description in terms of the gamma distribution

In the meteorological context, size spectra are often fitted to a gamma distribution (Berry and Reinhardt, 1974; Geoffroy et al., 2010),

(43) |

Here, the factor , with being the gamma function, is included so is normalized such that . In addition to , it has and as independent parameters, which are all functions of time. These are the basic parameters of the three-moment bulk scheme (Naumann and Seifert, 2016b). An advantage of using the gamma distribution lies in the fact that all moments can be calculated analytically. Thus, one may ask for which values of , and do the moments of the gamma distribution agree with those obtained here. For given values of and , we have

(44) |

As stated above, the value of is given by the zeroth moment. In Table LABEL:TamomText we give the values of together with selected normalized moments as well as the resulting parameters and . We also give , which has units of length and can be compared with the normalized moments .

Note that the value of decreases with time and approaches . At the same time, increases and reaches at the last time. Although has units of length and tends to give an indication about the cutoff of the distribution, its value is still five times smaller than that of and thus far away from the maximum droplet radius. This reinforces us in regarding as a useful measure of the largest droplets.

It turns out that the resulting profiles of the gamma distribution capture the broad tail of the distribution remarkably. This is shown in Figure 5, where we compare with the actual size spectra. It is obvious that the actual size distribution shows additional bumps, notably at . This is completely missed when the parameters of the gamma distribution are computed from the full data set. Note also that we have employed a double-logarithmic representation in Figure 5. To model the bump at , one could again use the gamma distribution, but now with moments that are based on in a restricted size range, , for example. This type of approach has been used by Seifert and Beheng (2001), where they describe the size spectra of cloud droplets and raindrops with different distributions.

Naumann and Seifert (2016b) also found reasonable agreement between numerical obtained size spectra and the corresponding gamma distribution. They used the third and sixth moments of the distribution, but in that case the corresponding expressions for and are more complicated. We show in Appendix F that the results do not change much when using and to compute the parameters of the gamma distribution.

Although the agreement between simulated size spectra and the gamma distribution turns out to be reasonable, it only works for . For smaller values of , the function is no longer normalizable, i.e., the integral over diverges for when . Most size spectra observed in meteorology have positive values of (Naumann and Seifert, 2016b). This is mainly because of the effect of evaporation, which is here neglected. Evaporation would lead to a depletion of for small values of and could lead to spectra that are more typical of a gamma distribution. In our case, we have an approximate fall-off over nearly two orders of magnitude, followed by an exponential cutoff. This is why the obtained from the moments turns out to be so close to . Thus, although the usefulness of the gamma distribution is still being debated (Khain et al., 2015), it can actually be remarkably good provided one allows for small negative values of .

### 3.4 Inhomogeneous collection in a straining flow

Spatial variation in the flow leads to local concentrations and thus to large peak values of that shorten the collection time (Saffman and Turner, 1956). Before studying the turbulent case, we consider first collectional growth in a steady two-dimensional (2-D) divergence-free straining flow. The straining flow is numerically inexpensive and easy to control and analyze compared with turbulence.

#### 3.4.1 Pure collection

We consider first the case of pure collection.
In Figure 6 we show the
time evolution of for the swarm model with collection
scheme II at different grid
resolutions ranging from to meshpoints.
Surprisingly, grows more slowly as we increase
the mesh resolution of the swarm model.
Given that the swarm models seem to converge toward the Eulerian
model, we are confronted with the question of what causes the growth
of in the swarm model to slow down at higher mesh resolution.
In this connection, we must emphasize that by default we use
the CIC algorithm to evaluate the gas properties
at the position of each Lagrangian particle.
As explained in section 2.6, the position of the shepherd particle
has no purpose other than to determine in which grid cell the swarm resides.
It is therefore *n*ot better to use any kind
of interpolation in determining the value of the fluid variables at the
position of the swarm, but rather to use NGP mapping.
This will play an important role, as will be discussed now.
For the sake of solving equations (7) and
(8), the use of the CIC algorithm is perfectly valid,
but this would only be relevant for a direct Lagrangian tracking algorithm.
This can be understood by realizing that in the special case of particles
with vanishingly small inertia, the particles will follow their local
fluid cell, and hence, two particles will in the real world never collide.
However, if the CIC scheme is used for equations (7) and
(8), two swarms residing at different positions within the
same grid cell may have different velocities, and hence,
equation (13) may yield a collection.

Since the swarms are filling the entire volume of the grid cell, this means that the two swarms will have different velocities and exist in the same volume, and hence, the swarms may collide. The larger grid cells yield potentially larger velocity differences between the particles, which explains why the collectional growth is larger for the coarser resolutions. When NGP mapping is adopted, the artificial speedup disappears, as shown in the inset of Figure 6.

However, the discrepancy between Lagrangian and Eulerian particle descriptions is still strong for collectional growth in the straining flow as shown in the inset of Figure 6. This is because that in a steady flow, the particles will end up near the vertices of converging flow vectors and will therefore be much more concentrated in the swarm model than what is possible to represent in the Eulerian model. This is evident by comparing the distribution of superparticles belonging to a certain radius (here ) with the corresponding distribution function in the Eulerian model; see Figure 7.

#### 3.4.2 Combined condensation and collection

When both condensation and collection play a role, it is no longer
possible to define a unique timescale, and the solution depends on
both and .
We consider here the straining flow using ,
and ,
which yields .
We investigate the role that particle viscosity and Brownian diffusion
play in simulations using the Eulerian model.
The Brownian motion of the particles is usually small, so the
particle diffusion coefficient in equation (26) should be finite, but small.
Since it is assumed that the particle flows are relatively dilute, there
should be very little interaction between the different particle fluids,
except for the occasional collections. This implies that the
particle viscosity in equation (36) should be close to zero^{3}^{3}3Note that the
particle viscosity represents the coupling between the particle fluids –
not the drag coupling between the particles and the gas phase..
For the Smoluchowski approach, both
and have to be made large in order to stabilize
the simulations in spatially extended cases.
It turns out that the values of these diffusion coefficients
have a surprisingly strong effect on the solutions,
which is shown in Figure 8.
This could be due to the fact that the viscosity between the particle
fluids diffuses the momentum of the particles and thereby modifies
the collection rate.

Comparing now with the swarm approach, which avoids artificial viscosity and enhanced Brownian diffusion altogether, we see from Figure 8 that Eulerian and Lagrangian approaches agree with each other at early times (). After , both swarm and the Eulerian models follow the same trend in the sense that the evolution of shows a bump. The bump occurs earlier for the swarm model than the Eulerian model. In the extreme case that the artificial viscosity in the Eulerian model were zero, the evolution of , as obtained from the swarm model, may come closer that of the Eulerian model. However, owing to the absence of a pressure term for particles, discontinuities would develop in the Eulerian model that destabilize the code if the viscosity and Brownian diffusion are too small. Again, this may be a strong argument in favor of using the swarm model for studying the collectional growth of cloud droplets.

To relate the speed of evolution in Figure 8 to , we plot in the inset of panel (a) the inverse of its unscaled value, , as a function of time. On average, we have . It is comparable to and both are long compared with . The relevant quantity is the scaled value, , which is much larger . This may suggest that the speed of growth is not governed by the spatially averaged kernel, but by its value weighted toward regions where the concentration is high.

We recall that growth of cloud droplets driven by pure collections in the straining flow depends on the models (Eulerian and Lagrangian models; see detailed comparisons in section 3.4.1). This suggests that condensation has a “regularizing” effect in that it makes the overall evolution of much less dependent on the initial conditions and other model details. This is due to the fact that the condensation process with constant positive supersaturation value leads to narrow size spectra of cloud droplets.

Another interesting aspect is the bump in the evolution of the mean radius. At first glance it seems counterintuitive that can actually decrease during some time interval. In Appendix G we consider an example of four particles, two large ones and two small ones. If two small ones collide, we still have the two large ones, but only 3 particles in total after the collection, so the average radius increases from 1/2 to 2/3. On the other hand, if two large ones collide, we are still left with the two small ones and one particle whose radius has only grown by a factor of (the radius scales with the mass to the 1/3 power). The average radius is therefore , which is less than the original mean radius, which is half the radius of the large ones.

### 3.5 Growth of droplets in 2-D turbulence

Turbulence is generally believed to help bridging the size gaps in both cloud droplet and planetesimal formation. In this section, pure turbulence-generated collections are simulated using both the Eulerian and Lagrangian models. We consider a 2-D squared domain of side length at a resolution of meshpoints, with viscosity (which is about 50 times the physical value for air), average forcing wavenumber , i.e., , and a root-mean-square velocity , resulting in a Reynolds number of . Our choice of corresponds to forcing at large scales that are not yet too large to be affected by constraints resulting form the Cartesian geometry. The rate of energy dissipation per unit volume is and the turnover time is . For the Lagrangian model, we use NGP mapping while for the Eulerian model we adopt artificial viscosity and enhanced Brownian diffusivity for the particles ().

#### 3.5.1 Size spectra

Figure 9 shows the comparison of size spectra for the swarm and Eulerian models in 2-D turbulence. The agreement of the spectra as well as the high moments for both schemes is good. Except for the latest time, the corresponding gamma distribution agrees reasonably well with the simulated size spectra, as shown by the black dashed curves in Figure 9; see selected moments and the parameters for the corresponding gamma distribution in Table LABEL:TamomTurbText. Here we also give the parameters and of the corresponding gamma distribution. Again, becomes negative at the last time, but, unlike the case with pure gravity (Table LABEL:TamomText), its value is still far away from . Furthermore, is now 10 times smaller than and does therefore not represent the maximum droplet radius.

We recall that good agreement is also observed in the case with gravity. Hence, we conclude that the gamma distribution can reasonably well represent the collectional size spectra. This may provide a strong argument in favor of using the gamma distribution when modelling cloud microphysics.

#### 3.5.2 Other numerical aspects

It is worth noting that the MBR convergence of the Smoluchowski equation depends on the flow pattern. While gravitational collection is rather sensitive to MBR (see Appendix E), it is much less sensitive for the straining flow and converges at in turbulence.

We emphasize that for the swarm model, the interpolation scheme of the tracked swarms does affect the results, but this does not seem to be the case for turbulence. Turbulence continues to mix particles all the time while the straining flow tends to sweep up particles into predetermined locations that do not change. We may therefore conclude that the restriction on the interpolation scheme depends on the spatio-temporal properties of the flow. Nevertheless, a high-order interpolation is not strictly applicable to the swarm model.

It is worth noting that in the case with pure gravity, the Eulerian model is rather sensitive to the presence or absence of the term. This is neither the case for turbulence nor for the straining flow as will be discussed in Appendix B.

#### 3.5.3 Comparison of computational cost

Comparing the Lagrangian and Eulerian models in Figure 9, it is worth noting that the Lagrangian one is clearly superior to the Eulerian one in terms of CPU time for simulating the collectional process in 2-D turbulence. A similar conclusion was drawn by Shima et al. (2009), who found the Lagrangian model to be computationally faster than the Eulerian one. We compare the computational cost between Eulerian and Lagrangian models using the 2-D turbulence runs 10A and 10B (runs in Figure 9), which have comparable accuracy; see Table 1 for details of these runs. The Lagrangian model with superparticles covers s in physical time within 24 hours of wall-clock time on 512 CPUs, while the Eulerian model with 53 mass bins covers only s in physical time within 24 hours wall-clock time on 1024 CPUs. This example demonstrates that the Lagrangian model is roughly ten times more efficient than a comparable Eulerian one.

#### 3.5.4 Combined condensation and collection

The combined condensational and collectional growth in turbulence is investigated as well. Again, the results are similar to the case with pure collectional growth due to the fact that the condensation process in the present study with constant supersaturation is homogeneous. In future studies, the supersaturation should be calculated self-consistently and the effects of turbulence on the condensational growth should be considered, similar to what was done previously (Kumar et al., 2014; Sardina et al., 2015).

## 4 Conclusion

The combined collectional and condensational growth of cloud droplets is studied in numerical simulations where the gas phase is solved on a mesh, while the particle phase is approximated by a point particle approach and is treated either by an Eulerian or a Lagrangian formalism. It is found that the Lagrangian approach agrees well with the analytic solution of condensational growth. By contrast, the Eulerian approach requires high resolution in the number of mass bins to avoid artificial speedup of the growth rate, which agrees with previous findings (Ohtsuki et al., 1990; Dra̧żkowska et al., 2014). It is worth noting that the MBR dependency is closely related to the temporal and spatial properties of the flow. The dependency is strongest for gravity [ in equation (1)], less strongly for the straining flow, and weak for turbulence.

A detailed comparison of the collectional size spectra between the Lagrangian and Eulerian models demonstrates consistency between the two, especially when both condensation and collection are included. This suggests that condensation has a regularizing effect and makes the overall evolution of the mean radius less dependent on details such as the precise form of the initial condition or discretization errors that might affect the early evolution. However, the evolution of the mean radius, i.e., the ratio of the two lowest (first and zeroth) moments of the size distribution function, is a rather insensitive measure of particle growth. This is also seen in the fact that the mean particle radius often increases by not much more than a factor of three, while the size distribution can become rather broad and even millimeter-sized particles can be produced within a relatively short time. The mean particle radius is also not the most relevant diagnostics in that it does not characterize properly the growth of the largest particles. In fact, as we have shown in Appendix G, the mean radius actually decreases when two large particles collide. This is somewhat counterintuitive, but actually quite natural. When two very small particles collide, the sum of all radii does basically not change, but the number of particles decreased by one, so the average increases. By contrast, when two large particles collide, the particle number again decreases by one, but the sum of the radii decreases from 2 to , so the average also decreases.

Remarkably, we found that the simulated tails of the size spectra agree fairly well with those obtained from the gamma distribution for negative . More importantly, the agreement is good both for gravity and turbulence cases.

When studying pure collection, the Eulerian approach yields satisfactory results only when the mass bins are sufficiently fine. Furthermore, for collections in the case of a straining flow, it is found that the Eulerian approach requires artificially large viscosity and Brownian diffusivity for keeping the resulting shocks in the particle fluid resolved. Because of this, it seems that for future studies of the effect of turbulence on condensational and collectional growth of particles, the Lagrangian swarm approach would be most suitable. However, several precautions have to be taken. First, the symmetric collection scheme II (Shima et al., 2009) is to be preferred because it shows less scatter in the mean radius than the asymmetric scheme I. This is because in scheme I the particle number is adjusted to keep the total mass in the swarm constant. Second, when interpolation of the gas properties at the position of each Lagrangian particle is invoked (for example the CIC algorithm or the triangular shaped cloud scheme), both collection schemes yield artificially increased collection rates. This is because two swarms within the same grid cell may always collide since the interpolation of the fluid velocity results in a velocity difference between the two swarms. This causes a speedup of the collection rate already at early times. At higher grid resolution, the interpolated velocity differences are smaller, which reduces the collectional growth. Therefore, it is best to map the gas properties to just the nearest grid point, which is found to yield converged results even at low resolution.

A shortcoming of the Eulerian model is that self-collections are strictly impossible. This should be mitigated by using finer mass bins, but it turns out that finer mass bins do not change the collection rate at early times, but rather decrease it at later times. This indicates that the contribution of self-collections to the collection rate is relatively small; see Appendix H for details.

The discrepancy between Lagrangian and Eulerian particle descriptions is particularly strong in the time-independent straining flow. This is because particles tend to be swept into extremely narrow lanes, which leads to high concentrations that can never be achieved with the Eulerian approach, in which sharp gradients must be smeared out by invoking artificial viscosity and large Brownian diffusivity. On the other hand, we are here primarily interested in turbulent flows that are always time-dependent, which limits the amount of particle concentration that can be achieved in a given time. In that case, the discrepancies between Eulerian and Lagrangian approaches are smaller at early times, but there are still differences in the evolution of the mean radius at late times. This can easily be caused by changes in the relative importance of collections of large and small particles. This is confirmed by the fact that the size distribution spectra in the turbulent case are more similar for Lagrangian and Eulerian approaches than in the straining flow.

Our present work neglects local and temporal changes in the supersaturation. In future studies, we will take into account that the supersaturation increases (decreases) as a fluid parcel rises (falls) and that droplet condensation (evaporation) act as sinks (sources) of supersaturation. We would then be able to account for the fact that the total water content should remain constant and that the supersaturation would become progressively more limited as water droplets grow by condensation. Another important shortcoming is our assumption of perfect collection efficiency, which resulted in artificially rapid cloud droplets growth. Alleviating these restrictions will be important tasks for future work. Furthermore, we have here only considered 2-D turbulence. Extending our work to 3-D is straightforward, but our conclusions regarding the comparison of different schemes should carry over to 3-D.

###### Acknowledgements.

We thank Nathan Kleeorin, Dhrubaditya Mitra, and Igor Rogachevskii for useful discussions. We also thank the anonymous referees for constructive comments and suggestions that lead to substantial improvements in the manuscript. This work was supported through the FRINATEK grant 231444 under the Research Council of Norway, the Swedish Research Council grant 2012-5797, and the grant “Bottlenecks for particle growth in turbulent aerosols” from the Knut and Alice Wallenberg Foundation, Dnr. KAW 2014.0048. The simulations were performed using resources provided by the Swedish National Infrastructure for Computing (SNIC) at the Royal Institute of Technology in Stockholm. This work utilized the Janus supercomputer, which is supported by the National Science Foundation (award number CNS-0821794), the University of Colorado Boulder, the University of Colorado Denver, and the National Center for Atmospheric Research. The Janus supercomputer is operated by the University of Colorado Boulder. The source code used for the simulations of this study, the Pencil Code, is freely available on https://github.com/pencil-code/. The input files as well as some of the output files of the simulations listed in Table 1 are available under http://www.nordita.org/~brandenb/projects/SwarmSmolu_numerics/.## Appendix A Upwinding scheme for a nonuniform mesh

In the presence of condensation alone, the evolution equation for as a function of radius and time is given by

(A.1) |

where and is given by equation (6). Thus, we have

(A.2) |

where is assumed independent of ; see equation (13.14) of Seinfeld and Pandis (2006). It can be seen from the form of the analytic solution that there will be a discontinuity at , which is numerically difficult to handle. In particular, it is difficult to ensure the positivity of . For these reasons, a low-order upwind scheme is advantageous. Furthermore, expanding the RHS of equation (A.2) using the quotient rule,

(A.3) |

it is obvious that the first term in isolation would lead to exponential growth of proportional to , which must be partially canceled by the second term. If the cancellation is numerically imperfect, will indeed grow exponentially, which tends to occur in regions where , i.e., where should vanish. For nonuniform mesh spacing, with , 2, …, , the first-order upwind scheme can be written as