Quantitative evaluation of numerical integration schemes for Lagrangian particle dispersion models

A rigorous methodology for the evaluation of integration schemes for Lagrangian particle dispersion models (LPDMs) is presented. A series of one-dimensional test problems are introduced, for which the Fokker-Planck equation is solved numerically using a finite-difference discretisation in physical space, and a Hermite function expansion in velocity space. Numerical convergence errors in the Fokker-Planck equation solutions are shown to be much less than the statistical error associated with 5 a practical-sized ensemble (N = 10) of LPDM solutions, hence the former can be used to validate the latter. The test problems are then used to evaluate commonly used LPDM integration schemes. The results allow for optimal time-step selection for each scheme, given a required level of accuracy. The following recommendations are made for use in operational models. First, if computational constraints require the use of moderate to long time-steps it is more accurate to solve the random displacement model approximation to the LPDM, rather than use existing schemes designed for long time-steps. Second, useful gains in 10 numerical accuracy can be obtained, at moderate additional computational cost, by using the relatively simple ‘small-noise’ scheme of Honeycutt.


Introduction
State-of-the-art Lagrangian particle dispersion models (LPDMs hereafter), for example FLEXPART (Stohl et al., 2005) and NAME (Jones et al., 2007), are key scientific tools for the study of the long-range transport and dispersal of the transport of atmospheric trace gases and aerosols.Applications are diverse, e.g.establishing the relationship between emissions of pollutants and air quality downstream Cassiani et al. (2013), aerosol dispersal following volcanic eruptions (Devenish et al., 2011;D'Amours et al., 2010), modelling of nuclear accident scenarios (Stohl et al., 2012), and determination of constraints on chemical emissions via inverse modelling (Seibert and Frank, 2004;Stohl et al., 2010).More fundamentally, LPDMs can be used to address key scientific questions concerning the nature of transport in the atmosphere (Legras et al., 2005;Berthet et al., 2007), including how transport might be influenced by a changing climate.
Mathematically speaking, LPDMs are formulated as stochastic differential equations (SDEs hereafter).(It is notable that it is possible to include jump processes (Platen and Liberati, 2010) as a representation of non-local convective parameterisations (Forster et al., 2007), but we will not be concerned here with this possibility.)Although the numerical analysis of solution techniques for SDEs (e.g.Kloeden and Platen, 1992;Milstein and Tretyakov, 2004) is a mature subject in mathematics, LPDMs have not, generally speaking, exploited developments in the subject, and are typically formulated using numerical schemes adapted from those used for ordinary differential equations (see e.g.Stohl et al., 2005).Validation of LPDMs has focussed instead on direct comparison with observational data (Stohl et al., 1998;Ryall and Maryon, 1998).Our contention is that observational comparison, while clearly a necessary aspect of model development, will be insufficient if any uncertainty exists concerning the accuracy of the numerical solution of the underlying equations.
The aim of the present work, therefore, is to introduce a rigorous framework for the testing and evaluation of numerical schemes for LPDMs.The framework is based on a standard one-dimensional dispersion model problem Published by Copernicus Publications on behalf of the European Geosciences Union.(Rodean, 1996;Wilson and Sawford, 1996) modelling the vertical dispersion of air parcels in the atmospheric boundary layer (ABL hereafter).Vertical profiles of turbulent statistics representative of both stable and neutral conditions will be considered, and the LPDM equations will be of the "wellmixed" class (Thomson, 1987), meaning that long time probability distribution of the solutions (the invariant measure of the SDEs) is given by a pre-specified "atmospheric" distribution (taken here to be uniform in physical space and Gaussian in velocity space).Hence the model problem, while idealised, captures key elements of the physics of dispersion in the stable and neutral ABL.The convective case, in which the vertical velocity statistics are non-Gaussian (see e.g.Cassiani et al., 2015), will require a separate test case and is not discussed here.
Our approach to evaluating a given LPDM numerical scheme is to cross-validate its performance against a numerical solution of the corresponding Fokker-Planck equation (FPE hereafter; see e.g.Gardiner, 2009).The FPE describes the time evolution of the probability density function (pdf) of the stochastic process, and is formulated in position-velocity space, and so in the context of the current problem of dispersion in one spatial dimension is a partial differential equation in 2 + 1 dimensions.Note that in three spatial dimensions in which the FPE is a 6 + 1 dimensional PDE, it will be computationally impractical in most circumstances to obtain accurate solutions to the FPE, and consequently LPDMs will be the only practical tool to solve the problem.
A solution method based on a Hermite function expansion is introduced in order to obtain accurate solutions of the FPE with computational efficiency.Evaluation of the LPDM scheme proceeds by a comparison of pdfs in appropriate error norm, where the LPDM pdf is generated from an ensemble of solutions, using the kernel density method (e.g.Silverman, 1986;Wand and Jones, 1994).The performances of various schemes are evaluated, as a function of time step t, including the textbook (basic) Euler-Maruyama scheme, the second-order and third-order weak Runge-Kutta scheme of Platen (see Sect. 15.1 of Kloeden and Platen, 1992), the "small-noise" second-order Runge-Kutta method of Honeycutt (Honeycutt, 1992), the "long time-step" scheme used operationally in FLEXPART (Stohl et al., 2005), and a suggested improvement to this last scheme.
The outline of the work is as follows.In Sect.2, the SDEs describing the evolution of particle trajectories in the LPDM are introduced, together with the corresponding FPE.A numerical solution scheme for the FPE is described and solutions are obtained and benchmarked for a number of test cases.In Sect.3, the methodology for using the FPE solution to assess specific numerical schemes for the LPDM is presented, and in Sect. 4 this methodology is then applied to specific schemes discussed above.In Sect. 5 the consequences of our findings are discussed and conclusions are drawn.

The model problem formulated as an LPDM
Consider a horizontally homogeneous turbulent ABL of uniform density, with a vertical velocity distribution that is Gaussian with zero mean and standard deviation σ w (z), and which has a Lagrangian decorrelation timescale τ (z).The canonical stochastic differential equation model (e.g.Rodean, 1996;Wilson and Sawford, 1996) for onedimensional vertical dispersion in the ABL is Here W t and Z t are the vertical velocity and height of a given air parcel.Both are stochastic variables, with each individual realisation determined by that of the Brownian (or Wiener) process B t .Further σ w = σ w (Z t ) and τ = τ (Z t ) are the values of σ w (z) and τ (z) local to the parcel.In operational LPDMs, such as FLEXPART, appropriate vertical profiles for σ w (z) and τ (z) are specified based on empirical fits to observations of different ABL conditions, as will be discussed below.The equation set ( 1) is typically augmented with reflecting boundary conditions at the Earth's surface and at the ABL top (see Thomson et al., 1997, for a detailed discussion of the top boundary condition).For definiteness, for our test-case runs, the initial velocity for Eq.(1) at t = 0 is sampled from a normal distribution W 0 ∼ N (0, σ 2 w (z 0 )) and, for ease of comparison to the FPE results below, the initial position is sampled from a distribution Z 0 ∼ N (z 0 , σ 2 z ) centred on an initial height z 0 with standard deviation σ z .
For the purposes of numerical solution, it is more convenient (e.g.Sect.3.1 of Rodean, 1996) to use Ito's lemma to express Eq.(1) in terms of the variables t = W t /σ w (Z t ) and Z t , leading to The simpler form (2) is exactly equivalent to Eq. (1).Moreover, the FPE of Eq. (2) has a considerably simpler form than the corresponding FPE of Eq. (1), a fact which will prove useful below.
It is simplest to view Eq. ( 2) as a non-dimensional equation, given that in particular t is already a non-dimensional variable.Hanna, 1982).The non-dimensional parameter = u * /f h is a boundary layer Rossby number (the value = 0.8 is taken in the test case).For the purposes of numerical stability (see text), in practice the modified profiles σw (z) and τ (z) are used, where is chosen to avoid singular behaviour at the boundaries (z b = 0.05).
To specify our test-case problems it is necessary to choose suitable (non-dimensional) profiles for σ w (z) and τ (z).Here we choose to focus on three such profiles, two of which are widely used (Hanna, 1982;Stohl et al., 2005) empirical fits to observed statistics in stable and neutral conditions respectively.The third has τ (z) constant and a linear profile for σ w (z), and is used to demonstrate a new LPDM scheme introduced below.The details of the profiles used are given in Table 1 and are plotted in Fig. 1.In practice, the exact profiles suggested by Hanna (1982) are modified slightly, as detailed in Table 1, to avoid singular behaviour at the ABL top and bottom.This is necessary because in Hanna's original profiles either σ w → 0 or τ → 0 as z → 0, 1, with neither type of behaviour being physical.
In Sect. 4 large ensembles of numerical solutions of Eq. (2) will be calculated using different numerical integration schemes.The accuracy of each numerical scheme, as a function of time step t, will be assessed by comparison with the corresponding solution of the FPE, to be detailed next.

The model problem formulated as an FPE
Following the standard procedure in stochastic calculus (e.g.Sect. 3.4.1 of Gardiner, 2009), the FPE which describes the time evolution of the probability density p(ω, z, t) of ( t , Z t ) in Eq. ( 2) can be obtained as Explicitly, here ω = w/σ w .The initial conditions consistent with those given in Eq. ( 2) are (for σ z 1 and z 0 not near the boundaries) The FPE Eq. (3) also requires boundary conditions at z = 0, 1 which are consistent with the reflecting boundary conditions for the LPDM.The boundary conditions consistent with reflection are which in probabilistic terms is equivalent to the reflection condition t → − t being applied at the boundaries.Wilson et al. (1993) found that this perfect reflection algorithm is exactly consistent with the "well-mixed constraint" in homogeneous Gaussian turbulence (see also the appendix of Sect.11 of Rodean, 1996).Equations ( 3)-( 5) constitute a well-defined initial-value problem which is suitable for numerical solution.An important quantity obtained from the solution p(ω, z, t) is the physical concentration of parcels given by (In general, tracer concentrations and the marginal probability given in Eq. 6 can differ by a normalisation constant.)The concentration c(z, t) will be our main benchmark quantity in Sect. 4 below.
3 Numerical solution of the FPE

The Hermite expansion for the FPE
The non-dimensionalised FPE Eq. ( 3) is a hypo-elliptic differential equation defined on R × [0, 1].Our approach to its numerical solution is to seek a solution based on the following Hermite polynomial expansion: In statistics, this expansion is also known as the Gram-Charlier series of Type A (see p. 23 of Barndorff-Nielsen and Cox, 1989).Here the functions C k (z, t) denote the projection, at the vertical level and time (z, t), of p(ω, z, t) onto the (probabilists') Hermite function He k (ω)e −ω 2 /2 / √ 2π where He k (ω) is the Hermite polynomial defined by Notice that it follows that the particle concentration Eq. ( 6) satisfies c(z, t) = C 0 (z, t).
Before inserting the expansion Eq. ( 7) into the FPE Eq. ( 3) it is helpful to rewrite the FPE in the form In this form the Hermite function identity Eq. (A2) can be used to evaluate the first term on the right-hand side.Further, the second and third terms on the right-hand side can be simplified using the derivative and recursion formulae for Hermite polynomials (Eqs.A5, A6).After some working the result is (using the convention Using the orthogonality property of Hermite functions (A3) it follows that The system (11) constitutes an infinite system of coupled 1 + 1 dimensional first-order partial differential equations for the coefficients C k .For a numerical solution this series can be truncated as we describe below.The initial conditions for Eq. ( 11) are easily obtained from (4) using the orthogonality property, The boundary conditions can be obtained using the symmetry He k (ω) = (−1) k He k (−ω).Substituting the expansion (7) into the boundary condition (5), it follows that and consequently It may seem surprising that the even equations have no boundary condition and the odd equations take two boundary conditions.However, as the system (11) consists of firstorder PDEs it is clear that the total number of boundary conditions will be correct, provided that the series is truncated at k = K odd.
It is worth noting that the series (11) can also be truncated at K = 0 by using an (approximate) quasi-stationary balance in the k = 1 equation of the form which results in the diffusion equation Geosci.Model Dev., 9, 2441-2457, 2016 www.geosci-model-dev.net/9/2441/2016/It is well known (e.g.Sect.3.5 of Thomson, 1987) that the LPDM Eq. ( 1) can be approximated by a random walk ("random displacement" or RDM) model Equation ( 16) is simply the Fokker-Planck equation of the RDM model (Eq.17), with the diffusivity κ of the RDM being κ = σ 2 w τ .Note that the RDM model can be derived formally from the LPDM in the distinguished limit of a short decorrelation time, σ w → ∞, τ → 0 with σ 2 w τ = κ finite (see e.g.Sect.6.3 of Rodean, 1996).It is much easier to obtain accurate solutions of Eq. ( 17), compared to Eq. ( 1) at relatively large time steps; hence, an interesting question concerns when exactly it is preferable to solve Eq. ( 17) rather than Eq. ( 1).This question is best answered by quantifying the difference between the solution of Eqs. ( 16) and ( 11) and using this difference as a benchmark for assessing the errors in LPDM calculations, as will be done in Sect. 4 below.

The numerical method and benchmark solutions for the FPE
Based on the analysis above, Eq. ( 3) can be solved numerically by integrating the system (11) with boundary conditions (Eq.14), truncated at k = K odd.Our approach is to use a standard finite-difference discretisation with N z grid points, equally spaced with z = 1/N z , on a staggered cell-centred grid (i.e.z i = (i − 1/2) z, for i = 1, . .., N z ) in order to apply the boundary conditions at z = 0, 1 systematically.The details of the implementation of the boundary conditions are described in Appendix B. The set ( 11) is stiff and a naive solution method would have the time step t bounded above by t Min z τ (z)/K, i.e. the timescale of exponential decay of the highest Hermite function mode.However, considerably longer time steps can be used if an exponential time-stepping scheme is chosen.Our choice is the Exponential Time-Differencing fourthorder Runge-Kutta (ETDRK4) scheme of Kassam and Trefethen (2005), with the "linear" operator in that scheme taken to be the first term on the right-hand side of Eq. ( 11) only, because it is this first term that is responsible for the stiffness of Eq. ( 11).
To obtain our benchmark solutions of Eq. ( 11) and therefore Eq. (3), tests of the convergence of the solutions as both t and z are decreased and K is increased have been performed.For all three case studies, it was found to be adequate to take K = 19 to obtain fully converged solutions, because the Hermite series was found to converge rapidly, i.e. |C 19 | 10 −16 everywhere in the domain.Comparison of a sequence of solutions with z = 1/N z with N z = 2 7 , 2 8 , . .., 2 12 revealed quadratic convergence with z as expected for our scheme.Figure 2 shows the relative error E j (t), with reference to the next-highest resolution solution, in the L 2 norm for the mean concentration c(z, t) at fixed times, for the two test cases.That is, where C j 0 (z, t) denotes the solution with N z = 2 j .Quadratic convergence is evident from the slope of the graphs in Fig. 2. For example, typical numerical errors at N z = 2 12 (highest resolution) are E 12 (t 1 ) = 9.7 × 10 −5 (stable boundary layer at t 1 = h/u * ) and 1.3 × 10 −4 (neutral boundary layer at t 1 = 3 h/u * ) respectively.The numerical accuracy above is sufficient for benchmarking our LPDM solutions, because the statistical error associated with reasonable-sized ensembles (N = 10 6 ) of the LPDM is of order E(t 1 ) ≈ 10 −2 , as will be discussed below.
Figure 3 shows snapshots of the particle concentration c(z, t) for each of the three FPE benchmark solutions described above.The left panel shows the constant τ case, the middle panel shows the stable ABL case, and the right panel the neutral ABL.In all three cases particles are initialised close to z = z 0 = 1/2 and disperse to become well mixed throughout the ABL at late times.The neutral and stable cases differ in that mixing is rather more rapid (in terms of the dimensional timescale h/u * ) for the stable case compared to the neutral case.Also, in the neutral case, mixing is relatively slow towards the top of the ABL where the amplitude of turbulent fluctuations decays exponentially.

Evaluation of numerical schemes for LPDMs
In this section, a range of textbook, commonly used and new numerical schemes for LPDMs will first be introduced and then evaluated using the FPE solutions described above as a benchmark.The task is somewhat simplified because the equation set ( 2) is time-independent (autonomous).Note that it may be necessary to modify some of the schemes described below if an ABL with time-dependent statistics is to be modelled with the same formal accuracy.Note that, in the terminology of SDE numerical schemes (Kloeden and Platen, 1992), we are able to use "weak" schemes (convergent in probability) in addition to "strong" schemes (convergent in path), because we are primarily interested in the concentration of particles, which can be obtained from the pdf p(ω, z, t).The rate of convergence of a scheme, as measured by quantities which depend on the pdf such as the concentration c(z, t), with respect to the time step t is known as its "weak" order (see e.g.Chap.9 of Kloeden and Platen, 1992).The weak order is the relevant measure of comparison between schemes for our study, and should not be confused with the "strong" order of a scheme, which refers to the rate of convergence of solution paths with respect to specific stochastic realisations.
It is important to note, however, that it is by no means obvious that a given scheme will attain its formal weak order when solving Eq. ( 2).This is because the assumptions under which the weak order of each scheme is derived are not met in the case of Eq. (2) because of the reflection boundary conditions.It is therefore necessary to solve Eq. (2) explicitly to assess each scheme.

LPDM numerical schemes
Tables 2-3 summarise the SDE numerical schemes to be investigated.The first, most obvious scheme to test is the Euler-Maruyama (E-M) scheme (Maruyama, 1955), i.e. the simplest and lowest order time-stepping scheme for SDEs.Next, as with ordinary differential equations (ODEs), it is possible to construct schemes with higher orders of formal accuracy in the spirit of Runge-Kutta schemes for ODEs.
Here we test the performance of Platen's "explicit order 2.0 weak scheme" (EXPLICIT 2.0) and "explicit order 3.0 weak scheme" (EXPLICIT 3.0) (see Chap. 15 of Kloeden and Platen, 1992).In common with schemes for ODEs, higher order schemes are somewhat more complicated to implement, and are more computationally expensive per time step t.The advantage, however, is that the schemes have weak order t 2 (EXPLICIT 2.0) and t 3 (EXPLICIT 3.0) compared to t for E-M.
A single candidate from a second class of schemes, the socalled "small-noise" schemes, to be investigated is the HON-SRKII scheme of Honeycutt (1992).Small-noise schemes typically have the same weak order ( t) as E-M (see e.g. the discussion in Chap. 3 of Milstein and Tretyakov, 2004), but the schemes are designed so that the leading-order error depends on the "noise amplitude" in the equation, which in many practical situations is sufficiently small that higherorder convergence is observed in practice (at least for intermediate length time steps; see discussion below).The HON-SRKII scheme will be shown below to converge with global error ∼ t 2 in this intermediate time-step regime.
A third class of schemes to be investigated is designed to work well with long time steps.Such schemes are of interest operationally, because the practical advantages of calculating large ensembles efficiently are thought to outweigh the disadvantage of loss of accuracy due to time-stepping errors.The FLEXPART model (Stohl et al., 2005), for example, switches between using E-M and a long time-stepping scheme due to Legg and Raupach (1982, LEGGRAUP).It is of some interest to verify that long time-stepping schemes such as LEG-GRAUP do indeed outperform E-M at operationally relevant values of t.In fact, in Appendix C we review the derivation of the LEGGRAUP scheme, and show that additional care is Table 2.The LPDM numerical schemes investigated in Sect. 4.Here t is the time step, B n ∼ N (0, t), n ∼ N (0, 1) and σ i = σ w (Z i ), and τ n = τ (Z n ).The drift function is denoted by F i = − i /τ (Z i ) + σ w (Z i ) where i = n, µ.

Scheme
Algorithm Reference and notes Kloeden and Platen (1992) Legg and Raupach ( 1982) needed to couple the velocity and position equations.A corrected scheme (LONGSTEP) is derived in Appendix C and is then compared with the schemes listed above in Sect.4.2.
The method used to compare the results from a particular scheme, at fixed time step t, to the particle concentration c(z, t) obtained from the numerical solution of the FPE, is as follows.First, a large ensemble (typically N = 10 6 ) of trajectories is calculated using the scheme under investigation.Next, the density of particles c is reconstructed from the resulting ensemble {Z (i) t , i = 1, . .., N } using kernel density estimation Here h b > 0 is a (small) smoothing parameter known as the bandwidth, and "image terms" refer to contributions from the images of trajectories, introduced to satisfy the boundary conditions.The function K(•) is the kernel function, and is non-negative with zero mean and has unit integral.Here we use a Gaussian kernel.Details, including how the optimal bandwidth h b = h * is chosen in practice, are given in Appendix D.

Scheme Algorithm
Reference and notes Kloeden and Platen (1992) where φ = ζ, ρ and . L 2 error (Eq.20) as a function of non-dimensional time step t u * / h for the constant τ = 0.1 test case with N = 10 6 ensemble integrated at time t = h/u * .The LONGSTEP scheme (purple diamonds) gives the best results in this case.Blue lines of slopes 1, 2, and 3 are plotted for reference.
resentative computational costs, in Table 4 the relative cost measured with reference to the E-M scheme is shown for our calculations.Following best practice in large operational calculations (see e.g.Stohl et al., 2005), the random numbers used to simulate the Wiener processes are pre-calculated so the costs of their generation are not included in the comparison.Another possible computational saving comes from the use of variable time steps.To test whether or not a significant computational saving is easily attainable, we have made some calculations in which t ∝ τ (the local Lagrangian decorrelation time).For each scheme tested, the use of variable time steps was found to lead to a computational saving of a factor of around 2 to 3 compared to fixed time steps, with the schemes otherwise performing as detailed below.More details on variable time-stepping schemes will be given elsewhere.

Results
The main results, showing the performance of the six schemes described in Tables 2-3 over a wide range of time steps t, are shown in Figs.4-6.Figures 4-6 detail the results for the constant τ test case, the stable ABL test case and the neutral ABL test case respectively (see Table 1).In each Table 4. Computational clock times, measured relative to the E-M scheme, for all of the schemes detailed in Tables 2 and 3.The calculations are for N = 10 6 trajectories, with time step t = 10 −3 h/u * and integration time h/u * .The computational times are obtained by taking the average of times elapsed in seconds from several simulations, coded in MATLAB, on a MacBook Pro with no other programs running.

Scheme
Relative computational time E-M 1.0 EXPLICIT 2.0 2.0 EXPLICIT 3.0 5.8 HON-SRKII 1.9 LEGGRAUP 1.2 LONGSTEP 1.5 figure, the L 2 error (Eq.20) is plotted as a function of nondimensional time step t u * / h.Logarithmic scales are used so that lines of constant slope correspond to the observed order of the schemes.Blue lines with slopes 1, 2, and 3 are plotted for reference.The statistical error, which is the lowest possible error that can be measured for a given scheme, is plotted as a solid black line in each panel.Also plotted in Figs.4-6, as a dotted black line, is the L 2norm difference c − C 0 2 between the concentration field c(z, t) obtained from the solution of the FPE Eq. ( 3) and C 0 (z, t) obtained from the diffusion Eq. ( 16).The dotted black line marks an important boundary on each panel.If the time step t is such that the error of a given scheme lies above this line, then it is preferable to solve the RDM Eq. ( 17) in place of Eq. (2), because (at fixed t) the numerical error for the former is more easily controlled.
Figure 4 shows results for the constant τ test case at time t = h/u * (see Fig. 3 and Table 1 for details).The lowest order schemes, LEGGRAUP (blue circles) and E-M (black squares), are seen to realise their formal weak order t.EX-PLICIT 2.0 (red hexagons) and HON-SRKII (green solid triangles) have weak order t 2 , whereas EXPLICIT 3.0 (blue triangles) has weak order t 3 , as expected.The best performing scheme for this particular case is the new scheme LONGSTEP (purple diamonds) derived in Appendix C. The rationale for LONGSTEP is that there is a conceptual error in the derivation of LEGGRAUP, which results in its performance at large t being no better than E-M.When this error is corrected in LONGSTEP, the performance is better than even EXPLICIT 3.0.LONGSTEP in effect uses exact solutions of the LPDM equations for constant τ and linear σ w , meaning that if the same calculations had been performed in an infinite domain, the numerical error would be zero.In the constant τ test case, errors in LONGSTEP arise only from the reflection boundary conditions at z = 0, 1. How- ever, LONGSTEP does not fare well in the remaining two test cases to be described next.
Figure 5 shows results for the stable ABL test case at intermediate time t = h/u * (upper panel) and at late time t = 4 h/u * (lower panel), when the concentration is almost well mixed across the ABL (see Fig. 3).The results are similar to those of the constant τ case, except LONGSTEP (purple diamonds) now performs as poorly as E-M.Both E-M and LONGSTEP outperform LEGGRAUP.It was not found to be possible to obtain solutions for EXPLICIT 3.0 using time steps longer than t = 0.02 h/u * because of problems with the reflective boundary conditions.
Figure 6 shows the results for the neutral ABL case at intermediate time t = 3 h/u * (upper panel) and at late time t = 12 h/u * (lower panel).In this case the performance of LONGSTEP and LEGGRAUP are comparable, but with the E-M scheme performing better than both, except at very long time steps where LEGGRAUP has slightly better accuracy at long time steps.As for the previous test cases, the EXPLICIT 3.0 (blue triangles) scheme gives the lowest errors (weak order t 3 ), and EXPLICIT 2.0 (red hexagons) along with HON-SRKII (green solid triangles) perform consistently well with weak order t 2 .
To give an impression of where the particle concentration errors are accumulating, Fig. 7 shows snapshots of particle density c(z, t) for the stable ABL case, at t = h/u * .
Results are shown for each scheme when a long time step t = 0.05 h/u * is used (left panel) and a moderate time step t = 0.007 h/u * (right panel).The errors in the long timestep case are large and are largely due to issues with the reflection of trajectories at the surface (z = 0).Numerical accuracy requires that t τ , which is evidently violated close to the boundary where τ (z) is small (see Fig. 1).Errors due to reflection are particularly acute for the higher order schemes (such as EXPLICIT 2.0 and HON-SRKII) that require the treatment of an intermediate step(s).See the discussion in Appendix B for how this step is implemented.The stable boundary layer case, where τ decays most rapidly near the z = 0 boundary, is the case which appears to be the most sensitive to the treatment of reflection there.

Conclusions
The main contribution of this paper is to introduce a protocol for the quantitative assessment of SDE numerical schemes, applied to the problem of dispersion in an idealised atmospheric boundary layer, as modelled by LPDMs.Accurate solutions of the Fokker-Planck equation (FPE, Eq. 3) are used to benchmark the distribution obtained from an ensemble of LPDM solutions obtained using a particular scheme with a fixed time step t.By using the FPE solution, our pro- tocol avoids the possibility of the LPDMs exhibiting spurious convergence to an incorrect distribution as t → 0 (e.g. by a poor treatment of reflection boundary conditions), and the FPE provides independent verification of the correctness of a specific implementation.
The convergence results obtained in our model test problems are valuable because, due to the importance of reflection of particles from the surface and top of the boundary layer, it is not possible to rely on the formal convergence rates of different SDE schemes (as given by e.g.Kloeden and Platen, 1992).All of the schemes tested attain their formal convergence rates at early times in the model test problem, i.e. before reflection becomes important, and thereafter are limited to an extent by the details of how reflection is implemented (see Appendix B for discussion).
Our results allow the following recommendations to be made, for consideration by operational modellers.
-For our test problems, the best results with respect to accuracy as a function of t were obtained with the EXPLICIT 3.0 weak order t 3 scheme.However, this scheme is time-consuming to implement and more expensive per step compared to the weak order t 2 schemes investigated, so the gains associated with it are marginal.A good compromise between ease-ofimplementation, flexibility and accuracy is the "smallnoise" scheme of Honeycutt (1992, here HON-SRKII).Formally, the weak order of HON-SRKII is just t, i.e. the same as Euler-Maruyama.However, the scheme designed so that at fixed t, in the limit of small noise, the weak error scales with t 2 (e.g.Chap. 3 of Milstein and Tretyakov, 2004).Although the boundary layer dispersion problems examined here are not formally "smallnoise" problems, our results show clearly that they behave as such in a practical implementation.As a consequence HON-SRKII scheme performs at least as well as the formally weak order t 2 scheme EXPLICIT 2.0 (which in fact has a very similar implementation for the specific LPDM problem we have examined here).
-The "long-step" scheme due to Legg and Raupach (1982, here LEGGRAUP), which is used operationally for global integrations of trajectories in FLEXPART (for example), should be avoided.LEGGRAUP does not significantly outperform Euler-Maruyama at any time step for any of the three profiles we have studied.The reason for this is a conceptual error in its derivation, which we have corrected here in the development of a new scheme -LONGSTEP -see C. LONGSTEP performs very well in the case of τ (z) = constant, but no better than LEGGRAUP for other τ (z) profiles; hence, we do not recommend it for operational use either.
-Global calculations often require the use of long time steps for reasons of computational efficiency.For such calculations, we recommend switching to the random displacement model (Eq.17) rather than solving the LPDM Eq. ( 2).The reason for this recommendation is apparent in Figs.4-6, where the numerical error for all of the schemes investigated is seen to exceed the difference between RDM and LPDM solutions when the time step t 0.02 h/u * .Given that the unit of time in our non-dimensionalisation is T = h/u * , where h = 100 − 1000 m is boundary layer height and u * = 0.1−1 m s −1 is surface friction velocity, for a typical T ≈ 1000 s errors will be minimised by using the RDM whenever a time step t 20 s is required.
Naturally, the recommendations above are based only on the limited set of schemes which we have studied.It is to be hoped that the protocol and test cases introduced here will be helpful to other researchers developing and testing novel www.geosci-model-dev.net/9/2441/2016/Geosci.Model Dev., 9, 2441-2457, 2016 methods for LPDMs.A key challenge in such development will be the careful treatment of reflection boundary conditions, including their generalisation to more complex physical situations (e.g.Wilson and Flesch, 1993;Thomson et al., 1997;Wilson and Yee, 2007).

Code availability
The MATLAB source code of the FPE solver can be found online via GitHub and by searching for the repository "MRE FPE solver" (https://github.com/nhramli/MRE-FPEsolver.git).

Figure 5 .
Figure 5. L 2 error (Eq.20) as a function of non-dimensional time step t u * / h for the stable ABL test case integrated at intermediate time t = h/u * (a) and at late time t = 4 h/u * (b).From left to right, blue lines of slopes 1, 2, and 3 are plotted for reference.

Figure 6 .
Figure 6.L 2 error (Eq.20) as a function of non-dimensional time step t u * / h for the neutral ABL test case integrated at intermediate time t = 3 h/u * (a) and at late time t = 12 h/u * (b).From left to right, blue lines of slopes 1 and 2 are plotted for reference.

Table 1 .
The natural non-dimensionalisation has length, velocity, and timescales of ABL height h, surface friction velocity u * , and h/u * respectively.Under this non-Vertical profiles of vertical velocity fluctuations σw (z) (a) and vertical velocity Lagrangian decorrelation time τ (z) (b) used in the test-case problems (see Table1).The dimensions for σw and τ are frictional velocity u * and h/u * respectively, where h is the ABL height.The non-dimensional profiles of σ w (z) and τ (z) suitable for (i) a constant τ profile, (ii) a stable ABL, and (iii) a neutral ABL (e.g. Figure 7. Snapshots of reconstructed particle density c(z, t) for the stable ABL case at time t = h/u * , shown at each scheme.(a) When long time step t = 0.05 h/u * is used and errors due to boundary conditions dominate.(b) When moderate time step t = 0.007 h/u * is used.