Explicit stochastic advection algorithms for the regional-scale particle-resolved atmospheric aerosol model WRF-PartMC (v1.0)

This paper presents the development of a stochastic particle method to simulate advection in regional-scale models with a particle-resolving aerosol representation. The new method is based on finite-volume discretizations with the flux terms interpreted as probabilities of particle transport between grid cells. We analyze the method in 1D and show that the stochastic particle sampling during transport injects energy at high spatial frequencies, which can be partially compensated for with the choice of a dissipative odd-order finite-volume scheme. We then apply the stochastic third- and fifth-order advection algorithms with monotonic limiters in WRF-PartMC, using idealized and realistic wind fields in 2D and 3D. In all cases we observe the expected convergence rates of the stochastic particle method to the finite-volume solution as the number of computational particles is increased. This work enables the use of particle-based aerosol models on the regional scale.

1 Introduction

Aerosol particles influence the climate system as cloud condensation nuclei (CCN), as ice nucleating particles, and as scatterers and absorbers of radiation (Masson-Delmotte et al., 2021). Estimating the magnitude of the aerosol impact on climate requires not only the information of bulk aerosol composition and size distribution, but also the information of the aerosol mixing state (Riemer et al., 2019), i.e., the way the chemical species are distributed across the particle population (Winkler, 1973). The aerosol mixing state can vary between a fully external mixture, where each particle contains only one chemical species which can differ between different particles, and a fully internal mixture, where each particle is composed of the same mixture of species. In reality, the mixing state is in between these two extreme cases (Bondy et al., 2018; O'Brien et al., 2015; Ye et al., 2018; Healy et al., 2014). Furthermore, many physical and chemical processes change the mixing state during the aerosol's lifetime in the atmosphere (Li et al., 2016). Representing these processes in models poses large challenges but is needed to predict the aerosol climate impact (Bauer et al., 2013; Fierce et al., 2017).

Atmospheric three-dimensional chemical transport models or Earth system models utilize a variety of aerosol representations that differ in their levels of detail. These can be categorized into bulk approaches (Koch, 2001; Tegen and Miller, 1998), modal modeling approaches (Whitby and McMurry, 1997), and sectional modeling approaches (Seigneur et al., 1986). These methods have in common that they do not fully resolve the mixing state of the aerosol but instead use a priori assumptions. For example, modal models assume that each mode is internally mixed, while different modes can differ in the set of species that they track. Sectional models capture the size dependence of aerosol composition, but within one section only the average aerosol composition is known. These approaches can be refined by introducing additional modes (Bauer et al., 2008; Liu et al., 2012, 2016), additional one-dimensional sectional distributions (Jacobson, 2002; Zhang et al., 2014), or additional dimensions to the bin structure itself (Matsui et al., 2013; Matsui, 2016; Zhu et al., 2015; Ching et al., 2016), where each dimension represents one species or group of species. Comparing these more sophisticated types of models against versions that use more simplified mixing-state representations shows that mixing-state approximations impact the estimation of optical and CCN properties and contribute to the structural and parametric model uncertainties. For example, Zhu et al. (2016) performed simulations with a sophisticated mixing-state-aware model (SCRAMS) for the region of Paris, France. Different mixing-state treatments caused differences in aerosol water uptake, which propagated into differences in aerosol optical depths of up to 70 %. Lee et al. (2016) carried out simulations with a mixing-state-resolving (source-oriented) version of WRF-Chem for the region of the Californian Central Valley. They found a decrease in the ratio of CCN to total aerosol number concentration from 94 % with an internal mixture assumption to 80 % with a more detailed source-oriented mixture. Furthermore, the range of uncertainties can depend on the degree to which mixing state is represented. This was shown by Matsui et al. (2018), who quantified the sensitivity of the present-day BC direct radiative effect due to uncertainties in emission size distributions. They found that the uncertainty is 5–7 times larger when the BC mixing state is sufficiently resolved compared to a simplified model representation where an internal mixture is assumed.

It is important to note that the storage requirements for multi-dimensional bin structures grow exponentially with the number of species (the curse of dimensionality). Therefore, in practice, the multi-dimensional bin approach is limited to two or three dimensions, whereas the composition space of the atmospheric aerosol contains tens or even hundreds of species. Hence, although this model approach carries more detail than 1D bin structures, it is still not able to resolve the mixing state fully.

In contrast to the above-mentioned distribution-based methods, particle-resolved methods provide a different approach to representing the atmospheric aerosol (Riemer et al., 2009; Shima et al., 2009; Grabowski et al., 2019). They use a collection of discrete computational particles, where each particle can be thought of as a vector that stores the masses of each aerosol species and other particle attributes (e.g., information about particle shape or particle source) and that evolves over the course of the simulation. Aerosol mixing state is therefore intrinsically resolved and does not require any ad hoc assumptions. Furthermore, it is straightforward to add more attributes to the particles as this does not result in an exponential increase in storage. Instead, it scales linearly with the number of particles. Particle methods are therefore beneficial for problems where high-dimensional data are involved as they break the curse of dimensionality.

In this paper we describe the development of stochastic advection algorithms that enable the particle-resolved aerosol model PartMC to be used on the regional scale, embedded within the Weather Research and Forecast (WRF) model. While we only present the development of stochastic advection schemes based on the finite-volume methods in WRF, the methodology described here is applicable to any finite-volume scheme or transport scheme such as corner-transport upwind (Colella, 1990; LeVeque, 2002) or flux-form semi-Lagrangian (Lin and Rood, 1996, 1997) that can be found in other host models. This paper builds on previous work of developing the stochastic, particle-resolved PartMC-MOSAIC box model (Riemer et al., 2009; DeVille et al., 2011; Curtis et al., 2016; DeVille et al., 2019) and the one-dimensional single-column model WRF-PartMC-MOSAIC-SCM (Curtis et al., 2017). These modeling tools have been used to investigate the black carbon aging process (Riemer et al., 2010; Fierce et al., 2015), to quantify the role of the mixing state in determining CCN concentration (Ching et al., 2012, 2016, 2017) and aerosol optical properties (Fierce et al., 2016; Yao et al., 2022), and to determine structural uncertainty in more approximate aerosol models (Fierce et al., 2017; Zheng et al., 2021).

The methods for particle transport due to turbulent diffusion described in Curtis et al. (2017) and for mean wind advection described in this paper are based on the idea that the movement of particles between grid cells is represented by stochastic sampling. Importantly, the particle position within the grid cell is not tracked. This concept is therefore distinct from the particle-based Lagrangian techniques in the cloud modeling community (Heus et al., 2010; Arabas et al., 2015; Grabowski et al., 2018). These methods are similar to ours in that they also explicitly simulate microphysical processes on a population of computational particles (called super-particles in the cloud physics community), each representing a large number of real particles. However, they are different in that they simulate transport by tracking the super-droplet positions within the Eulerian grid.

There are advantages and disadvantages to each method. First, a stochastic algorithm can be constructed analogously to the finite-volume transport schemes used in numerical weather models and chemical transport models, as we will show in this paper. This is beneficial for direct comparisons of different aerosol representations, which is one of our main motivations for developing particle-resolved aerosol models on the regional scale. Second, stochastic methods are more easily implemented in models that rely on different numerical grid structures, because they are based on the discretizations of the host model on the host grid. Lastly, stochastic methods for transport are computationally less expensive than tracking and updating particle positions throughout the simulation. However, stochastic transport algorithms have the disadvantage of numerical diffusion, similar to finite-volume methods. This is in contrast to Lagrangian particle tracking methods that are inherently free of numerical diffusion.

The contribution of this paper is the development of stochastic transport algorithms for advection that remove the modeling limitations of the single-column model (WRF-PartMC-MOSAIC-SCM) by enabling a fully three-dimensional model to allow particle-resolved simulations on the regional scale (WRF-PartMC). WRF-PartMC is a tool for error quantification and benchmarking of traditional chemistry-transport models (e.g., WRF-Chem or CMAQ) that apply simplified aerosol-mixing-state representation, without the advection schemes being a potential source of differences.

The paper is structured as follows. Section 2 develops the stochastic particle advection method. Section 3 analyzes a series of four numerical experiments of increasing complexity, ranging from simple one-dimensional test cases with constant, uniform wind fields to a simulation with complex terrain and evolving meteorological fields. Section 4 summarizes our work. See Table G1 for a list of symbols used throughout the paper.

2 Stochastic particle transport scheme

This section describes the spatial and temporal discretization of the advection equation and then explains the stochastic sampling algorithm for the use in particle-resolved models. We will present the detailed derivation for one spatial dimension. The generalization to three dimensions in space is straightforward and for brevity will not be explicitly written out, although see Sect. 2.3 for notes on implementation details. In this study, we adopted the advection methods implemented by the host model WRF, rather than exploring alternative approaches. This choice ensures that future comparisons between WRF-PartMC and other aerosol representations in WRF-Chem will be fair and consistent. The WRF-PartMC model was developed using WRFv3.9.1.1, coupled with chemistry for gas scalar transport and PartMC-MOSAIC for gas and aerosol chemistry, with an additional interface for simulating stochastic particle transport.

2.1 Spatial and temporal discretization

The one-dimensional advection equation of a scalar quantity with (number) concentration n(x,t) can be written as

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial n(x,t)}{\partial t}}=-u{\displaystyle \frac{\partial n(x,t)}{\partial x}},\end{array}$$

where u>0 is the velocity of the advecting wind field (assumed to be constant in time and uniform in space here), x is the spatial coordinate, and t is time.

We discretize this equation spatially as

$$\begin{array}{}\text{(2)}& {\displaystyle \frac{\partial n_i\left(t\right)}{\partial t}}=-{\displaystyle \frac{\mathrm{1}}{\mathrm{\Delta }x}}\left(f_{i+\frac{\mathrm{1}}{\mathrm{2}}}\right(t)-f_{i-\frac{\mathrm{1}}{\mathrm{2}}}(t\left)\right),\end{array}$$

where Δx is the grid spacing in the x coordinate, and f_i+½(t) and f_i−½(t) are the fluxes through the right and left grid cell boundaries of grid cell i at time t, respectively, for $i=\mathrm{0},\mathrm{\dots },N_{\mathrm{x}}-\mathrm{1}$.

The fluxes can be spatially discretized to different orders, with the WRF schemes of orders 1 to 6 written as (Wicker and Skamarock, 2002; Shu, 2009)

$$\begin{array}{ll}\text{(3)}& {\displaystyle }{f}_{i-\frac{\mathrm{1}}{\mathrm{2}}}^{\text{1st}}& {\displaystyle }=u{n}_{i-\mathrm{1}},\text{(4)}& {\displaystyle }{f}_{i-\frac{\mathrm{1}}{\mathrm{2}}}^{\text{2nd}}& {\displaystyle }={\displaystyle \frac{u}{\mathrm{2}}}(n_i+n_{i-\mathrm{1}}),\text{(5)}& {\displaystyle }{f}_{i-\frac{\mathrm{1}}{\mathrm{2}}}^{\text{3rd}}& {\displaystyle }={\displaystyle \frac{u}{\mathrm{6}}}(\mathrm{2}{n}_i+\mathrm{5}{n}_{i-\mathrm{1}}-n_{i-\mathrm{2}}),\text{(6)}& {\displaystyle }{f}_{i-\frac{\mathrm{1}}{\mathrm{2}}}^{\text{4th}}& {\displaystyle }={\displaystyle \frac{u}{\mathrm{12}}}(-n_{i+\mathrm{1}}+\mathrm{7}{n}_i+\mathrm{7}{n}_{i-\mathrm{1}}-n_{i-\mathrm{2}}),{\displaystyle }{f}_{i-\frac{\mathrm{1}}{\mathrm{2}}}^{\text{5th}}& {\displaystyle }=\\ \text{(7)}& {\displaystyle }& {\displaystyle \frac{u}{\mathrm{60}}}(-\mathrm{3}{n}_{i+\mathrm{1}}+\mathrm{27}{n}_i+\mathrm{47}{n}_{i-\mathrm{1}}-\mathrm{13}{n}_{i-\mathrm{2}}+\mathrm{2}{n}_{i-\mathrm{3}}),{\displaystyle }{f}_{i-\frac{\mathrm{1}}{\mathrm{2}}}^{\text{6th}}& {\displaystyle }=\\ \text{(8)}& {\displaystyle }& {\displaystyle \frac{u}{\mathrm{60}}}(n_{i+\mathrm{2}}-\mathrm{8}{n}_{i+\mathrm{1}}+\mathrm{37}{n}_i+\mathrm{37}{n}_{i-\mathrm{1}}-\mathrm{8}{n}_{i-\mathrm{2}}+n_{i-\mathrm{3}}),\end{array}$$

and similarly for the fluxes through the other boundary, f_i+½. We will explore the effect of using different orders of discretization in the context of stochastic particle-based advection in Sect. 3.1 and 3.2.

For the temporal discretization, we use a third-order Runge–Kutta method, analogous to the approach in WRF (Wicker and Skamarock, 2002), where the concentration at time ℓ+1 is calculated from the values at time ℓ as

$$\begin{array}{}\text{(9)}& {\displaystyle }{n}_i^{†}& {\displaystyle }=n_i^{\mathrm{\ell }}-{\displaystyle \frac{\mathrm{\Delta }t}{\mathrm{3}}}{\displaystyle \frac{\mathrm{1}}{\mathrm{\Delta }x}}(f_{i+\frac{\mathrm{1}}{\mathrm{2}}}^{\mathrm{\ell }}-f_{i-\frac{\mathrm{1}}{\mathrm{2}}}^{\mathrm{\ell }}),\text{(10)}& {\displaystyle }{n}_i^{††}& {\displaystyle }=n_i^{\mathrm{\ell }}-{\displaystyle \frac{\mathrm{\Delta }t}{\mathrm{2}}}{\displaystyle \frac{\mathrm{1}}{\mathrm{\Delta }x}}(f_{i+\frac{\mathrm{1}}{\mathrm{2}}}^{†}-f_{i-\frac{\mathrm{1}}{\mathrm{2}}}^{†}),\text{(11)}& {\displaystyle }{n}_i^{\mathrm{\ell }+\mathrm{1}}& {\displaystyle }=n_i^{\mathrm{\ell }}-\mathrm{\Delta }t{\displaystyle \frac{\mathrm{1}}{\mathrm{\Delta }x}}(f_{i+\frac{\mathrm{1}}{\mathrm{2}}}^{††}-f_{i-\frac{\mathrm{1}}{\mathrm{2}}}^{††}),\end{array}$$

where Δt is the time step. The fluxes $f_{i+\frac{\mathrm{1}}{\mathrm{2}}}^{†}$, $f_{i-\frac{\mathrm{1}}{\mathrm{2}}}^{†}$, $f_{i+\frac{\mathrm{1}}{\mathrm{2}}}^{††}$, and $f_{i-\frac{\mathrm{1}}{\mathrm{2}}}^{††}$ are calculated as given in Eqs. (3)–(8), using the concentrations $n_i^{†}$ and $n_i^{††}$, respectively.

Next, we will show how the discretized equations defined above are translated to specify the probabilities of particles moving between grid cells. Although we are presenting the method using the particular discretizations above, it is straightforward to derive stochastic versions of other spatial and temporal discretizations in the same way.

2.2 Stochastic sampling

To transform the method from Sect. 2.1 to a stochastic particle method, we consider a set of $N_i^{\mathrm{\ell }}$ particles in grid cell i at time step ℓ. In reality, each particle will have an exact spatial location with a well-defined x coordinate and will be moving with constant velocity u. However, in our stochastic method we will not track this per-particle spatial location and instead only track the set of particles in each grid cell. This is equivalent to the usual finite-volume method of tracking the concentration in each grid cell, except that we are now sampling the concentration with a finite set of particles, allowing us to capture the high-dimensional variation in particle properties.

Note that the full WRF-PartMC model implementation explicitly tracks each particle in each grid cell so that it can store additional information about each particle (particle diameter, chemical constituents, etc.). In the following exposition we will not explicitly track the particles but instead will only track the number of particles, $N_i^{\mathrm{\ell }}$, in each grid cell. Section 2.4 contains further comments on translating the count-based scheme to a true per-particle method.

To advect the particles, we relate the number of particles in each grid cell to the concentration in that grid cell, using

$$\begin{array}{}\text{(12)}& n_i^{\mathrm{\ell }}={\displaystyle \frac{N_i^{\mathrm{\ell }}}{V}},\end{array}$$

where V is the computational sampling volume within each grid cell. We can think of V as controlling the “resolution” of the particle sampling, and it will generally be much smaller than the true grid cell volume.

Having computed the values of $n_i^{\mathrm{\ell }}$ for each grid cell, we then compute the finite-volume fluxes $f_{i+\frac{\mathrm{1}}{\mathrm{2}}}^{††}$ through each boundary from Eq. (11). This tell us that the average number of particles that should cross boundary $i+\frac{\mathrm{1}}{\mathrm{2}}$ is

$$\begin{array}{}\text{(13)}& {\stackrel{\mathrm{‾}}{F}}_{i+\frac{\mathrm{1}}{\mathrm{2}}}^{\mathrm{\ell }}=V{\displaystyle \frac{\mathrm{\Delta }t}{\mathrm{\Delta }x}}{f}_{i+\frac{\mathrm{1}}{\mathrm{2}}}^{††}.\end{array}$$

We interpret this probabilistically to mean that each of the $N_i^{\mathrm{\ell }}$ particles has a probability of

$$\begin{array}{}\text{(14)}& p_{i+\frac{\mathrm{1}}{\mathrm{2}}}^{\mathrm{\ell }}={\displaystyle \frac{{\stackrel{\mathrm{‾}}{F}}_{i+\frac{\mathrm{1}}{\mathrm{2}}}^{\mathrm{\ell }}}{N_i^{\mathrm{\ell }}}}\end{array}$$

of crossing the boundary and leaving grid cell i. We then sample the number of particles that actually cross the boundary using a binomial distribution with $N_i^{\mathrm{\ell }}$ trials and probability $p_{i+\frac{\mathrm{1}}{\mathrm{2}}}^{\mathrm{\ell }}$ to give the discrete particle flux across the boundary to be

$$\begin{array}{}\text{(15)}& F_{i+\frac{\mathrm{1}}{\mathrm{2}}}^{\mathrm{\ell }}=\text{Binom}(N_i^{\mathrm{\ell }},p_{i+\frac{\mathrm{1}}{\mathrm{2}}}^{\mathrm{\ell }}).\end{array}$$

Finally, we update the number of particles in each grid cell according to

$$\begin{array}{}\text{(16)}& N_i^{\mathrm{\ell }+\mathrm{1}}=N_i^{\mathrm{\ell }}-F_{i+\frac{\mathrm{1}}{\mathrm{2}}}^{\mathrm{\ell }}+F_{i-\frac{\mathrm{1}}{\mathrm{2}}}^{\mathrm{\ell }}.\end{array}$$

Note that this method obviously conserves the total number of discrete particles, because the $F_{i+\frac{\mathrm{1}}{\mathrm{2}}}^{\mathrm{\ell }}$ particles that leave grid cell i will all be transferred to grid cell i+1. In addition, because the mean of a binomial distribution is equal to the number of trials times the probability of success, we see that the average value of $N_i^{\mathrm{1}}$ is exactly equal to $V{n}_i^{\mathrm{1}}$ for the first time step. However, because the next time step will start from the stochastically sampled discrete value $N_i^{\mathrm{1}}$, the average value of $N_i^{\mathrm{2}}$ will not be exactly equal to the average value of $V{n}_i^{\mathrm{2}}$.

Since probabilities larger than 1 are not meaningful, the time step needs to be chosen such that the probability (14) is less than or equal to 1. As a result, WRF-PartMC may need to take somewhat smaller time steps than required by the finite-difference advection in WRF.

2.3 Three-dimensional advection

The above derivation is for a one-dimensional domain, but the extension to three dimensions is straightforward. In three dimensions, we have a set of $N_{i,j,k}^{\mathrm{\ell }}$ particles in grid cell $(i,j,k)$ at time step ℓ. The fluxes are then defined as $f_{i+\frac{\mathrm{1}}{\mathrm{2}},j,k}^{\mathrm{\ell }}$, $f_{i,j+\frac{\mathrm{1}}{\mathrm{2}},k}^{\mathrm{\ell }}$, and $f_{i,j,k+\frac{\mathrm{1}}{\mathrm{2}}}^{\mathrm{\ell }}$ for the fluxes through the three positive boundaries of grid cell $(i,j,k)$, respectively. The fluxes through the other boundaries are defined similarly. The fluxes are then computed from a 3D finite-volume discretization, but with the concentrations n_i replaced by the number of particles $N_{i,j,k}^{\mathrm{\ell }}$. This then yields probabilities of particles crossing each boundary by extending Eq. (14) from i to three dimensions $(i,j,k)$. However, we now have three different probabilities, one for each boundary, corresponding to the three different directions in which particles can move. The time step should be chosen so that the sum of these probabilities is at most 1. We then sample the number of particles that move in each direction using a multinomial distribution with $N_{i,j,k}^{\mathrm{\ell }}$ trials and probabilities $p_{i+\frac{\mathrm{1}}{\mathrm{2}},j,k}^{\mathrm{\ell }}$, $p_{i,j+\frac{\mathrm{1}}{\mathrm{2}},k}^{\mathrm{\ell }}$, and $p_{i,j,k+\frac{\mathrm{1}}{\mathrm{2}}}^{\mathrm{\ell }}$ for the three directions, respectively. See Curtis et al. (2017) for a detailed description of the multinomial sampling algorithm. Finally, the number of particles in each grid cell is updated by extending Eq. (16) from one dimension (i) to three dimensions $(i,j,k)$.

2.4 Explicit tracking of individual particles

Much of the power of a particle-based aerosol model is the ability to track the chemical composition and potential morphology of individual particles, as is done by the PartMC (Riemer et al., 2009) model, which explicitly tracks a set of particles ${\mathrm{\Pi }}_{i,j,k}^{\mathrm{\ell }}$ in grid cell $(i,j,k)$ at time step ℓ. To apply the stochastic advection algorithm of Sect. 2.2 to such a case, the stochastic fluxes can be computed using the total number, $N_{i,j,k}^{\mathrm{\ell }}$, of particles in each grid cell to give $F_{i+\frac{\mathrm{1}}{\mathrm{2}},j,k}^{\mathrm{\ell }}$, $F_{i+\frac{\mathrm{1}}{\mathrm{2}},j,k}^{\mathrm{\ell }}$, and $F_{i+\frac{\mathrm{1}}{\mathrm{2}},j,k}^{\mathrm{\ell }}$ as the number of particles that will cross each boundary. However, rather than simply updating the particle counts using the fluxes, we uniformly randomly sample $F_{i+\frac{\mathrm{1}}{\mathrm{2}},j,k}^{\mathrm{\ell }}$ particles from the set ${\mathrm{\Pi }}_{i,j,k}^{\mathrm{\ell }}$ to move across the boundary and similarly for the other two directions. This approach is used in the WRF-PartMC model.

2.5 Monotonicity

It is advisable in WRF-Chem simulations to use monotonic, positive-definite advection schemes (Wang et al., 2009; Chapman et al., 2009). WRF advection schemes without limiters have the tendency to overshoot as well as locally produce unrealistically low values. This is particularly problematic for chemical variables that have strong gradients due to the heterogeneity of emissions. The host WRF model features only a fifth-order scheme for monotonic limiters. However due to the high computational expense of WRF-PartMC and the required domain decomposition to adequately meet that expense, we implemented third-order advection with monotonic limiters in WRF. This implementation utilized the existing third-order positive-definite scheme in WRF and applied the same limiter as used in the fifth-order monotonic scheme term (Skamarock, 2006; Wang et al., 2009). We focused on third- and fifth-order advection schemes because they combine good accuracy with some numerical dissipation at high spatial frequencies to suppress stochastic oscillations, as we will see in Sect. 3.1 and 3.2.

2.6 Mixing ratio versus concentration

Many 3D atmospheric models such as WRF track the aerosol mixing ratio q (units of particles kg⁻¹) rather than the number concentration n (units particles m⁻³) because this removes the need to adjust the tracer for changes in air density. However, the stochastic advection algorithm described above is based on the number concentration. In WRF-PartMC we convert the aerosol number concentration to a mass mixing ratio via $q=n/\mathit{\rho }$, compute mixing-ratio fluxes using WRF's finite-volume discretization, convert these back to number-concentration fluxes by multiplying by ρ, and then sample the stochastic particle transport with Eqs. (14)–(16).

2.7 Variable sampling volumes and grid cell sizes

In Sect. 2.2 we assumed that the sampling volume V is a constant. However, in the WRF-PartMC model the sampling volume is allowed to vary in space and time. This is done by defining a set of $V_{i,j,k}^{\mathrm{\ell }}$ volumes in each grid cell $(i,j,k)$ at time step ℓ, and this allows the “particle resolution” to be adaptive to increase the accuracy while minimizing the computational cost. In such simulations a target number of particles per grid cell, N_p, is chosen to be a fixed value, and the sampling volume is then adapted using a halving/doubling procedure to maintain the actual number of particles per grid cell close to N_p (Riemer et al., 2009).

As described in detail in Curtis et al. (2017), the variable sampling volumes mean that the number of particles that move out of a grid cell (the particle loss) is no longer generally equal to the number of particles that move into the neighboring cell (the particle gain). Instead, if$F_{i+\frac{\mathrm{1}}{\mathrm{2}},j,k}^{\mathrm{\ell }}$ particles move out of grid cell $(i,j,k)$, then the number of particles that move into grid cell $(i+\mathrm{1},j,k)$ is scaled by the ratio of the sampling volumes. That is, the number of particles that move into grid cell $(i+\mathrm{1},j,k)$ is

$$\begin{array}{}\text{(17)}& F_{i+\frac{\mathrm{1}}{\mathrm{2}},j,k}^{\mathrm{\ell }}{\displaystyle \frac{V_{i+\mathrm{1},j,k}^{\mathrm{\ell }}}{V_{i,j,k}^{\mathrm{\ell }}}}.\end{array}$$

Similarly, if the grid cells have different physical volumes ${\text{Vol}}_{i,j,k}^{\mathrm{\ell }}$, then the above expression must be additionally scaled by the ratio ${\text{Vol}}_{i,j,k}^{\mathrm{\ell }}/{\text{Vol}}_{i+\mathrm{1},j,k}^{\mathrm{\ell }}$. The advection algorithm in WRF-PartMC implements this scaling following the method in Curtis et al. (2017), which uses a variance-minimizing sampling algorithm that first samples the larger of the particle loss and gain terms and then subsamples from this to determine the other term. A potential concern is that the repeated resampling due to varying computational volumes, grid cell volumes, and air densities may cause the high-dimensional information carried by particles (see Sect. 2.4) to degenerate into overly similar representations. For example, if the particles carry a diameter sampled from a size distribution, the repeated resampling may cause the particles to converge to a single diameter. In Sect. 3.4 we investigate this numerically and see that it is not a significant issue in practice.

2.8 Computational cost

Regarding the computational costs of the finite-volume, stochastic sampling, and Lagrangian particle tracking approaches, we consider a domain consisting of N_g grid cells and N_p computational particles per grid cell. The finite-volume method, which only depends on the number of grid cells, has a cost 𝒪(N_g). In contrast, the Lagrangian particle tracking and stochastic methods depend on both the number of grid cells and the number of particles. Therefore, these methods scale as 𝒪(N_g×N_p), but the Lagrangian method has a higher cost as each particle must be checked and updated. In contrast, the cost of the stochastic method depends on the number of particles that actually move from one grid cell to another, which is frequently only a small fraction of the total number.

Figure 1 One-dimensional soft-hat test case (Sect. 3.1): a single ensemble member of the stochastic solution is shown in blue for first- to sixth-order methods, the deterministic finite-volume solution is represented by the solid red line, and the analytical solution is shown as a dashed black line. The stochastic solution was simulated using N_p=10⁴ computational particles per grid cell.

2.9 Comparison to Lagrangian particle tracking

With particles transported by deterministic advection there is no variance in the final position of particles that start in the same initial position. However, when we quantize space and only store which grid cell a particle is in, we can no longer move particles to the exact position where they should be located. That is, we are forced to incur some error. In a classic bias–variance tradeoff, we could achieve zero variance by moving all collocated particles to the same new grid cell, but this would result in an incorrect average position of the particles and a large bias. Alternatively, as we do in this paper, we can move some particles and not others, resulting in the correct mean velocity (zero bias) at the expense of introducing variance in particle position. Consequently, some particles will move faster and some slower than the mean velocity. To quantify the magnitude of this effect, see the example in Sect. 3.2 and Fig. 7.

3 Numerical experiments

The total error of a stochastic transport scheme can be bounded by two error terms that can be evaluated independently: (1) the stochastic error between the stochastic solution and the finite-volume solution and (2) the deterministic error due to the space–time discretization of the finite-volume scheme. That is, for a stochastic solution n^stoc, a finite-volume solution n^FV, and an exact true solution n^true, we can write

$$\begin{array}{}\text{(18)}& \underset{\text{total error}}{\underbrace{\Vert {\mathit{n}}^{\text{stoc}}-{\mathit{n}}^{\text{true}}\Vert }}\le \underset{\text{stochastic error}}{\underbrace{\Vert {\mathit{n}}^{\text{stoc}}-{\mathit{n}}^{\text{FV}}\Vert }}+\underset{\text{deterministic error}}{\underbrace{\Vert {\mathit{n}}^{\text{FV}}-{\mathit{n}}^{\text{true}}\Vert }}.\end{array}$$

In this section we focus on the stochastic error. We do not consider the refinement of Δx→0 or Δt→0 as it is well understood how the finite-volume methods converge to the true solution (deterministic error goes to zero) in these limits.

We present numerical examples of increasing complexity and discuss their convergence properties as the number of computational particles increases. Section 3.1 presents a one-dimensional test case with a soft-hat initial condition advected by a constant, uniform wind velocity to quantify the convergence as the number of computational particles increases. Section 3.2 simplifies the 1D test case to a uniform initial condition to study how the order of the advection scheme impacts convergence. Section 3.3 presents an idealized two-dimensional test case for solid-body rotational flow developed within WRF-PartMC using monotonic advection schemes. Finally, Sect. 3.4 shows stochastic particle transport for a realistic model domain and with realistic and evolving meteorological fields as simulated by the WRF-PartMC model.

To quantify the accuracy of the stochastic particle-resolved transport algorithm described above, we use the relative root mean square error (RRMSE) between two solutions n and n^′ as given by

$$\begin{array}{}\text{(19)}& \text{RRMSE}(\mathit{n},{\mathit{n}}^{\prime })={\displaystyle \frac{\sqrt{{\sum }_i^{N_{\mathrm{x}}}(n_i-n_i^{\prime }{)}^{\mathrm{2}}}}{\sqrt{{\sum }_i^{N_{\mathrm{x}}}(n_i^{\prime }{)}^{\mathrm{2}}}}}.\end{array}$$

To determine the mean and confidence intervals for the RRMSE, we ran an ensemble of simulations with different random seeds for the stochastic sampling algorithm. The RRMSE was computed for each simulation run, and then the overall mean and standard deviation were calculated, with the standard deviation being used to determine the 95 % confidence interval for the RRMSE mean.

3.1 One-dimensional test case: soft hat advected by uniform wind

We begin with the 1D soft-hat test case from Wicker and Skamarock (2002), with initial condition

$$\begin{array}{}\text{(20)}& n(x,t=\mathrm{0})=\mathrm{0.5}+{\displaystyle \frac{\mathrm{0.5}}{\mathrm{1}+{\exp }^{\mathrm{80}\left(\right|x-\mathrm{0.5}|-\mathrm{0.15})}}}\end{array}$$

on the periodic domain $x\in [\mathrm{0},\mathrm{1}]$. Equation (20) was modified from the expression in Wicker and Skamarock (2002) to include a background concentration so that there are some particles everywhere throughout the domain. The uniform velocity field was u=1 (all quantities are taken as dimensionless here). For all presented results, the number of grid cells was N_x=50 and the time step was Δt=0.008, resulting in a Courant number of 0.4. The simulation duration was T=2, giving two full revolutions of the domain. Simulation results were produced for first- to sixth-order advection schemes with no limiters applied.

Figure 1 shows the solution to the soft-hat problem after two revolutions (t=2) for the finite-volume method and one ensemble member of the stochastic method. Considering the finite-volume solutions, we observe that the even-order methods (second, fourth, sixth) produce more oscillatory solutions than the odd-order methods. This disparity between the even- and odd-order methods also occurs for the stochastic method, where the even-order solutions contain significantly more high-frequency noise than the odd-order solutions. We will analyze this phenomenon in more detail in Sect. 3.2.

We now turn to understanding the convergence of the stochastic particle method as the number of particles is increased. Figure 2a shows the ensemble mean error for the particle solution compared to the finite-volume solutions for each order of advection; i.e., this is the error due to using a finite particle number but does not include any spatial discretization error as that is present for the stochastic and finite-volume methods. As the number of computational particles per grid cell, N_p, increased, the solution converged to the deterministic finite-volume solution. The rate of convergence for these stochastic methods is expected to be $\mathrm{1}/\sqrt{N_{\mathrm{p}}}$ due to the central limit theorem and is denoted by the dashed line with slope $-\mathrm{1}/\mathrm{2}$. The stochastic error is largest for the even-order methods and smallest for the first-order method. In Sect. 3.2 we will show that this is because the odd-order methods benefit from the damping of high-frequency noise, and the first-order method has the lowest stochastic error because it has the most damping.

Figure 2 One-dimensional soft-hat test case (Sect. 3.1). (a) Relative root mean square error (RRMSE) between the stochastic particle solution and the deterministic finite-volume solution for first- to sixth-order advection with varying number of computational particles. The dashed line shows the expected $\mathrm{1}/\sqrt{N_{\mathrm{p}}}$ convergence rate. (b) Relative root mean square error (RRMSE) between the particle solution and the analytical solution for first- to sixth-order advection with varying number of computational particles per grid cell at t=2. Error bars denote the 95 % confidence interval as determined from an ensemble of 25 simulations.

Figure 3 One-dimensional uniform concentration advected by uniform wind: a single ensemble member of the stochastic solution is shown in blue for first- to sixth-order methods, the deterministic finite-volume solution is represented by the solid red line, and the analytical solution is shown as a dashed black line. The stochastic solution was simulated using N_p=10⁴ computational particles per grid cell.

Figure 2b shows the ensemble mean error for the particle solution compared to the analytical solution, i.e., the total error (finite-volume error plus stochastic error). As the number of particles increases, the ensemble mean error approaches a constant value, which is the error due to the finite-volume discretization. No matter how many computational particles are used, the total error cannot become smaller than the error introduced by the finite-volume discretization. The finite-volume error decreases in magnitude as the order of the advection method increases. For small values of N_p, the stochastic error dominates the total error. The odd-order methods (third- and fifth-order) have lower stochastic error than the even-order methods, which results in the total error converging to the finite-volume error with fewer computational particles. In contrast, the first-order method has such a large finite-volume error, shown with the poor comparison of the finite-volume solution to the true solution in Fig. 1, and such a low stochastic error, shown in Fig. 2, that the finite-volume error dominates the total error immediately. In all cases the finite-volume error could also be reduced by decreasing the grid cell size, following the standard convergence analysis of finite-volume methods (Durran, 2010).

In summary, these results show that the stochastic error of particle-resolved advection converges as expected with the rate of $\mathrm{1}/\sqrt{N_{\mathrm{p}}}$. Conservative even-order schemes exhibit high-frequency oscillations in the finite-volume solution that are compounded by high-frequency noise from the stochastic sampling. For the dissipative odd-order schemes, numerical dissipation damps the high-frequency oscillations, as will be shown in Sect. 3.2. We therefore recommend the use of dissipative (odd-order) advection schemes. We also note that the stochastic advection scheme will be especially useful in open domains where we have an outflow boundary condition. In this case, the artificial noise injected by the stochastic sampling will be advected out of the domain and will not accumulate.

3.2 One-dimensional test case: uniform concentration advected by uniform wind

The difference observed in Fig. 1 between the even- and odd-order solutions is of course due to the amount of numerical dissipation in the methods, where the even-order methods are conservative, while the odd-order methods have numerical dissipation of energy at high spatial frequencies (see, e.g., Durran, 2010, Sect. 3.3.2–3.3.3). To understand how this interacts with the stochastic particle solution, we derived a simple explicit model for the power spectrum of the stochastic solution. The details of this derivation are given in Appendices C–E. Briefly, we considered the uniform initial condition

$$\begin{array}{}\text{(21)}& n(x,t=\mathrm{0})=\mathrm{1}\end{array}$$

on the periodic domain $x\in [\mathrm{0},\mathrm{1}]$ with the uniform velocity field u=1. We used a computational volume of V=10 000, which thus means the stochastic particle system started with N_p=10 000 particles per grid cell.

Figure 3 shows the solution to the case with uniform concentration and uniform wind after two revolutions for the finite-volume method and for one ensemble member of the stochastic method using the same parameters as for Fig. 1, i.e., N_x=50 grid cells and a time step of Δt=0.008 resulting in a Courant number of 0.4. The exact solution to this problem is clearly $n(x,t)=\mathrm{1}$ for all x and t, and the finite-volume solution yields this solution exactly. However, as the stochastic method moves particles from one grid cell to the next, it samples a per-grid random number that is uncorrelated between grid cells and thus injects energy at high spatial frequencies. As we will show with further analysis, this energy may then be dissipated by the numerical dissipation of the spatial discretization, and the question is whether the dissipation can effectively dampen the energy injection. To answer this question we will derive a simple model for the power spectrum of the stochastic solution.

We start by writing ${\widehat{n}}_k$ for the discrete Fourier transform (DFT) of n_i and recalling the classical fact that the power spectrum of the spatially discretized system evolves according to

$$\begin{array}{}\text{(22)}& P_k^{\mathrm{\ell }+\mathrm{1}}=\exp \left(A_k\right)P_k^{\mathrm{\ell }},\end{array}$$

where $P_k^{\mathrm{\ell }}=|{\widehat{n}}_k^{\mathrm{\ell }}{|}^{\mathrm{2}}$ is the power at wavenumber k and time step ℓ, and A_k is the amplification factor at wavenumber k (see Appendix A for details). Figure 4 shows the amplification factors for the different spatial discretizations. From this, we see that the even-order methods have an amplification factor of zero at all wavenumbers, meaning that these methods are exactly conservative. In contrast, the odd-order methods have negative amplification factors at higher wavenumbers, showing that these methods will dissipate high-spatial-frequency components.

Figure 4 One-dimensional uniform test case (Sect. 3.2). (a) Amplification factors, A_k, for the first- to sixth-order spatial discretizations. See Eqs. (B8)–(B13) for details. (b) Excitation term, E_k. See Eq. (D12) for details.

To understand the interaction between the stochastic sampling and the spatial discretization dissipation, Appendix D derives a recurrence relation for the power spectrum of an approximation to the stochastic solution:

$$\begin{array}{}\text{(23)}& {\stackrel{\mathrm{̃}}{P}}_k^{\mathrm{\ell }+\mathrm{1}}=\exp \left(A_k\right){\stackrel{\mathrm{̃}}{P}}_k^{\mathrm{\ell }}+E_k,\end{array}$$

where ${\stackrel{\mathrm{̃}}{P}}_k^{\mathrm{\ell }}$ is the power at wavenumber k and time step ℓ of the approximate stochastic solution $\stackrel{\mathrm{̃}}{N}$, A_k is the amplification factor of the spatial discretization, and E_k is a stochastic excitation term (see Appendix D for details). Figure 4b plots the excitation term, and we see that it injects energy at higher wavenumbers, due to the uncorrelated random noise in each grid cell from the stochastic transport. By comparing panels (a) and (b) in Fig. 4, we see that the negative amplification factors at high wavenumbers will tend to suppress the energy injection.

To improve the clarity of results, we discretized with N_x=50 grid cells, a time step of Δt=0.00125, and total time T=2 to give 1600 time steps for two revolutions, and all simulations were run without limiters. Figure 5 shows the power spectrum of the stochastic solution for third- and fourth-order advection after 1, 100, 400, and 1600 time steps, with the black lines showing the model prediction for the power spectrum (see Appendix E for details). Here we see the constant energy injection at high wavenumbers, with the fourth-order method steadily scaling up the power spectrum by the excitation term at each time step. In contrast, the third-order method has dissipation at higher wavenumbers which partially suppresses the injected energy and eventually reaches an equilibrium. This serves to suppress the high-frequency noise in the solution and explains the difference between the even- and odd-order stochastic solutions in Fig. 1.

Figure 5 One-dimensional uniform test case (Sect. 3.2): mean power spectra for third- and fourth-order advection methods after 1, 100, 400, and 1600 time steps. The overlaid black lines indicate the model prediction for each method at each time step. Each stochastic case was repeated 100 times to obtain the mean power spectra.

This point is further emphasized in Fig. 6 where all the stochastic methods are compared after 1600 time steps. The conservative even-order schemes all fall on the same curve, increasing in power at higher frequencies. For the dissipative odd-order methods, high frequencies are damped. As expected from Fig. 4, the damping was least pronounced for the fifth-order method and most pronounced for first order. In general, stochastic methods are less stable than their finite-volume counterparts as the stochastic noise injects energy on average. Conservative even-order methods are unconditionally unstable due to this noise injection, because the scheme itself will never damp any of this additional energy.

Figure 6 One-dimensional uniform test case (Sect. 3.2): power spectra for first- to sixth-order advection methods after 1600 time steps (two full revolutions of the system). The overlaid black lines indicate the model predictions for each method.

Finally, to study the effect of spatial quantization where some particles move faster and some slower, causing variance in particle velocity and position (Sect. 2.9), let us consider the following example. If we assume a constant solution at all times (as in Appendix C), then the probability that a particle moves k grid cells is $\text{Binom}(k;N_{\mathrm{t}},p)$, where N_t is the number of time steps and p is the probability of moving each step, which is equal to the Courant number. To investigate this, we refined the grid spacing and time step by a factor of 10 to be Δx=0.002 (N_x=500) and Δt=0.0008, which preserves the Courant number of $C=p=\mathrm{0.4}$ of the original simulation, and we took T=1 (N_t=1250) for one revolution. Then, using the binomial distribution, the mean number of grid cells moved in one revolution is $\mathit{\mu }=N_{\mathrm{t}}p=\mathrm{500}$, which is an exact approximation (zero bias), while the standard deviation is $\mathit{\sigma }=\sqrt{N_{\mathrm{t}}p(\mathrm{1}-p)}=\mathrm{17.3}$. This corresponds to a physical distance of $x_{\mathit{\sigma }}=\mathit{\sigma }\mathrm{\Delta }x=\mathrm{0.035}$. To understand the limiting behavior, we can use $N_{\mathrm{t}}=T/\mathrm{\Delta }t$ and $p=C=u\mathrm{\Delta }t/\mathrm{\Delta }x$ to rewrite x_σ as

$$\begin{array}{}\text{(24)}& {\displaystyle }{x}_{\mathit{\sigma }}=\mathrm{\Delta }x\sqrt{{\displaystyle \frac{T}{\mathrm{\Delta }t}}{\displaystyle \frac{u\mathrm{\Delta }t}{\mathrm{\Delta }x}}\left(\mathrm{1}-C\right)}=\sqrt{Tu(\mathrm{1}-C)\mathrm{\Delta }x}.\end{array}$$

Now consider refining the grid (Δx→0) and time step (Δt→0) while keeping constant the Courant number C, the final time T, and the velocity u. In this limit, we can see that x_σ→0, so that the numerical diffusion of particles caused by the stochastic method vanishes.

Figure 7 shows the numerical result of the diffusion after one revolution for the particles originating in grid cell 250 (in blue), with the analytical binomial model shown in red. During sampling, some particles will travel faster and some will travel slower, resulting in the binomial distribution of particles around the mean position.

Figure 7 One-dimensional uniform test case (Sect. 3.2) for the effect of spatial quantization on the stochastic solution. (a) Initial condition showing the uniform number concentration in all grid cells and the location of number concentration originating in grid cell 250 at about x=0.5. (b) Number concentration of particles originating from grid cell 250 after one revolution at time t=1 (blue points) with the analytical binomial model (solid red line). The vertical dashed line separates particles that moved too fast (to the right) and too slow (to the left).

3.3 Two-dimensional test case: Gaussian cone advected by solid-body rotational wind field

To test the schemes in 2D, we used a scalar advection problem modified from Wicker and Skamarock (2002), where a Gaussian cone is advected in a square domain by a prescribed solid-body rotation flow. Simulations were conducted using third- and fifth-order monotonic advection schemes. Figure 8 shows the initial conditions. The domain is 100 × 100 nondimensional units, and the velocity field is defined as $u(x,y)=-\mathit{\omega }\left(y-\mathrm{50}\right)$ and $v(x,y)=\mathit{\omega }(x-\mathrm{50})$, where $\mathit{\omega }=\frac{\mathrm{2}\mathit{\pi }}{\mathrm{628}}$. We took $\mathrm{\Delta }x=\mathrm{\Delta }y=\mathrm{1}$ and Δt=0.5, so that one full rotation requires 1256 time steps. The maximum Courant number was 0.5. The initial particle mixing ratio was given as

$$\begin{array}{}\text{(25)}& q(x,y)=max\left({\mathrm{10}}^{\mathrm{10}}\exp \left(-{\left({\displaystyle \frac{r}{r_{\mathrm{0}}}}\right)}^{\mathrm{2}}\right),{\mathrm{10}}^{-\mathrm{15}}\right),\end{array}$$

where $r=\sqrt{(x-\mathrm{50}{)}^{\mathrm{2}}+(y-\mathrm{75}{)}^{\mathrm{2}}}$ and r₀=6. The grid cell average values were constructed using 5×5 point Gaussian quadrature.

Figure 8 Two-dimensional test case (Sect. 3.3): (a) rotational wind field, (b) true solution of the Gaussian cone after one complete revolution, and (c) true solution of the red outlined region in panel (b).

Figure 9 shows the solution after one revolution for the region of interest with 100, 1000, and 10 000 computational particles per grid cell as well as the finite-volume solution. As the number of computational particles increased, the solution became less noisy and more similar to the finite-volume solution. This is quantified in Fig. 10, which shows the error of the particle solution compared to the finite-volume solution. As expected, the stochastic error of the Monte Carlo method has a rate of convergence of $\mathrm{1}/\sqrt{N_{\mathrm{p}}}$. As the convergence of finite-volume solutions to the analytical solution is well studied (Wicker and Skamarock, 2002), we do not include results showing convergence in Δx and Δy.

Figure 9 Two-dimensional test case (Sect. 3.3): stochastic particle solution for 100, 1000, and 10 000 computational particles per grid cell and the finite-volume solution after one revolution for the region shown in red in Fig. 8b.

Figure 10 Two-dimensional test case (Sect. 3.3): relative root mean square error (RRMSE) between the particle solution and the finite-volume solution for third- and fifth-order monotonic advection. Error bars indicate the 95 % confidence interval from 10 simulations. The black reference line indicates the theoretical convergence rate with slope $\mathrm{1}/\sqrt{N_{\mathrm{p}}}$.

3.4 Three-dimensional test case: plume transported by WRF-simulated meteorology

For this simulation we used WRF to fully simulate the meteorology, resulting in an evolving velocity field in 3D. We prescribed an idealized initial condition of particle mixing ratio and gas tracer mixing ratio for the model domain of Northern California. The gas tracer mixing ratio was used as a proxy for the solution of the finite-volume method. The domain comprised $\mathrm{170}\times \mathrm{160}\times \mathrm{40}$ grid cells, with $\mathrm{\Delta }x=\mathrm{\Delta }y=\mathrm{4}$ km and Δz increasing logarithmically from an average value of 55 m near the surface to 650 m near the top of the model domain. The model time step was set to Δt=20 s, ensuring that the sum of particle cell transfer probabilities did not exceed 1. For this case, an initial cloud of aerosol particle mixing ratio and gas mixing ratio was determined by

$$\begin{array}{ll}\text{(26)}& {\displaystyle }& {\displaystyle }{q}_{\text{grid}}(x,y,z)={\displaystyle }& {\displaystyle }max\left({\mathrm{10}}^{\mathrm{10}}\exp \left(-\left(\sqrt{\begin{array}{c}{\left(\frac{x-x_{\mathrm{0}}}{r_x}\right)}^{\mathrm{2}}+{\left(\frac{y-y_{\mathrm{0}}}{r_y}\right)}^{\mathrm{2}}\\ +{\left({\displaystyle \frac{z-z_{\mathrm{0}}}{r_z}}\right)}^{\mathrm{2}}\end{array}}\right)\right),{\mathrm{10}}^{-\mathrm{15}}\right),\end{array}$$

where $r_x=r_y=\mathrm{6}$ and r_z=4, and the cloud is centered at grid cell x₀=75, y₀=75, and z₀=1. Here $q_{\text{grid}}(x,y,z)$ is specified in grid coordinates (each grid cell is a square of size $\mathrm{1}\times \mathrm{1}\times \mathrm{1}$ grid units) before being transformed to physical coordinates for the simulation. Figure 11a shows the initial condition described by Eq. (26) at the lowest model layer. The initial condition was advected by the dynamic meteorology over a 12 h period beginning at 00:00 UTC on 7 June 2010 using a time step of Δt=20 s. Meteorological initial and boundary conditions were based on analyses from the National Centers for Environmental Prediction's North American Mesoscale Forecast System (NAM) model. The temporal evolution of the wind field is shown in Fig. 11b–d in increments of 6 h. Gases and particles are subject only to advection and do not experience turbulent diffusion or any removal processes. Gas and aerosol boundary conditions were prescribed from initial values given in Eq. (26). When flow enters the domain at a boundary grid cell, the prescribed value is applied. Conversely, when flow exits the domain, the boundary grid cell assumes a zero gradient condition, consistent with the host model WRF. Simulations were conducted using third- and fifth-order monotonic advection.

Figure 11 Three-dimensional test case (Sect. 3.4): (a) the initial condition and (b–d) snapshots of the wind velocity field at times t=0, 6, and 12 h in the lowest model layer.

Figure 12 Three-dimensional test case (Sect. 3.4): lowest layer mixing ratios after 12 h of simulation for 10, 100, and 1000 computational particles per grid cell and the deterministic finite-volume solution reference solution.

Figure 13 Three-dimensional test case (Sect. 3.4): convergence of the relative root mean square error (RRMSE) between the stochastic solution and the finite-volume solution as the number of computational particles per grid cell increases. Error bars show the 95 % confidence interval from an ensemble of five simulations.

Figure 12 shows the solution after 12 h for a varying number of computational particles per grid cell, with the finite-volume solution for comparison. The simulation with 10 particles per grid cell is noisy as expected, capturing only general features of the particle number mixing ratio. As the number of computational particles was increased, the particle number mixing-ratio field became smoother and similar to the finite-volume solution.

Figure 13 shows the convergence of the three-dimensional test case for third-order monotonic advection. As the number of computational particles increased, the error when compared to the finite-volume solution converged at the expected rate of $\mathrm{1}/\sqrt{N_{\mathrm{p}}}$. Due to the stochastic nature of the problem, monotonic limiters may be applied to the particle number mixing-ratio field that does not exist in the finite-volume solution. As a result, a perfect $\mathrm{1}/\sqrt{N_{\mathrm{p}}}$ convergence rate is not expected.

Figure 14 confirms that the stochastic solution converges to the finite-volume solution for the three-dimensional test case and that the variance decreases as the number of computational particles increases. For reference, Fig. 14a shows an x–y cross section of the mean mixing ratio in the lowest model layer at t=12 h. The mean mixing ratio was calculated by averaging the stochastic solution over five simulations using N_p=100 computational particles.

Figure 14 Three-dimensional test case (Sect. 3.4): (a) ensemble mean mixing ratio averaged over five simulations after 12 h for the lowest model layer with N_p=100 computational particles per grid cell; (b) vertical profile of the mixing ratio on a logarithmic scale for stochastic solutions of N_p=10, 100, and 1000 computational particles per grid cell at x=75 and y=75 at time t=12 h; (c) x transect at y=25 for stochastic solutions of N_p=10, 100, and 1000 computational particles per grid cell at t=12 h; and (d) time series at x=75 and y=25 for N_p=10, 100, and 1000 computational particles per grid cell. The finite-volume solutions for the profile, transect, and time series are denoted by black lines. Points show means of five simulations, and error bars denote the corresponding 95 % confidence intervals.

Figure 14b–d show different transects through the three-dimensional space and time. The star in Fig. 14a marks the location of the vertical mixing-ratio profile (log-scaled) in Fig. 14b and the time series shown in Fig. 14d. The red line denotes the transect shown in Fig. 14c. The finite-volume solution is compared to the ensemble mean of 10, 100, and 1000 computational particles with error bars denoting the 95 % confidence interval. As the number of particles increased, the variance decreased and the solution converged to the finite-volume solution.

In Sect. 2.7, we discussed sampling complexities due to different computational volumes, grid cell volumes, and air densities. When these quantities substantially differ in adjacent grid cells, it could lead to undersampling of rare particle types. In our three-dimensional example, the largest ratio in density was 1.29, and the largest grid cell volume ratio was 1.96. For most of the grid cells, these ratios were closer to 1, indicated by domain average ratios of 1.01 and 1.11, respectively, at t=12 h.

To investigate whether undersampling occurred in practice, we ran the same scenario but sampled the particle diameter (a 1D attribute carried by particles) from a log-normal size distribution so that rare large and small particles existed while most computational particles resided in the center of the size distribution. We then compared the final size distributions with the initial size distributions to determine to what extent the rare large and small particles were systematically lost due to undersampling.

Figure 15a shows the locations for the initial and final size distribution plots. The locations of the initial and final points were chosen so that the final point is downwind of the initial point. All grid cells were initialized with 100 computational particles drawn from a single log-normal mode, all with a constant geometric mean diameter and geometric standard deviation where only the magnitude of the distribution was adjusted. Figure 15b shows the normalized mean particle size distribution at the initial time and the final time at two single grid cells. Each distribution was averaged over five ensemble runs.

Figure 15 (a) Locations and patches where the normalized size distributions were computed. (b) Normalized size distributions at single grid cell locations shown in panel (a) using an ensemble of five simulations. (c) Mean normalized size distributions for the 15×15 subdomain patches shown in panel (a). Solid blue lines show the initial distributions computed at t=0 h and at the initial point/region, while dashed red lines show the final distributions computed at t=12 h and at the final point/region.

As we see from Fig. 15b, the size distribution at the final time was similar to that at the initial time, with some stochastic noise. To reduce the stochastic noise, Fig. 15c shows the normalized mean particle size distribution at the initial time and the final time for two 15×15 grid cell patches surrounding the points chosen for Fig. 15b. Here the normalized size distributions were nearly identical, indicating that the size distribution information was not lost in the sampling procedure.

4 Conclusions

In this paper we presented the development of a stochastic particle advection method and demonstrated its performance for particle-resolved atmospheric aerosol transport in the combined WRF-PartMC model. The method is based on finite-volume advection schemes but interprets the fluxes as probabilities of particle transport, which can then be stochastically sampled. We analyzed the method in the one-dimensional setting to show that the stochastic particle sampling injects noise at high spatial frequencies, and so the method performs best when using dissipative finite-volume discretizations, such as the third- and fifth-order schemes used in WRF.

We applied the new method in WRF-PartMC with the existing monotonic limiter for the fifth-order scheme and a new limiter for third order. We considered two test cases: a solid-body rotational wind field in 2D and an atmospherically relevant dynamic wind field over complex terrain in 3D. In both cases we observed the expected rates of convergence of the stochastic particle transport to the finite-volume solution as the number of computational particles per grid cell was increased. For these examples, significant stochastic noise was evident in simulations with 100 computational particles per grid cell, but stochastic noise was found to be less than 10 % for simulations with 1000 particles per grid cell. This is considered a reasonable number of computational particles for large-scale WRF-PartMC simulations, as these simulations typically use on the order of 10 000 computational particles to accurately capture properties of the aerosol mixing state (Gasparik et al., 2020).

The value of this work is to enable direct comparison of particle-resolved aerosol representations to models that use approximate aerosol representations with simplified assumptions regarding size and composition (e.g., internally mixed modes or bins). Because the stochastic particle method is based on the same finite-volume schemes used for the approximate representations, model comparisons can isolate the differences arising due to aerosol representation. Additionally, the new stochastic transport scheme allows the WRF-PartMC model to be used on the regional scale to quantify the impact of aerosol mixing state on climate-relevant aerosol properties, such as aerosol absorption and CCN concentration, and to compare these findings to existing studies (Matsui et al., 2013; Zhang et al., 2014; Zhu et al., 2016).

Appendix A: One-dimensional advection in the frequency domain

To understand the behavior of the 1D deterministic and stochastic numerical methods it is helpful to write them in the frequency domain. To do this, we start in this section by considering only the deterministic (finite-volume) case. We will then extend this to the stochastic case in the next section. We will use the vector notation

$$\begin{array}{}\text{(A1)}& {\displaystyle }\mathit{n}& {\displaystyle }=[n_{\mathrm{0}},n_{\mathrm{1}},\mathrm{\dots },n_{N_{\mathrm{x}}-\mathrm{1}}],\text{(A2)}& {\displaystyle }\mathit{f}& {\displaystyle }=[f_{\frac{\mathrm{1}}{\mathrm{2}}},f_{\mathrm{3}/\mathrm{2}},\mathrm{\dots },f_{N_{\mathrm{x}}-\frac{\mathrm{1}}{\mathrm{2}}}].\end{array}$$

We assume periodicity, so $n_i=n_{i+N_{\mathrm{x}}}$ and $f_{i-\frac{\mathrm{1}}{\mathrm{2}}}=f_{i-\frac{\mathrm{1}}{\mathrm{2}}+N_{\mathrm{x}}}$ for any i. Similarly, we encode the finite-difference stencils as vectors:

$$\begin{array}{}\text{(A3)}& {\displaystyle }{\mathit{r}}^{\text{1st}}& {\displaystyle }=[\mathrm{1},\mathrm{0},\mathrm{\dots },\mathrm{0}],\text{(A4)}& {\displaystyle }{\mathit{r}}^{\text{2nd}}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}[\mathrm{1},\mathrm{0},\mathrm{\dots },\mathrm{0},\mathrm{1}],\text{(A5)}& {\displaystyle }{\mathit{r}}^{\text{3rd}}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathrm{6}}}[\mathrm{5},-\mathrm{1},\mathrm{0},\mathrm{\dots },\mathrm{0},\mathrm{2}],\text{(A6)}& {\displaystyle }{\mathit{r}}^{\text{4th}}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathrm{12}}}[\mathrm{7},-\mathrm{1},\mathrm{0},\mathrm{\dots },\mathrm{0},-\mathrm{1},\mathrm{7}],\text{(A7)}& {\displaystyle }{\mathit{r}}^{\text{5th}}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathrm{60}}}[\mathrm{47},-\mathrm{13},\mathrm{2},\mathrm{0},\mathrm{\dots },\mathrm{0},-\mathrm{3},\mathrm{27}],\text{(A8)}& {\displaystyle }{\mathit{r}}^{\text{6th}}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathrm{60}}}[\mathrm{37},-\mathrm{8},\mathrm{1},\mathrm{0},\mathrm{\dots },\mathrm{0},\mathrm{1},-\mathrm{8},\mathrm{37}].\end{array}$$

This allows us to express the flux Eqs. (3)–(8) via a convolution:

$$\begin{array}{}\text{(A9)}& {\displaystyle }\mathit{f}& {\displaystyle }=u\mathit{r}*\mathit{n},\text{(A10)}& {\displaystyle }{f}_{i+\frac{\mathrm{1}}{\mathrm{2}}}& {\displaystyle }=u\sum _{j=\mathrm{0}}^{N_{\mathrm{x}}-\mathrm{1}}{r}_{i-j}{n}_j.\end{array}$$

Next, define the finite-difference stencil

$$\begin{array}{}\text{(A11)}& {\displaystyle }\mathit{d}{\displaystyle }=[\mathrm{1},-\mathrm{1},\mathrm{0},\mathrm{\dots },\mathrm{0}]\end{array}$$

so we can approximate the spatial derivative as

$$\begin{array}{}\text{(A12)}& {\displaystyle \frac{\partial \mathit{n}}{\partial x}}\approx {\displaystyle \frac{\mathrm{1}}{\mathrm{\Delta }x}}\mathit{d}*\mathit{n}.\end{array}$$

Using this we can write the spatially discretized advection Eq. (2) as

$$\begin{array}{}\text{(A13)}& {\displaystyle \frac{\partial \mathit{n}}{\partial t}}& {\displaystyle }=-{\displaystyle \frac{\mathrm{1}}{\mathrm{\Delta }x}}\mathit{d}*\mathit{f}\text{(A14)}& {\displaystyle }& {\displaystyle }=-{\displaystyle \frac{u}{\mathrm{\Delta }x}}\mathit{d}*\mathit{r}*\mathit{n}.\end{array}$$

We denote the discrete Fourier transform (DFT) using a hat, so $\widehat{n}=\mathcal{F}\left(n\right)$ and similarly for other variables, and recall that the DFT is given by

$$\begin{array}{}\text{(A15)}& {\widehat{n}}_k=\sum _{j=\mathrm{0}}^{N_{\mathrm{x}}-\mathrm{1}}{n}_j\exp (-i\mathrm{2}\mathit{\pi }jk/N_{\mathrm{x}}),\end{array}$$

where i is the imaginary unit. Taking the DFT of Eq. (A14) gives

$$\begin{array}{}\text{(A16)}& {\displaystyle \frac{\partial {\widehat{n}}_k}{\partial t}}=-{\displaystyle \frac{u}{\mathrm{\Delta }x}}{d}_k{r}_k{n}_k\end{array}$$

for each wavenumber k. The solution over one time step is then given by

$$\begin{array}{}\text{(A17)}& {\widehat{n}}_k^{\mathrm{\ell }+\mathrm{1}}=\exp (-C{\widehat{d}}_k{\widehat{r}}_k){\widehat{n}}_k^{\mathrm{\ell }},\end{array}$$

where C is the Courant number given by

$$\begin{array}{}\text{(A18)}& {\displaystyle }C{\displaystyle }={\displaystyle \frac{u\mathrm{\Delta }t}{\mathrm{\Delta }x}}.\end{array}$$

Composing ℓ time steps gives the solution at time step ℓ as

$$\begin{array}{}\text{(A19)}& {\widehat{n}}_k^{\mathrm{\ell }}=\exp (-\mathrm{\ell }C{\widehat{d}}_k{\widehat{r}}_k){\widehat{n}}_k^{\mathrm{0}}.\end{array}$$

To understand the numerical effect of the finite-difference approximation we can compute the evolution of the power spectrum of the solution. The power spectrum is given by

$$\begin{array}{}\text{(A20)}& P_k=|{\widehat{n}}_k{|}^{\mathrm{2}},\end{array}$$

and the evolution of the power spectrum over one time step is given by

$$\begin{array}{}\text{(A21)}& {\displaystyle }|{\widehat{n}}_k^{\mathrm{\ell }+\mathrm{1}}{|}^{\mathrm{2}}& {\displaystyle }={\widehat{n}}^{\mathrm{\ell }+\mathrm{1}}{\widehat{n}}_k^{(\mathrm{\ell }+\mathrm{1})*}\text{(A22)}& {\displaystyle }& {\displaystyle }=\left(\exp \right(-C{\widehat{d}}_k{\widehat{r}}_k\left){\widehat{n}}_k^{\mathrm{\ell }}\right)\left(\exp \right(-C{\widehat{d}}_k{\widehat{r}}_k){\widehat{n}}_k^{\mathrm{\ell }}{)}^{*}\text{(A23)}& {\displaystyle }& {\displaystyle }=\exp (-C{\widehat{d}}_k{\widehat{r}}_k)\exp (-C{\widehat{d}}_k^{*}{\widehat{r}}_k^{*}){\widehat{n}}_k^{\mathrm{\ell }}{\widehat{n}}_k^{\mathrm{\ell }*}\text{(A24)}& {\displaystyle }& {\displaystyle }=\exp (-\mathrm{2}C\text{Re}({\widehat{d}}_k{\widehat{r}}_k\left)\right)|{\widehat{n}}_k^{\mathrm{\ell }}{|}^{\mathrm{2}}.\end{array}$$

The energy amplification of the method is thus given by

$$\begin{array}{}\text{(A25)}& A_k=-\mathrm{2}C\text{Re}\left({\widehat{d}}_k{\widehat{r}}_k\right),\end{array}$$

and we can write the power spectrum evolution as

$$\begin{array}{}\text{(A26)}& P_k^{\mathrm{\ell }+\mathrm{1}}=\exp \left(A_k\right)P_k^{\mathrm{\ell }}.\end{array}$$

If A_k is zero then the method conserves the energy in wavenumber k, while negative values indicate that the method will dissipate energy with each time step.

Appendix B: DFT of finite-difference stencils

The DFT of the finite-difference stencils d and r is found by applying Eq. (A15) to Eq. (A11) and Eqs. (A3)–(A8). This gives

$$\begin{array}{}\text{(B1)}& {\widehat{d}}_k=\mathrm{1}-\exp (-i\mathrm{2}\mathit{\pi }k/N_{\mathrm{x}})\end{array}$$

and

$$\begin{array}{ll}\text{(B2)}& {\displaystyle }{\widehat{r}}_k^{\text{1st}}& {\displaystyle }=\mathrm{1},\text{(B3)}& {\displaystyle }{\widehat{r}}_k^{\text{2nd}}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(\exp \right(i\mathrm{2}\mathit{\pi }k/N_{\mathrm{x}})+\mathrm{1}),\text{(B4)}& {\displaystyle }{\widehat{r}}_k^{\text{3rd}}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathrm{6}}}\left(\mathrm{2}\exp \right(i\mathrm{2}\mathit{\pi }k/N_{\mathrm{x}})+\mathrm{5}-\exp (-i\mathrm{2}\mathit{\pi }k/N_{\mathrm{x}}\left)\right),{\displaystyle }{\widehat{r}}_k^{\text{4th}}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathrm{12}}}(-\exp (i\mathrm{2}\mathit{\pi }\mathrm{2}k/N_{\mathrm{x}})+\mathrm{7}\exp (i\mathrm{2}\mathit{\pi }k/N_{\mathrm{x}})\\ \text{(B5)}& {\displaystyle }& {\displaystyle }+\mathrm{7}-\exp (-i\mathrm{2}\mathit{\pi }k/N_{\mathrm{x}})),{\displaystyle }{\widehat{r}}_k^{\text{5th}}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathrm{60}}}(-\mathrm{3}\exp (i\mathrm{2}\mathit{\pi }\mathrm{2}k/N_{\mathrm{x}})+\mathrm{27}\exp (i\mathrm{2}\mathit{\pi }k/N_{\mathrm{x}})\\ {\displaystyle }& {\displaystyle }+\mathrm{47}-\mathrm{13}\exp (-i\mathrm{2}\mathit{\pi }k/N_{\mathrm{x}})\\ \text{(B6)}& {\displaystyle }& {\displaystyle }+\mathrm{2}\exp (-i\mathrm{2}\mathit{\pi }\mathrm{2}k/N_{\mathrm{x}})),{\displaystyle }{\widehat{r}}_k^{\text{6th}}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathrm{60}}}\left(\exp \right(i\mathrm{2}\mathit{\pi }\mathrm{3}k/N_{\mathrm{x}})-\mathrm{8}\exp (i\mathrm{2}\mathit{\pi }\mathrm{2}k/N_{\mathrm{x}})\\ {\displaystyle }& {\displaystyle }+\mathrm{37}\exp (i\mathrm{2}\mathit{\pi }k/N_{\mathrm{x}})+\mathrm{37}\\ {\displaystyle }& {\displaystyle }-\mathrm{8}\exp (-i\mathrm{2}\mathit{\pi }k/N_{\mathrm{x}})\\ \text{(B7)}& {\displaystyle }& {\displaystyle }+\exp (-i\mathrm{2}\mathit{\pi }\mathrm{2}k/N_{\mathrm{x}})).\end{array}$$

The amplification A_k of the above stencils can now be found by evaluating Eq. (A25) to give

$$\begin{array}{ll}\text{(B8)}& {\displaystyle }{A}_k^{\text{1st}}& {\displaystyle }=C(-\mathrm{2}+\mathrm{2}\cos (\mathrm{2}\mathit{\pi }k/N_{\mathrm{x}}\left)\right),\text{(B9)}& {\displaystyle }{A}_k^{\text{2nd}}& {\displaystyle }=\mathrm{0},\text{(B10)}& {\displaystyle }{A}_k^{\text{3rd}}& {\displaystyle }={\displaystyle \frac{C}{\mathrm{3}}}(-\mathrm{3}+\mathrm{4}\cos (\mathrm{2}\mathit{\pi }k/N_{\mathrm{x}})-\cos (\mathrm{2}\mathit{\pi }\mathrm{2}k/N_{\mathrm{x}}\left)\right),\text{(B11)}& {\displaystyle }{A}_k^{\text{4th}}& {\displaystyle }=\mathrm{0},{\displaystyle }{A}_k^{\text{5th}}& {\displaystyle }={\displaystyle \frac{C}{\mathrm{30}}}(-\mathrm{20}+\mathrm{30}\cos (\mathrm{2}\mathit{\pi }k/N_{\mathrm{x}})\\ \text{(B12)}& {\displaystyle }& {\displaystyle }-\mathrm{12}\cos (\mathrm{2}\mathit{\pi }\mathrm{2}k/N_{\mathrm{x}})+\mathrm{2}\cos (\mathrm{2}\mathit{\pi }\mathrm{3}k/N_{\mathrm{x}})),\text{(B13)}& {\displaystyle }{A}_k^{\text{6th}}& {\displaystyle }=\mathrm{0}.\end{array}$$

Appendix C: An approximate model for particle advection in 1D

We want to model the stochastic particle advection process as a deterministic advection process with some additional noise. We start by writing Eq. (15) as

$$\begin{array}{}\text{(C1)}& {\displaystyle }{F}_{i+\frac{\mathrm{1}}{\mathrm{2}}}^{\mathrm{\ell }}& {\displaystyle }=\text{Binom}(N_i^{\mathrm{\ell }},p_{i+\frac{\mathrm{1}}{\mathrm{2}}}^{\mathrm{\ell }})\text{(C2)}& {\displaystyle }& {\displaystyle }=\mathrm{E}\left[F_{i+\frac{\mathrm{1}}{\mathrm{2}}}^{\mathrm{\ell }}\right]+S_{i+\frac{\mathrm{1}}{\mathrm{2}}}^{\mathrm{\ell }}\text{(C3)}& {\displaystyle }& {\displaystyle }=p_{i+\frac{\mathrm{1}}{\mathrm{2}}}^{\mathrm{\ell }}{N}_i^{\mathrm{\ell }}+S_{i+\frac{\mathrm{1}}{\mathrm{2}}}^{\mathrm{\ell }}\text{(C4)}& {\displaystyle }& {\displaystyle }={\stackrel{\mathrm{‾}}{F}}_{i+\frac{\mathrm{1}}{\mathrm{2}}}^{\mathrm{\ell }}+S_{i+\frac{\mathrm{1}}{\mathrm{2}}}^{\mathrm{\ell }},\end{array}$$

where ${\stackrel{\mathrm{‾}}{F}}_{i+\frac{\mathrm{1}}{\mathrm{2}}}^{\mathrm{\ell }}$ is the deterministic mean flux and $S_{i+\frac{\mathrm{1}}{\mathrm{2}}}^{\mathrm{\ell }}$ is a zero-mean random variable representing the stochastic noise, given by

$$\begin{array}{}\text{(C5)}& S_{i+\frac{\mathrm{1}}{\mathrm{2}}}^{\mathrm{\ell }}=\text{Binom}(N_i^{\mathrm{\ell }},p_{i+\frac{\mathrm{1}}{\mathrm{2}}}^{\mathrm{\ell }})-{\stackrel{\mathrm{‾}}{F}}_{i+\frac{\mathrm{1}}{\mathrm{2}}}^{\mathrm{\ell }}.\end{array}$$

We approximate this stochastic noise by assuming that it is sampled from a constant uniform particle state with exactly $\stackrel{\mathrm{˘}}{N}$ particles per grid cell. From Eq. (12) we have

$$\begin{array}{}\text{(C6)}& \stackrel{\mathrm{˘}}{n}={\displaystyle \frac{\stackrel{\mathrm{˘}}{N}}{V}},\end{array}$$

and because the velocity u is constant and uniform the discretized flux is given by

$$\begin{array}{}\text{(C7)}& {\stackrel{\mathrm{˘}}{f}}^{††}=u\stackrel{\mathrm{˘}}{n}.\end{array}$$

From Eqs. (13) and (14) we then have

$$\begin{array}{}\text{(C8)}& {\displaystyle }\stackrel{\mathrm{˘}}{\stackrel{\mathrm{‾}}{F}}& {\displaystyle }=V{\displaystyle \frac{\mathrm{\Delta }t}{\mathrm{\Delta }x}}{\stackrel{\mathrm{˘}}{f}}^{††}\text{(C9)}& {\displaystyle }& {\displaystyle }=V{\displaystyle \frac{\mathrm{\Delta }t}{\mathrm{\Delta }x}}u{\displaystyle \frac{\stackrel{\mathrm{˘}}{N}}{V}}\text{(C10)}& {\displaystyle }& {\displaystyle }=C\stackrel{\mathrm{˘}}{N},\text{(C11)}& {\displaystyle }\stackrel{\mathrm{˘}}{p}& {\displaystyle }={\displaystyle \frac{\stackrel{\mathrm{˘}}{\stackrel{\mathrm{‾}}{F}}}{\stackrel{\mathrm{˘}}{N}}}\text{(C12)}& {\displaystyle }& {\displaystyle }=C.\end{array}$$

We can thus write the approximate stochastic noise by modifying Eq. (C5) to give

$$\begin{array}{}\text{(C13)}& {\displaystyle }{\stackrel{\mathrm{˘}}{S}}_i& {\displaystyle }=\text{Binom}(\stackrel{\mathrm{˘}}{N},\stackrel{\mathrm{˘}}{p})-\stackrel{\mathrm{˘}}{\stackrel{\mathrm{‾}}{F}}\text{(C14)}& {\displaystyle }& {\displaystyle }=\text{Binom}(\stackrel{\mathrm{˘}}{N},C)-C\stackrel{\mathrm{˘}}{N}.\end{array}$$

We want to write the approximate stochastic model in the frequency domain by taking a DFT. It is thus helpful to rewrite the equations in vector form, as we did in Sect. A. Similarly to Eqs. (A1) and (A2), we can write the particle counts $N_i^{\mathrm{\ell }}$ and particle fluxes $F_{i-\frac{\mathrm{1}}{\mathrm{2}}}^{\mathrm{\ell }}$ as vectors N^ℓ and F^ℓ, and we can also do the same for other variables such as the average particle flux ${\stackrel{\mathrm{‾}}{F}}_{i+\frac{\mathrm{1}}{\mathrm{2}}}^{\mathrm{\ell }}$ and probabilities $p_{i+\frac{\mathrm{1}}{\mathrm{2}}}^{\mathrm{\ell }}$.

Using the above vector notation and the difference stencil Eq. (A11), we can write the temporal update Eq. (16) as

$$\begin{array}{}\text{(C15)}& {\displaystyle }{N}_i^{\mathrm{\ell }+\mathrm{1}}& {\displaystyle }=N_i^{\mathrm{\ell }}-F_{i+\frac{\mathrm{1}}{\mathrm{2}}}^{\mathrm{\ell }}+F_{i-\frac{\mathrm{1}}{\mathrm{2}}}^{\mathrm{\ell }},\text{(C16)}& {\displaystyle }{N}^{\mathrm{\ell }+\mathrm{1}}& {\displaystyle }=N^{\mathrm{\ell }}+d*F^{\mathrm{\ell }}\text{(C17)}& {\displaystyle }& {\displaystyle }=N^{\mathrm{\ell }}+d*{\stackrel{\mathrm{‾}}{F}}^{\mathrm{\ell }}+d*S^{\mathrm{\ell }}.\end{array}$$

Taking the DFT now gives

$$\begin{array}{}\text{(C18)}& {\displaystyle }{\widehat{N}}_k^{\mathrm{\ell }+\mathrm{1}}& {\displaystyle }={\widehat{N}}_k^{\mathrm{\ell }}+{\widehat{d}}_k{\widehat{F}}_k^{\mathrm{\ell }}+{\widehat{d}}_k{\widehat{S}}_k^{\mathrm{\ell }}\text{(C19)}& {\displaystyle }& {\displaystyle }\approx \exp (-C{\widehat{d}}_k{\widehat{r}}_k){\widehat{N}}_k^{\mathrm{\ell }}+{\widehat{d}}_k{\widehat{S}}_k^{\mathrm{\ell }}\text{(C20)}& {\displaystyle }& {\displaystyle }\approx \exp (-C{\widehat{d}}_k{\widehat{r}}_k){\widehat{N}}_k^{\mathrm{\ell }}+{\widehat{d}}_k{\widehat{\stackrel{\mathrm{˘}}{S}}}_k.\end{array}$$

In Eq. (C19) we approximated the update of the deterministic component with the exact solution of the deterministic advection equation, as in Eq. (A17). That is, we approximated the Runge–Kutta time step update with the exact solution. We then approximated the update of the stochastic component in Eq. (C20) by using the approximate stochastic noise $\stackrel{\mathrm{˘}}{S}$.

Defining $\stackrel{\mathrm{̃}}{N}$ to be the solution of the approximate model, we can write the final approximate model from Eqs. (C20) and (C14) as

$$\begin{array}{}\text{(C21)}& {\displaystyle }{\widehat{\stackrel{\mathrm{̃}}{N}}}_k^{\mathrm{\ell }+\mathrm{1}}& {\displaystyle }=\exp (-C{\widehat{d}}_k{\widehat{r}}_k){\widehat{\stackrel{\mathrm{̃}}{N}}}_k^{\mathrm{\ell }}+{\widehat{d}}_k{\widehat{\stackrel{\mathrm{˘}}{S}}}_k,\text{(C22)}& {\displaystyle }{\stackrel{\mathrm{˘}}{S}}_i& {\displaystyle }=\text{Binom}(\stackrel{\mathrm{˘}}{N},C)-C\stackrel{\mathrm{˘}}{N}.\end{array}$$

We observe that the approximate stochastic noise has mean and variance given by

$$\begin{array}{}\text{(C23)}& {\displaystyle }\mathrm{E}\left[{\stackrel{\mathrm{˘}}{S}}_i\right]& {\displaystyle }=\mathrm{0},\text{(C24)}& {\displaystyle }\text{Var}\left[{\stackrel{\mathrm{˘}}{S}}_i\right]& {\displaystyle }=C(\mathrm{1}-C)\stackrel{\mathrm{˘}}{N},\end{array}$$

for all i. The initial condition for the approximate model is given by ${\stackrel{\mathrm{̃}}{N}}_i^{\mathrm{0}}=\stackrel{\mathrm{˘}}{N}$ for all i, which has DFT given by

$$\begin{array}{}\text{(C25)}& {\widehat{\stackrel{\mathrm{̃}}{N}}}_k^{\mathrm{0}}=N_{\mathrm{x}}\stackrel{\mathrm{˘}}{N}{\mathit{\delta }}_{k,\mathrm{0}}.\end{array}$$

Appendix D: Recurrence relations for the first and second moments of the approximate model

Our aim is to solve the approximate model Eqs. (C21) and (C22) analytically. Because the process is stochastic, we will solve for the first two moments of the particle counts $\stackrel{\mathrm{̃}}{N}$ in the frequency domain, and in this section we begin by deriving the appropriate recurrence relations.

Taking an expected value of Eq. (C21) gives the following recurrence relation for the first moment:

$$\begin{array}{}\text{(D1)}& {\displaystyle }\mathrm{E}\left[{\widehat{\stackrel{\mathrm{̃}}{N}}}_k^{\mathrm{\ell }+\mathrm{1}}\right]& {\displaystyle }=\mathrm{E}\left[\exp \right(-C{\widehat{d}}_k{\widehat{r}}_k){\widehat{\stackrel{\mathrm{̃}}{N}}}_k^{\mathrm{\ell }}+{\widehat{d}}_k{\widehat{\stackrel{\mathrm{˘}}{S}}}_k]\text{(D2)}& {\displaystyle }& {\displaystyle }=\exp (-C{\widehat{d}}_k{\widehat{r}}_k)\mathrm{E}\left[{\widehat{\stackrel{\mathrm{̃}}{N}}}_k^{\mathrm{\ell }}\right]+{\widehat{d}}_k\mathrm{E}\left[{\widehat{\stackrel{\mathrm{˘}}{S}}}_k\right]\text{(D3)}& {\displaystyle }& {\displaystyle }=\exp (-C{\widehat{d}}_k{\widehat{r}}_k)\mathrm{E}\left[{\widehat{\stackrel{\mathrm{̃}}{N}}}_k^{\mathrm{\ell }}\right],\end{array}$$

where we used the fact that the stochastic noise has zero mean.

Next we obtain a recurrence relation for the second moment of the particle counts. We use Eq. (C21) to compute

$$\begin{array}{ll}{\displaystyle }\mathrm{E}& {\displaystyle }\left[{\widehat{\stackrel{\mathrm{̃}}{N}}}_k^{\mathrm{\ell }+\mathrm{1}}{\widehat{\stackrel{\mathrm{̃}}{N}}}_k^{(\mathrm{\ell }+\mathrm{1})*}\right]\\ {\displaystyle }& {\displaystyle }=\mathrm{E}\left[\right(\exp (-C{\widehat{d}}_k{\widehat{r}}_k){\widehat{\stackrel{\mathrm{̃}}{N}}}_k^{\mathrm{\ell }}+{\widehat{d}}_k{\widehat{\stackrel{\mathrm{̃}}{S}}}_k)\\ \text{(D4)}& {\displaystyle }& {\displaystyle }\times \left(\exp \right(-C{\widehat{d}}_k{\widehat{r}}_k){\widehat{\stackrel{\mathrm{̃}}{N}}}_k^{\mathrm{\ell }}+{\widehat{d}}_k{\widehat{\stackrel{\mathrm{˘}}{S}}}_k{)}^{*}]{\displaystyle }& {\displaystyle }=\exp (-C{\widehat{d}}_k{\widehat{r}}_k-C{\widehat{d}}_k^{*}{\widehat{r}}_k^{*})\mathrm{E}\left[{\widehat{\stackrel{\mathrm{̃}}{N}}}_k^{\mathrm{\ell }}{\widehat{\stackrel{\mathrm{̃}}{N}}}_k^{\mathrm{\ell }*}\right]\\ \text{(D5)}& {\displaystyle }& {\displaystyle }+{\widehat{d}}_k\exp (-C{\widehat{d}}_k^{*}{\widehat{r}}_k^{*})\mathrm{E}\left[{\widehat{\stackrel{\mathrm{̃}}{N}}}_k^{\mathrm{\ell }}{\widehat{\stackrel{\mathrm{˘}}{S}}}_k^{*}\right]\text{(D6)}& {\displaystyle }& {\displaystyle }+\exp (-C{\widehat{d}}_k{\widehat{r}}_k){\widehat{d}}_k\mathrm{E}\left[{\widehat{\stackrel{\mathrm{˘}}{S}}}_k{\widehat{\stackrel{\mathrm{̃}}{N}}}_k^{\mathrm{\ell }*}\right]+{\widehat{d}}_k{\widehat{d}}_k^{*}\mathrm{E}\left[{\widehat{\stackrel{\mathrm{˘}}{S}}}_k{\widehat{\stackrel{\mathrm{˘}}{S}}}_k^{*}\right]{\displaystyle }& {\displaystyle }=\exp (-\mathrm{2}C\text{Re}({\widehat{d}}_k{\widehat{r}}_k\left)\right)\mathrm{E}\left[{\widehat{\stackrel{\mathrm{̃}}{N}}}_k^{\mathrm{\ell }}{\widehat{\stackrel{\mathrm{̃}}{N}}}_k^{\mathrm{\ell }*}\right]\\ \text{(D7)}& {\displaystyle }& {\displaystyle }+{\widehat{d}}_k{\widehat{d}}_k^{*}\mathrm{E}\left[{\widehat{\stackrel{\mathrm{˘}}{S}}}_k{\widehat{\stackrel{\mathrm{˘}}{S}}}_k^{*}\right].\end{array}$$

In the final step above we used the fact that the approximate stochastic noise, ${\widehat{\stackrel{\mathrm{˘}}{S}}}_k$, has zero mean and is uncorrelated with the current solution, ${\widehat{\stackrel{\mathrm{̃}}{N}}}_k$, because the noise is sampled from a fixed distribution, Eq. (C22), at each time step. This means that $\mathrm{E}\left[{\widehat{\stackrel{\mathrm{̃}}{N}}}_k^{\mathrm{\ell }}{\widehat{\stackrel{\mathrm{˘}}{S}}}_k^{*}\right]=\mathrm{E}\left[{\widehat{\stackrel{\mathrm{˘}}{S}}}_k{\widehat{\stackrel{\mathrm{̃}}{N}}}_k^{\mathrm{\ell }*}\right]=\mathrm{0}$, and so the cross terms vanish.

To compute the expected value of the squared magnitude of the stochastic noise, $\mathrm{E}\left[{\widehat{\stackrel{\mathrm{˘}}{S}}}_k{\widehat{\stackrel{\mathrm{˘}}{S}}}_k^{*}\right]$, we use Eq. (F7) and the statistics Eqs. (C23) and (C24) of the stochastic noise to obtain

$$\begin{array}{}\text{(D8)}& {\displaystyle }\mathrm{E}\left[{\widehat{\stackrel{\mathrm{˘}}{S}}}_k{\widehat{\stackrel{\mathrm{˘}}{S}}}_k^{*}\right]& {\displaystyle }=N_{\mathrm{x}}^{\mathrm{2}}\left|\mathrm{E}\right[{\stackrel{\mathrm{˘}}{S}}_{\mathrm{0}}]{|}^{\mathrm{2}}{\mathit{\delta }}_{k,\mathrm{0}}+N_{\mathrm{x}}\text{Var}[{\stackrel{\mathrm{˘}}{S}}_{\mathrm{0}}]\text{(D9)}& {\displaystyle }& {\displaystyle }=N_{\mathrm{x}}C(\mathrm{1}-C)\stackrel{\mathrm{˘}}{N}.\end{array}$$

Substituting this into Eq. (D7) gives the recurrence relation

$$\begin{array}{ll}{\displaystyle }\mathrm{E}\left[\right|{\widehat{\stackrel{\mathrm{̃}}{N}}}_k^{\mathrm{\ell }+\mathrm{1}}{|}^{\mathrm{2}}]& {\displaystyle }=\exp \left(A_k\right)\mathrm{E}\left[\right|{\widehat{\stackrel{\mathrm{̃}}{N}}}_k^{\mathrm{\ell }}{|}^{\mathrm{2}}]\\ \text{(D10)}& {\displaystyle }& {\displaystyle }+\left|{\widehat{d}}_k{|}^{\mathrm{2}}{N}_{\mathrm{x}}C\right(\mathrm{1}-C)\stackrel{\mathrm{˘}}{N},\end{array}$$

where we have also used the amplification factor A_k given by Eq. (A25).

Define the power at wavenumber k by

$$\begin{array}{}\text{(D11)}& {\stackrel{\mathrm{̃}}{P}}_k^{\mathrm{\ell }}=\mathrm{E}\left[\right|{\widehat{\stackrel{\mathrm{̃}}{N}}}_k^{\mathrm{\ell }}{|}^{\mathrm{2}}]\end{array}$$

and the excitation as

$$\begin{array}{}\text{(D12)}& E_k=\left|{\widehat{d}}_k{|}^{\mathrm{2}}{N}_{\mathrm{x}}C\right(\mathrm{1}-C)\stackrel{\mathrm{˘}}{N},\end{array}$$

where we can evaluate

$$\begin{array}{}\text{(D13)}& |{\widehat{d}}_k{|}^{\mathrm{2}}=\mathrm{2}-\mathrm{2}\cos (\mathrm{2}\mathit{\pi }k/N_{\mathrm{x}}).\end{array}$$

Using the above expressions we can write the final recurrence relation for the second moment (the power) as

$$\begin{array}{}\text{(D14)}& {\stackrel{\mathrm{̃}}{P}}_k^{\mathrm{\ell }+\mathrm{1}}=\exp \left(A_k\right){\stackrel{\mathrm{̃}}{P}}_k^{\mathrm{\ell }}+E_k.\end{array}$$

The first term on the right-hand side represents the evolution of the second moment due to the discretized advection scheme, which may preserve the second moment or dissipate it depending on the scheme. This first term is identical to the evolution of the power for the semi-discretization Eq. (A26). The second term on the right-hand side represents a constant injection of variance (energy) due to the stochastic noise.

Appendix E: Analytical solution for the moments of the approximate model

In Appendix D we derived the recurrence relations for the first and second moments of the approximate model. In this section we solve these recurrence relations analytically. Starting with the first moment, the recurrence relation Eq. (D3) has the solution

$$\begin{array}{}\text{(E1)}& {\displaystyle }\mathrm{E}\left[{\widehat{\stackrel{\mathrm{̃}}{N}}}_k^{\mathrm{\ell }}\right]& {\displaystyle }=\exp (-\mathrm{\ell }C{\widehat{d}}_k{\widehat{r}}_k)\mathrm{E}\left[{\widehat{\stackrel{\mathrm{̃}}{N}}}_k^{\mathrm{0}}\right]\text{(E2)}& {\displaystyle }& {\displaystyle }=\exp (-\mathrm{\ell }C{\widehat{d}}_k{\widehat{r}}_k)N_{\mathrm{x}}\stackrel{\mathrm{˘}}{N}{\mathit{\delta }}_{k,\mathrm{0}}\text{(E3)}& {\displaystyle }& {\displaystyle }=\exp (-\mathrm{\ell }C{\widehat{d}}_{\mathrm{0}}{\widehat{r}}_{\mathrm{0}})N_{\mathrm{x}}\stackrel{\mathrm{˘}}{N}{\mathit{\delta }}_{k,\mathrm{0}}\text{(E4)}& {\displaystyle }& {\displaystyle }=N_{\mathrm{x}}\stackrel{\mathrm{˘}}{N}{\mathit{\delta }}_{k,\mathrm{0}}\text{(E5)}& {\displaystyle }& {\displaystyle }=\mathrm{E}\left[{\widehat{\stackrel{\mathrm{̃}}{N}}}_k^{\mathrm{0}}\right],\end{array}$$

where we used the initial condition Eq. (C25) and the fact that ${\widehat{d}}_{\mathrm{0}}=\mathrm{0}$. From this we see that the first moment of the approximate model is constant in time and thus equal to its initial condition. We can write this as

$$\begin{array}{}\text{(E6)}& \mathrm{E}\left[{\stackrel{\mathrm{̃}}{N}}_k^{\mathrm{\ell }}\right]=\stackrel{\mathrm{˘}}{N},\end{array}$$

for all i and ℓ. We thus see that the mean of the approximate model is identical to the solution Eq. (A19) of the deterministic spatial semi-discretization Eq. (A16), which is also constant for a uniform initial condition. That is, the approximate model mean is exactly the same as the exact time integration of the finite-volume discretization, which is consistent with the observation that in Fig. 1 the particle solution oscillates around the finite-volume solution.

To solve the recurrence relation for the second moment Eq. (D14) we first recall that the linear first-order recurrence relation

$$\begin{array}{}\text{(E7)}& z^{\mathrm{\ell }+\mathrm{1}}=a{z}^{\mathrm{\ell }}+b\end{array}$$

for $a\in [\mathrm{0},\mathrm{1}]$ has the solution

$$\begin{array}{}\text{(E8)}& z^{\mathrm{\ell }}=\left\{\begin{array}{ll}{a}^{\mathrm{\ell }}{z}^{\mathrm{0}}+(\mathrm{1}-a^{\mathrm{\ell }})z^{\mathrm{\infty }}& \text{ if }a<\mathrm{1},\\ z^{\mathrm{0}}+\mathrm{\ell }b& \text{ if }a=\mathrm{1},\end{array}\right.\end{array}$$

where z⁰ is the initial condition and z^∞ is the steady-state solution in the decaying case, given by

$$\begin{array}{}\text{(E9)}& z^{\mathrm{\infty }}={\displaystyle \frac{b}{\mathrm{1}-a}}.\end{array}$$

Applying this to Eq. (D14) gives

$$\begin{array}{}\text{(E10)}& {\stackrel{\mathrm{̃}}{P}}_k^{\mathrm{\ell }}=\left\{\begin{array}{ll}\exp \left(\mathrm{\ell }{A}_k\right){\stackrel{\mathrm{̃}}{P}}_k^{\mathrm{0}}+(\mathrm{1}-\exp (\mathrm{\ell }{A}_k\left)\right){\stackrel{\mathrm{̃}}{P}}_k^{\mathrm{\infty }}& \\ & \mathrm{\text{if }}{A}_k<\mathrm{0},\\ {\stackrel{\mathrm{̃}}{P}}_k^{\mathrm{0}}+\mathrm{\ell }{E}_k& \mathrm{\text{if }}{A}_k=\mathrm{0},\end{array}\right.\end{array}$$

where the limiting moments are

$$\begin{array}{}\text{(E11)}& {\displaystyle }{\stackrel{\mathrm{̃}}{P}}_k^{\mathrm{0}}& {\displaystyle }=N_{\mathrm{x}}^{\mathrm{2}}{\stackrel{\mathrm{˘}}{N}}^{\mathrm{2}}{\mathit{\delta }}_{k,\mathrm{0}},\text{(E12)}& {\displaystyle }{\stackrel{\mathrm{̃}}{P}}_k^{\mathrm{\infty }}& {\displaystyle }={\displaystyle \frac{E_k}{\mathrm{1}-\exp \left(A_k\right)}},\end{array}$$

using Eq. (C25). To evaluate the above expression we need the amplification factor Eqs. (B8)–(B13), the excitation Eq. (D12), and the squared magnitude Eq. (D13).

Appendix F: Power identity for vectors with i.i.d. random components

Consider a vector of independent and identically distributed (i.i.d.) random variables z_i for $i=\mathrm{0},\mathrm{\dots },N_{\mathrm{x}}-\mathrm{1}$. We want to compute the expected value of the squared magnitude of the DFT of this vector, i.e., $\mathrm{E}\left[\right|{\widehat{z}}_k{|}^{\mathrm{2}}]$ for each wavenumber k.

We start by observing that $\mathrm{E}\left[z_i{z}_j^{*}\right]=\mathrm{E}\left[z_i\right]\mathrm{E}\left[z_j\right]=\left|\mathrm{E}\right[z_{\mathrm{0}}]{|}^{\mathrm{2}}$ for i≠j because the random variables are independent. We also have $\mathrm{E}\left[z_i{z}_i^{*}\right]=\mathrm{E}\left[\right|z_i{|}^{\mathrm{2}}]=\mathrm{E}[\left|z_{\mathrm{0}}{|}^{\mathrm{2}}\right]=\text{Var}\left[z_{\mathrm{0}}\right]+\left|\mathrm{E}\right[z_{\mathrm{0}}]{|}^{\mathrm{2}}$. We can thus write

$$\begin{array}{}\text{(F1)}& \mathrm{E}\left[z_i{z}_j^{*}\right]=\left|\mathrm{E}\right[z_{\mathrm{0}}]{|}^{\mathrm{2}}+\text{Var}[z_{\mathrm{0}}]{\mathit{\delta }}_{i,j},\end{array}$$

for all i,j, where δ_i,j is the Kronecker delta.

We can now compute the expected value of the squared magnitude of the DFT:

$$\begin{array}{ll}\text{(F2)}& {\displaystyle }\mathrm{E}\left[\right|{\widehat{z}}_k{|}^{\mathrm{2}}]& {\displaystyle }=\mathrm{E}\left[{\widehat{z}}_k{\widehat{z}}_k^{*}\right]{\displaystyle }& {\displaystyle }=\mathrm{E}\left[\left(\sum _{j=\mathrm{0}}^{N_{\mathrm{x}}-\mathrm{1}}{z}_j\exp (-i\mathrm{2}\mathit{\pi }jk/N_{\mathrm{x}})\right)\right.\\ \text{(F3)}& {\displaystyle }& {\displaystyle }\times \left.\left(\sum _{\mathrm{\ell }=\mathrm{0}}^{N_{\mathrm{x}}-\mathrm{1}}{z}_{\mathrm{\ell }}^{*}\exp (i\mathrm{2}\mathit{\pi }\mathrm{\ell }k/N_{\mathrm{x}})\right)\right]\text{(F4)}& {\displaystyle }& {\displaystyle }=\sum _{j=\mathrm{0}}^{N_{\mathrm{x}}-\mathrm{1}}\sum _{\mathrm{\ell }=\mathrm{0}}^{N_{\mathrm{x}}-\mathrm{1}}\mathrm{E}\left[z_j{z}_{\mathrm{\ell }}^{*}\right]\exp (-i\mathrm{2}\mathit{\pi }(j-\mathrm{\ell })k/N_{\mathrm{x}}){\displaystyle }& {\displaystyle }=\sum _{j=\mathrm{0}}^{N_{\mathrm{x}}-\mathrm{1}}\sum _{\mathrm{\ell }=\mathrm{0}}^{N_{\mathrm{x}}-\mathrm{1}}\left(\right|\mathrm{E}\left[z_{\mathrm{0}}\right]{|}^{\mathrm{2}}+\text{Var}\left[z_{\mathrm{0}}\right]{\mathit{\delta }}_{j,\mathrm{\ell }})\\ \text{(F5)}& {\displaystyle }& {\displaystyle }\times \exp (-i\mathrm{2}\mathit{\pi }(j-\mathrm{\ell })k/N_{\mathrm{x}}){\displaystyle }& {\displaystyle }=\left|\mathrm{E}\right[z_{\mathrm{0}}\left]{|}^{\mathrm{2}}\sum _{j=\mathrm{0}}^{N_{\mathrm{x}}-\mathrm{1}}\sum _{\mathrm{\ell }=\mathrm{0}}^{N_{\mathrm{x}}-\mathrm{1}}\exp \right(-i\mathrm{2}\mathit{\pi }(j-\mathrm{\ell })k/N_{\mathrm{x}})\\ {\displaystyle }& {\displaystyle }+\text{Var}\left[z_{\mathrm{0}}\right]\sum _{j=\mathrm{0}}^{N_{\mathrm{x}}-\mathrm{1}}\sum _{\mathrm{\ell }=\mathrm{0}}^{N_{\mathrm{x}}-\mathrm{1}}{\mathit{\delta }}_{j,\mathrm{\ell }}\\ \text{(F6)}& {\displaystyle }& {\displaystyle }\times \exp (-i\mathrm{2}\mathit{\pi }(j-\mathrm{\ell })k/N_{\mathrm{x}}).\end{array}$$

Consider the first term in the above expression. When k=0 the sum is $N_{\mathrm{x}}^{\mathrm{2}}$ and when k≠0 the sum is zero because the inner sum consists of N_x complex numbers that are spaced around the unit circle in a symmetric fashion. Now consider the second term. This collapses to ${\sum }_{j=\mathrm{0}}^{N_{\mathrm{x}}-\mathrm{1}}\exp (-i\mathrm{2}\mathit{\pi }(j-j)k/N_{\mathrm{x}})=N_{\mathrm{x}}$ for all k. We thus have the final expression

$$\begin{array}{}\text{(F7)}& \mathrm{E}\left[\right|{\widehat{z}}_k{|}^{\mathrm{2}}]=N_{\mathrm{x}}^{\mathrm{2}}|\mathrm{E}\left[z_{\mathrm{0}}\right]{|}^{\mathrm{2}}{\mathit{\delta }}_{k,\mathrm{0}}+N_{\mathrm{x}}\text{Var}\left[z_{\mathrm{0}}\right].\end{array}$$

We see that the power spectrum consists of a uniform component that depends on the variance of the random variable and a DC component that depends on the mean of the random variable.

Appendix G: Symbols used in this paper

Table G1 lists the symbols used in this paper.

Table G1 Symbols used in this paper.