IL-GLOBO (1.0) - integrated Lagrangian particle model and Eulerian general circulation model GLOBO: development of the vertical diffusion module

Abstract. The development and validation of the vertical diffusion module of IL-GLOBO, a Lagrangian transport model coupled online with the Eulerian general circulation model GLOBO, is described. The module simulates the effects of turbulence on particle motion by means of a Lagrangian stochastic model (LSM) consistently with the turbulent diffusion equation used in GLOBO. The implemented LSM integrates particle trajectories, using the native σ-hybrid coordinates of the Eulerian component, and fulfils the well-mixed condition (WMC) in the general case of a variable density profile. The module is validated through a series of 1-D offline numerical experiments by assessing its accuracy in maintaining an initially well-mixed distribution in the vertical. A dynamical time-step selection algorithm with constraints related to the shape of the diffusion coefficient profile is developed and discussed. Finally, the skills of a linear interpolation and a modified Akima spline interpolation method are compared, showing that both satisfy the WMC with significant differences in computational time. A preliminary run of the fully integrated 3-D model confirms the result only for the Akima interpolation scheme while the linear interpolation does not satisfy the WMC with a reasonable choice of the minimum integration time step.


Introduction
Global-(or hemispheric-) scale transport is recognised as an important issue in air pollution and climate change studies.Pollutants can travel across continents and have an influence even far from their source (see, for recent examples, Fiore et al., 2011;Yu et al., 2013).Moreover, transport of volcanic emissions (e.g. the recent Eyjafjallajökull eruption) or accidental hazardous releases (like the Fukushima and Chernobyl nuclear accidents) are also important at the global scale.
The natural framework for the description of tracer transport inflows is the Lagrangian approach (see, for example, the seminal works by Taylor, 1921, andRichardson, 1926).In the Lagrangian framework, the tracer transport is described by integrating the kinematic equation of motion for fluid "particles" in a given flow velocity field, provided by, e.g. a meteorological model.The turbulent motion unresolved by Eulerian equations for averaged quantities (in the Reynolds or volume-filtered sense) can be accounted for by including a stochastic component into the kinematic equation.
The stochastic component can be added to the particle position to give the Lagrangian equivalent of the Eulerian advection-diffusion equation.This kind of model is usually called a random displacement model (RDM) and is suitable for dispersion over long timescales.When the stochastic component is added to the velocity, the model is usually called a random flight model (RFM), which is more suitable for shorter time dispersion.In both cases, the stochastic model formulation has to be consistent with some basic physical requirements (Thomson, 1987(Thomson, , 1995)).
Various Lagrangian transport models exist which can be used at the global scale.Some are designed specifically for the description of atmospheric chemistry (Reithmeier and Sausen, 2002;Wohltmann and Rex, 2009;Pugh et al., 2012, see, e.g.), while others focus on the transport of tracers.In the latter class, two of the most widely used models are Published by Copernicus Publications on behalf of the European Geosciences Union.
FLEXPART (FLEXible PARTicle dispersion model) (Stohl et al., 2005) and HYSPLIT (Hybrid Single Particle Lagrangian Integrated Trajectory Model) (Draxler and Hess, 1998), which are highly flexible and can be easily used in a variety of situations.Both are compatible with different input types (usually provided by meteorological services like the European Centre for Medium-Range Weather Forecasts (ECMWF)), relying on their own parameterisation for fields not available from the meteorological model output.Models of this kind are suited for both forward and backward dispersion studies.
An alternative approach is to couple the Eulerian and Lagrangian parts online.On one hand, this makes the Eulerian fields available to the Lagrangian model at each Eulerian time step, increasing the accuracy for temporal scales shorter than the typical meteorological output interval.On the other hand, it also allows the consistent parameterisation of processes in the Eulerian and Lagrangian frameworks (e.g. the vertical dispersion in the boundary layer).Moreover, where the considered tracer may have an impact on meteorology (e.g. on radiation or cloud microphysics), online integration provides a natural way to include these effects (Baklanov et al., 2014).Online coupling also ensures the consistency of a mixed Eulerian-Lagrangian analysis of the evolution of atmospheric constituents (e.g.water or pollutants) along a trajectory (Sodemann et al., 2008;Real et al., 2010, see, e.g.).Malguzzi et al. (2011) recently developed a new global numerical weather prediction model, named GLOBO, based on a uniform latitude-longitude grid.The model is an extension to the global scale of the Bologna Limited Area Model (BO-LAM) (Buzzi et al., 2004), developed and employed during the early 90s.GLOBO is used for daily forecasting at the Institute of Atmospheric Sciences and Climate of the National Research Council of Italy (ISAC-CNR) and is also used to produce monthly forecasts.Online integration with BOLAM family models has already yielded interesting results in the development of the meteorology and composition model BOLCHEM (BOLam + CHEMistry) (Mircea et al., 2008).Considering that experience, the GLOBO model constitutes the natural basis for the further development of an integrated Lagrangian model.
In the following, the development of the vertical diffusion module is presented, focusing in particular on its compliance with basic theoretical requirements (Thomson, 1987(Thomson, , 1995, the well-mixed condition, see ) in connection with different numerical issues.In Sect. 2 the theoretical basis of the model formulation is given, while Sect. 3 describes different aspects of the numerical implementation.Finally, the model verification is presented and discussed in Sect. 4.

Lagrangian stochastic model formulation
In application to dispersion in turbulent flows, Lagrangian stochastic models (LSMs), Markovian at order M(M = 0, 1, . ..), are described by a set of stochastic differential equations (SDEs).The equation for the Mth order derivative M is dX where i and j indicate the components and X (k) i is the kthorder time-derivative of the Lagrangian Cartesian coordinate component X i ≡ X (0) i .Coefficients a i and b ij are called drift and Wiener coefficients, respectively.The remaining equations of the set (1 ≤ k ≤ M) are described by dX (2) The set of equations is equivalent to the Fokker-Planck equation: where A k i = 1 for k < M and A k i = a i for k = M, x i is the Eulerian equivalent of X i and K ij ≡ b ik b j k /2 (Thomson, 1987).Equation (3) describes the evolution of the probability density function p(x (0) , . .., x (M) , t), where 3 ).For the evolution of (X (0) , . .., X (M) ) to be approximated by a Markov process, the time correlation of the variable X (M+1) has to be much shorter than the characteristic evolution time of X (M) .If the model has to describe the evolution of dispersion at time t τ , where τ is the correlation time of turbulent velocity fluctuations, the process is well captured at order M = 0.When shorter times are considered, as in the case of dispersion from a single point source before the Taylor (1921) diffusive regime occurs (t ≤ τ ), order M must be increased to 1.The model of lowest order (M = 0) is referred to as random displacement model (RDM) and is sufficiently accurate to describe the transport and mixing of particles at a time and space resolution typical of a global model.
The correct formulation of a RDM in a variable density flow was first obtained by Venkatram (1993) and then refined and generalised by Thomson (1995) and is briefly recalled here.Equation ( 3) is valid for the probability density function p of particle position with the initial condition p((x), t)| t=t 0 = p((x), t 0 ).Since the ensemble average concentration c is proportional to p, Eq. (3) can be rewritten as If c ∝ ρ at some time t , where ρ is the ensemble average of air density, then for all t > t the two quantities must remain proportional.This condition, called well-mixed condition (WMC) after Thomson (1987), implies that ρ is also a solution of Eq. ( 4).Substituting c with ρ in Eq. ( 4) and using the continuity equation where u i is the density weighted mean velocity, defined as (Thomson, 1995): the following expression is obtained: Then, integrating both sides and rearranging gives where the non-uniqueness implied by the integration is removed considering that in the well-mixed state, the mixing ratio flux must be proportional to u i ρ .Substituting Eq. ( 8) into Eq.( 4) gives the equivalent of Eq. ( 2) in Thomson (1995).
At the coarse resolution typical of global models, vertical motions can be considered decoupled from the horizontal ones.Therefore, only the vertical coordinate x 3 ≡ z (and X 3 ≡ Z in Lagrangian terms) need to be considered.In this case, the RDM reduces to a single differential stochastic equation where w ≡ u 3 and K ≡ K 33 .

Numerical implementation of the vertical diffusion module
In its final form, IL-GLOBO is designed to be a fully online integrated model (or at least an online-access model, according to Baklanov et al., 2014), where the different components share the same "view" of the atmosphere, i.e. use the same discretisation, parameterisations, etc.The development of the vertical diffusion module is based on this principle.

Vertical coordinate
Within IL-GLOBO, the Lagrangian equations are integrated in the same coordinate system used in the Eulerian model.This choice maintains the consistency between the Lagrangian and Eulerian components and reduces the interpolation errors and computational cost.GLOBO uses a hybrid vertical coordinate system in which the terrain-following coordinate σ (0 < σ < 1) smoothly tends, with height above the ground, to a pressure coordinate P , according to where P 0 is a reference pressure (typically 1000 hPa), P S is the surface pressure and α is a parameter that gives the classical σ coordinate for α = 1 (Phillips, 1957).The parameter α depends on the model orography and, therefore, on resolution.It is limited by the condition ∂ ∂p σ ≥ 0 that results in the relationship: which is satisfied by the typical setting α = 2, used for a wide range of resolutions in GLOBO applications (Malguzzi et al., 2011).
The vertical Lagrangian coordinate is identified by , corresponding to the vertical coordinate σ , and is connected to the Lagrangian vertical position Z above the ground through Eq. ( 10) and the hydrostatic relationship.In the meteorological component, the height above the ground z is a diagnostic quantity that can be derived from the geopotential through z(σ ) = ( (σ )− g )g −1 , where g is the geopotential at the height of roughness length.Since the determination of the different terms in Eq. ( 9) involves discrete Eulerian fields and their numerical derivatives, the choice of employing σ also has the advantage of making interpolation straightforward and consistent with the Eulerian part.
Because σ (z) is not linear (σ is not a Cartesian coordinate system), the stochastic chain rule (see, e.g.Kloeden and Platen, 1992, p. 80) must be used to derive the correct form of Eq. ( 9) for , giving where ω is the vertical velocity in the σ coordinate system and z is the Cartesian vertical coordinate.The last term in square brackets stems from the Itô-Taylor expansion of order dW 2 , which must be included for the correct description at order dt (Gardiner, 1990, p. 63).

Discretisation and interpolation
The GLOBO prognostic variables are computed on a Lorenz (1960) vertical grid: all the quantities are on "integer" levels σ i , except vertical velocity, turbulent kinetic energy and mixing length and, consequently, diffusion coefficients, located at "semi-integer" levels σ h i (see Fig. 1).In typical applications, the GLOBO vertical grid is regularly spaced in σ (Malguzzi et al., 2011), although it is possible to use a variable grid spacing, as in its limited-area version BOLAM (Buzzi et al., 1994).for the first order derivative.Following the same considerations made for ρ, the derivatives of σ with respect of z are computed from relationships similar to Eq. ( 13) and Eq. ( 14).

255
For the highly varying K profiles, two different methods are tested, the first with two variants.The first method interpolates the function linearly at the particle position, and uses finite differences derivatives.In the first variant (labeled D), the first order two-points derivative is computed and kept 260 constant between two grid points.In order to give a smoother description of the derivatives, a variant (labeled D ′ ) is also tested in which the three-points centered derivative is computed and interpolated linearly at the particle position.For D ′ , the values of first order derivative at the lowest boundary 265 is computed as: This is assumed because K is expected to be linear near the surface, according to Monin-Obukhov similarity theory where: for the neutral case, with proper modifications for diabatic cases.
The second method (labeled A) is based on the Akima (1991) cubic spline.For each interval it considers the pre-275 vious and the next two adjacent intervals (for a total number of 6 grid points) to compute the coefficients of the interpolating cubic polynomial.This algorithm reduces the number of oscillations in the interpolating function compared to regular cubic splines and enforces the linearity when 4 points 280 are collinear (Akima, 1991).Using this property, a linear profile near the ground is imposed to the interpolating function by adding two fictitious points below the ground that are collinear with the two lower grid points of the domain.In addition, to ensure the positivity of the interpolating functions, 285 the local algorithm of Fischer et al. (1991) is used, which also preserves the continuity of first order derivatives.

Integration scheme and time-step selection
The most common integration scheme for SDE in atmospheric transport models is the Euler-Maruyama forward 290 scheme: The coefficients a and b come from Equation ( 12).The Euler-Maruyama forward scheme is the simplest strong Taylor approximation and turns out to be of order of strong con-295 vergence γ = 0.5 (Kloeden and Platen, 1992, p. 305).By a rather simple modification of the Euler-Maruyama scheme, i.e. adding the term: where b ′ is the first-order derivative of b, the Milstein scheme 300 is obtained, which is of order of strong convergence γ = 1.
It is worth noting that the strong order γ = 1 of the Milstein scheme corresponds to the strong order γ = 1 of the Euler deterministic scheme.Therefore, Milstein can be regarded as the correct generalization of the deterministic Euler scheme 305 (Kloeden and Platen, 1992, p. 345).The additional term uses only already computed quantities involved in the determination of the drift term of Equation ( 12).Preliminary idealized tests do not show any appreciable accuracy improvement with respect to the Euler-Maruyama scheme.However, because they confirm the negligible extra computational cost of this method, the Milnstein scheme will be used to integrate the model.
In the meteorology component of IL-GLOBO, the Eulerian equations are solved with a macro time-step ∆T , which 315 depends basically on the horizontal resolution due to the limitations imposed by the Courant number.Other timesteps are involved in the Eulerian part but are not relevant here.In typical implementations, ∆T ranges from 432 s for 362 × 242 point resolution (used for monthly forecasts 1 ) to 320 150 s for 1202 × 818 point resolution (used for high resolution weather forecasts 2 ).The macro time-step is taken as the upper limit for the solution of Equation ( 12).The time-step needed to reach the required accuracy depends on the quantities involved in determining the various elements in Equa-325 tion (17).
1 http://www.isac.cnr.it/dinamica/projects/forecastdpc/month 2 http://www.isac.cnr.it/dinamica/projects/forecasts/glob Figure 1.Schematic representation of field value distributions between integer (continuous lines) and semi-integer (dashed lines) levels in the GLOBO model.being a continuous coordinate, the quantities needed to compute the terms of Eq. ( 12) must be interpolated from the Eulerian fields given at discrete levels.The computation of first-and second-order derivatives of Eulerian model quantities is also required in the implementation of the LSM.Interpolation and derivation algorithms can influence both the accuracy and the computational cost of the Lagrangian model and thus require careful assessment.
For density ρ and geopotential , linear interpolation and central differences derivative are used assuming that those fields are regular enough.At the lower boundary, it is required that which implies for the first-order derivative.Following the same considerations made for ρ, the derivatives of σ with respect of z are computed from relationships similar to Eqs. ( 13) and ( 14).
For the highly varying K profiles, two different methods are tested, the first with two variants.The first method interpolates the function linearly at the particle position and uses finite differences derivatives.In the first variant (labelled D), the first-order two-point derivative is computed and kept constant between two grid points.In order to give a smoother description of the derivatives, a variant (labelled D ) is also tested in which the three-point centered derivative is computed and interpolated linearly at the particle position.
For D , the values of the first-order derivative at the lowest boundary are computed as This is assumed because K is expected to be linear near the surface, according to Monin-Obukhov similarity theory where for the neutral case, with proper modifications for diabatic cases.
The second method (labelled A) is based on the Akima (1991) cubic spline.For each interval it considers the previous and the next two adjacent intervals (for a total number of six grid points) to compute the coefficients of the interpolating cubic polynomial.This algorithm reduces the number of oscillations in the interpolating function compared to regular cubic splines and enforces the linearity when four points are collinear (Akima, 1991).Using this property, a linear profile near the ground is imposed to the interpolating function by adding two fictitious points below the ground that are collinear with the two lower grid points of the domain.In addition, to ensure the positivity of the interpolating functions, the local algorithm of Fischer et al. (1991) is used, which also preserves the continuity of first-order derivatives.

Integration scheme and time-step selection
The most common integration scheme for SDE in atmospheric transport models is the Euler-Maruyama forward scheme: The coefficients a and b come from Eq. ( 12).The Euler-Maruyama forward scheme is the simplest strong Taylor approximation and turns out to be of the order of strong convergence γ = 0.5 (Kloeden and Platen, 1992, p. 305).By a rather simple modification of the Euler-Maruyama scheme, i.e. adding the term: where b is the first-order derivative of b, the Milstein scheme is obtained, which is of the order of strong convergence γ = 1.It is worth noting that the strong order γ = 1 of the Milstein scheme corresponds to the strong order γ = 1 of the Euler deterministic scheme.Therefore, Milstein can be regarded as the correct generalisation of the deterministic Euler scheme (Kloeden and Platen, 1992, p. 345).The additional term uses only already-computed quantities involved in the determination of the drift term of Eq. ( 12).Preliminary idealised tests do not show any appreciable accuracy improvement with respect to the Euler-Maruyama scheme.However, because they confirm the negligible extra computational cost of this method, the Milstein scheme will be used to integrate the model.

Geosci
In the meteorology component of IL-GLOBO, the Eulerian equations are solved with a macro time step T , which depends basically on the horizontal resolution due to the limitations imposed by the Courant number.Other time steps are involved in the Eulerian part but are not relevant here.In typical implementations, T ranges from 432 s for 362 × 242 point resolution (used for monthly forecasts1 ) to 150 s for 1202×818 point resolution (used for high-resolution weather forecasts2 ).The macro time step is taken as the upper limit for the solution of Eq. ( 12).The time step needed to reach the required accuracy depends on the quantities involved in determining the various elements in Eq. ( 17).
First, a straightforward constraint is that the time step must satisfy the relationship (Wilson and Yee, 2007, see, e.g.), which expresses the requirement that the average root mean square step length must be much smaller than the scale of the variations of K.This gives rise to a limitation that is consistent with the surfacelayer behaviour of the diffusion coefficient, Eq. ( 16).The condition expressed by Eq. ( 19) makes t 1 vanish for z → 0. Such behaviour ensures that the WMC is satisfied theoretically, but clearly poses problems for numerical implementation (Ermak and Nasstrom, 2000;Wilson and Yee, 2007).However, in the application of a global model, where particles can be distributed throughout the troposphere, this problem affects only a small fraction of particles in the vicinity of the surface.Therefore, it can be dealt with by selecting a t min small enough for the solution to be within the accepted error and, at the same time, large enough to not impact the overall computational cost.
In addition to Eq. ( 19), another constraint is needed to account also for the presence of maxima in the K profile, which must be present if one considers the whole atmosphere.At maxima (or minima), Eq. ( 19) gives an unlimited t 1 , which is not suitable for the integration of the model as it could cause the trajectory to cross the maximum (or minimum), with a significant change in K(z) associated to a change in ∂ z K sign.To avoid this problem, a further constraint is introduced, based on the normalised second-order derivative, which gives an estimation of the width of the maximum.The constraint reads The above equation has the property of limiting t 2 according to the sharpness of the K peak.
Taking the minimum among T , t 1 and t 2 (and replacing with = C T in Eqs.19 and 20), gives where the parameter C T quantifies the "much less" condition and, therefore, must be at least 0.1 or smaller.
Figure 2 shows the application of Eq. ( 21) for a K profile representative of GLOBO (see Sect. 4) and a C T = 0.01.The t decreases in the presence of K gradients thanks to condition ( 19), and is limited around the K maximum (where ∂K/∂σ = 0) by condition (20).The maximum of t = T is attained at higher levels.
It should be kept in mind that the method is based on local quantities and may fail if strong variations of K occur in one time step along the particle path.To overcome this problem, an additional constraint is used to make the algorithm nonlocal (or "less local").Using the t 0 computed at the particle position at time t, two other time steps ( t + and t − ) are evaluated at the positions: The minimum t among t 0 , t + and t − is then used to advance the particle position t+ t .

Boundary conditions
The necessary boundary condition for the conservation of the probability (and therefore of the mass) is the reflective boundary (Gardiner, 1990, p. 121).Wilson and Flesch (1993) show that the elastic reflection ensures the WMC if the integration time step is small enough.However, in cases of nonhomogeneous K, numerical implementation requires that t vanishes as the particle approaches the boundary.For models that focus on near-surface dispersion, the time step needed to achieve the required accuracy can become very small.Ermak and Nasstrom (2000) describe a theoretically well-founded method to speed up (roughly by a factor of 10) simulations of this kind.
In the case of IL-GLOBO, it will be shown that the elastic reflection condition at σ = 1, coupled with the adaptive time-step algorithm described in Sect.3.3, can ensure a good approximation of the solution while maintaining affordable the computational cost.

Model verification: the well-mixed condition
In order to verify the vertical diffusion module of IL-GLOBO, a series of experiments was performed with a 1-D version of the code and then tested in a preliminary version of the full 3-D model.19), the blue line the contribution of Eq. ( 20), and the black line the combined condition (Eq.21, with ∆T = 432 s and CT = 0.01).
the problem, an additional constraint is used to make the algorithm non-local (or less local).Using the ∆t 0 computed at 375 the particle position at time t, two other time-step (∆t + and ∆t − ) are evaluated at the positions: The minimum ∆t among ∆t 0 , ∆t + and ∆t − is then used to advance the particle position Σ t+∆t .380

Boundary conditions
The necessary boundary condition for the conservation of the probability (and therefore of the mass) is the reflective boundary (Gardiner, 1990, p. 121).Wilson and Flesch (1993) show that the elastic reflection ensures the WMC if 385 the integration time-step is small enough.However, in cases of non-homogeneous K, numerical implementation requires that ∆t vanishes as the particle approaches the boundary.For models that focus on near surface dispersion, the time-step needed to achieve the required accuracy can become very 390 small.Ermak and Nasstrom (2000) describe a theoretically well founded method to speed-up (roughly by a factor of 10) simulations of this kind.
In the case of IL-GLOBO, it will be shown that the elastic reflection condition at σ = 1, coupled with the adaptive time 395 step algorithm described in Section 3.3, can ensure a good approximation of the solution while maintaining affordable the computational cost.

Model verification: the well-mixed condition
In order to verify the vertical diffusion module of IL-400 GLOBO, a series of experiments was performed with a 1-D version of the code and then tested in a preliminary version  19), the blue line the contribution of Eq. ( 20) and the black line the combined condition (Eq.21, with T = 432 s and C T = 0.01).
of 362 × 242 cells and 50 vertical levels evenly spaced in σ ) starting at 11 March 2011 00:00 UTC.After 36 h of simulation (12:00 UTC), averaging of σ = const surfaces were performed for K, ρ and , obtaining vertical profiles as a function of σ .Fields of ρ and were averaged over the whole domain.As far as K is concerned, averages were performed for latitudes between +60 and −60 • North in daytime (longitudes between −45 and +45 • East) and night-time (longitudes between +135 and −135 • East) conditions, over land and sea separately.The most intense K profile is selected which corresponds to the daytime conditions over land.Profiles of ρ and z are rather smooth and regular over space and time, while K displays large variability.The profiles were fitted with analytical functions derived by combining the hydrostatic equation and the perfect gas law.The following analytical expressions were used: and with T 0 = 288.0K, ρ 0 = 1.2 kg m −3 and = −0.007K m −1 .As a consequence of the hydrostatic perfect gas assumption, by expressing the density ρ in sigma vertical units ρ σ = ρ dz dσ and using Eqs.( 24) and ( 23), the following constant value is obtained: Figure 3 shows the GLOBO-average profiles and their fitting functions for the density ρ and the geopotential height g −1 as function of σ .As far as the K profile is concerned, the function ρ and z are rather smooth and regular over space and time, while K displays a large variability.The profiles were fitted with analytical functions derived combining the hydrostatic equation and the perfect gas law.The following analytical expressions were used: and: with T 0 = 288.0K, ρ 0 = 1.2 kgm − 3 and Γ = −0.007K m −1 .As a consequence of the hydrostatic perfect gas assump-425 tion, by expressing the density ρ in sigma vertical units (ρ σ = ρ dz dσ ) and using Equations ( 24) and ( 23), the following constant value is obtained: Figure 3 shows the GLOBO averaged profiles and their fit- As far as th  used to account for the specific K features: it should display a linear behaviour near the surface, must tend to zero near the boundary layer top 3 and, therefore, must display a maximum at some height.In Eq. ( 26), A = 0.29 ms −1 was first determined according to average surface-layer properties (the first GLOBO vertical level), and corresponds to a friction velocity u * 0.7 ms −1 .Then, the other two parameters were allowed to vary to fit the average profile giving B = 1.3 × 10 −3 m −1 and C = 1.6.
Although the above profile is representative of the typical GLOBO diffusivity, real profiles can be remarkably less regular, creating challenging conditions for the model.For this reason, a profile was selected among those showing isolated strong maxima near the ground.This is typical of strong convective conditions just after sunrise.Fitting Eq. ( 26) to this second profile gives A = 0.3 ms −1 , B = 4.0 × 10 −3 m −1 and C = 4.5.Figure 4 reports the GLOBO "average" and "peaked" K profiles as function of σ .

Determination of the optimal setting for the adaptive time-step selection algorithm
The first series of experiments concerns the optimisation of the adaptive scheme for t, i.e. the selection of the best suited value for the coefficient C T in Eq. ( 21).The 'average' profile is shown in red, while the 'peaked' profile is shown in green.The functional form of both profiles is described by Eq. ( 26).
As far as the K profile is concerned, the function is used to account for the specific K features: it should dis-435 play a linear behavior near the surface, must tend to zero near the boundary layer top 3 and, therefore, must display a maximum at some height.In Equation ( 26), A = 0.29 ms −1 was first determined according to average surface-layer properties (the first GLOBO vertical level), and corresponds to 440 a friction velocity u * ≃ 0.7 ms −1 .Then, the other two parameters were let to vary to fit the average profile giving B = 1.3 × 10 −3 m −1 and C = 1.6.
Although the above profile is representative of the typical GLOBO diffusivity, real profiles can be remarkably less regular, creating challenging conditions for the model.For this reason, a profile was selected among those showing isolated strong maximum near the ground.This is typical of strong convective conditions just after sunrise.Fitting Equation (26) on this second profile gives A = 0.3 ms −1 , B = 450 4.0 × 10 −3 m −1 and C = 4.5.Figure 4 reports the GLOBO 'averaged' and 'peaked' K profiles as function of σ.

Determination of the optimal setting for the adaptive time-step selection algorithm
The first series of experiments concerns the optimization of 455 the adaptive scheme for ∆t, i.e., the selection of the best suited value for the coefficient C T in Equation 21.Simulations were performed in flow conditions described by Equations ( 23), ( 24) and ( 26), distributing particles with number concentration proportional to ρ.For the WMC to be 460 3 In GLOBO, K also accounts for a part of the instability generated by moist convection and therefore it may not vanish at the boundary layer top.The "average" profile is shown in red, while the "peaked" profile is shown in green.The functional form of both profiles is described by Eq. ( 26).
Simulations were performed in flow conditions described by Eqs. ( 23), ( 24) and ( 26), distributing particles with number concentration proportional to ρ.For the WMC to be satisfied, this distribution must remain constant as the time evolves.Equation ( 12) was integrated for 4 × 10 5 particles and for 200 macro time steps, each 432 s long, for a total of T = 86 400 s = 24 h.The actual time step used is given by Eq. ( 21) with the additional lower limit t min = 0.01.Simulations were performed using 12 cores of an Intel Xeon machine.Since the initial condition was already well-mixed (C ∝ ρ), the simulation time was considered sufficient to assess the skill of the model in satisfying the WMC.At the end of the simulation, final concentration profiles were computed in "σ volume", i.e. c(σ ) = N(σ )( σ ) −1 , where N(σ ) is the number of particles between σ and σ + σ .The skill of the model in reproducing the WMC was evaluated using the root mean square error (RMSE) of the final normalised concentration profile with respect to the normalised density profile (derived using Eq.25).
Figure 5 reports the different profiles of concentration after 24 h of simulation computed using different values of C T .The shaded region represents the interval between 3 standard deviations from the expected value.RMSE values for each simulation are reported in Table 1 along with the computation time.The RMSE error becomes comparable to the statistical error for C T = 0.01, which is selected as the optimal value.In order to evaluate the possible dependency of C T on the number of particles, two additional sets of runs were performed with 10 5 and 16 × 10 5 particles that correspond to halving and doubling, respectively, the statistical error of the base experiment.Results are reported in Fig. 6 which shows that, in the considered range, the optimal C T is quite independent of the number of particles.
It is worth noting that the time-step selection algorithm with the proper choice of C T ensures that the WMC is also satisfied at the reflective boundary too, as mentioned in Sect.3.4.satisfied, this distribution must remain constant as the time evolves.Equation ( 12) was integrated for 4 × 10 5 particles and for 200 macro time-steps, each 432 s long, for a total of T = 86400 s = 24 h.The actual time-step used is given by Equation ( 21) with the additional lower limit ∆t min = 0.01.

465
Simulations were performed using 12 cores of an Intel Xeon machine.Since the initial condition was already well-mixed (C ∝ ρ), the simulation time was considered sufficient to assess the skill of the model in satisfying the WMC.At the end of the simulation, final concentration profiles were computed 470 in "σ volume", i.e., c(σ) = N (σ)(∆σ) −1 , where N (σ) is the number of particles between σ and σ + ∆σ.The skill of the model in reproducing the WMC was evaluated using the root mean square error (RMSE) of the final normalized concentration profile with respect to the normalized density profile 475 (derived using Equation 25).
Figure 5 reports the different profiles of concentration after 24 hours of simulation computed using different values of C T .The shaded region represents the interval between 3 standard deviations from the expected value.RMSE values 480 for each simulation are reported in Table 1 along with the computation time.The RMSE error becomes comparable to the statistical error for C T = 0.01, which is selected as the optimal value.In order to evaluate the possible dependency of C T on the number of particles, two additional sets of runs 485 were performed with 10 5 and 16 × 10 5 particles that correspond to halving and doubling, respectively, the statistical error of the base experiment.Results are reported in Figure 6 which shows that, in the considered range, the optimal C T is quite independent of the number of particles.

490
It is worth noting that the time-step selection algorithm with the proper choice of C T ensures that the WMC is also satisfied at the reflective boundary too, as mentioned in Section 3.4.satisfied, this distribution must remain constant as the time evolves.Equation ( 12) was integrated for 4 × 10 5 particles and for 200 macro time-steps, each 432 s long, for a total of T = 86400 s = 24 h.The actual time-step used is given by Equation ( 21) with the additional lower limit ∆t min = 0.01.

465
Simulations were performed using 12 cores of an Intel Xeon machine.Since the initial condition was already well-mixed (C ∝ ρ), the simulation time was considered sufficient to assess the skill of the model in satisfying the WMC.At the end of the simulation, final concentration profiles were computed 470 in "σ volume", i.e., c(σ) = N (σ)(∆σ) −1 , where N (σ) is the number of particles between σ and σ + ∆σ.The skill of the model in reproducing the WMC was evaluated using the root mean square error (RMSE) of the final normalized concentration profile with respect to the normalized density profile 475 (derived using Equation 25).
Figure 5 reports the different profiles of concentration after 24 hours of simulation computed using different values of C T .The shaded region represents the interval between 3 standard deviations from the expected value.RMSE values 480 for each simulation are reported in Table 1 along with the computation time.The RMSE error becomes comparable to the statistical error for C T = 0.01, which is selected as the optimal value.In order to evaluate the possible dependency of C T on the number of particles, two additional sets of runs 485 were performed with 10 5 and 16 × 10 5 particles that correspond to halving and doubling, respectively, the statistical error of the base experiment.Results are reported in Figure 6 which shows that, in the considered range, the optimal C T is quite independent of the number of particles.

490
It is worth noting that the time-step selection algorithm with the proper choice of C T ensures that the WMC is also satisfied at the reflective boundary too, as mentioned in Section 3.4.

Evaluation of the interpolation algorithms 495
In the subsequent set of experiments, the model skill in reproducing the WMC was evaluated for the interpolation techniques D, D ′ and A described in Section 3.2.
In the first experiment, the analytical fields described by Equations ( 23), ( 24) and ( 26) with the parameters of the 'av-500 erage' diffusivity profile were resampled on a 50 point regular grid.This provides a discrete version of the experiment described in the previous section, with the same vertical resolution of the GLOBO original fields.
The particle number, initial distribution and simulation 505 time are the same as in the experiment described in section 4.1.The integration time-step is selected using the local algorithm.The time-step selection algorithm requires the computation of the second order derivative of K, which is not possible for the D interpolation scheme.Therefore, it is estimated using finite differences of the first order derivative.The results of this experiment are shown in Figure 7.In the upper panel, the integration time-step profiles of the three simulations and the Akima interpolated diffusion coefficient profile, are displayed.The lower panel shows the normalized 515 distribution of the particle after 24 hours of simulation along with the expected value.Table 2 displays the integration time and RMSE obtained for the various experimental settings.
The time-step profiles are similar, except for the A profile around the region of maximum of K, where it shows strong 520 variations and, on the average, is longer than the others.Looking at the distribution of particles (lower panel), it re-

Evaluation of the interpolation algorithms
In the subsequent set of experiments, the model skill in reproducing the WMC was evaluated for the interpolation techniques D, D and A described in Sect.3.2.
In the first experiment, the analytical fields described by Eqs. ( 23), ( 24) and ( 26) with the parameters of the "average" diffusivity profile were resampled on a 50-point regular grid.This provides a discrete version of the experiment described in the previous section, with the same vertical resolution of the GLOBO original fields.
The particle number, initial distribution and simulation time are the same as in the experiment described in Sect.4.1.The integration time step is selected using the local algorithm.The time-step selection algorithm requires the computation of the second-order derivative of K, which is not possible for the D interpolation scheme.Therefore, it is estimated using finite differences of the first-order derivative.The results of this experiment are shown in Fig. 7  sults that simulations with A and D interpolation algorithms both satisfy the WMC within the statistical limit, while the simulation with the D ′ algorithm fails to maintain the well 525 mixed state, in particular near the ground.Additional experiments (not reported) show that in order to obtain a well mixed solution with D ′ , resolution must be doubled, at least.The problem is probably related to the definition of derivatives of K between grid points.In fact, although D ′ computes 530 derivatives at higher order of approximation than D, they are not consistent with a linear variation of K.Although the use of D ′ can be appropriate for slowly varying and monotone functions like ρ and z, it turns out to be unsuitable for the more complex K profile which, in addition, affects both the 535 Wiener stochastic term and the drift term.For these reasons, the D ′ interpolation scheme is not used in the following experiments.
The second experiment concerns the 'peaked' profile.In this case, the K profile is used directly, without the resam-540 pling of the fitting function.Simulations with A and D algorithms were performed with both local and non-local timestep selection algorithm.Figure 8 reports the time-step and concentration profiles, while execution times and RMSEs are shown in Table 3.Although the integration time-step profiles 545 look very similar for the local and non-local algorithms, the small differences have large impact on the results: the local algorithm strongly fails in reproducing the WMC for both  interpolation schemes, especially for D. Conversely, the nonlocal algorithm turns to be effective in selecting the ap-550 propriate time-step even in presence of strong gradients and isolated maxima.This is reflected on its higher computational cost (see Table 3).

Implementation on the 3-D model
A preliminary test of the algorithms on the 3-D model has 555 been performed.The interpolation algorithm has been implemented in a simplified quasi-1-D form, where the diffusion coefficient has been considered to be horizontally constant between grid points.IL-GLOBO uses the same parallelization of GLOBO, with particle exchanged between processes 560 at each macro time-step.Particles are first advected horizontally for a macro time-step using their deterministic velocity, and then 'diffused' in the vertical according to Equation (12).
After 12 h of spinup, 5 × 10 5 particles are released with a vertical distribution proportional to the average density pro-565 file, and randomly and homogeneously distributed in the horizontal.Particle statistics are computed after 24 h from the release.
A and D interpolation algorithm were tested using the nonlocal time-step selection.It is found that, while interpolation 570 scheme A maintains the WMC reasonably (RMSE=0.024), the time-step selection algorithm for scheme D requires extremely short time-steps ( ≪ ∆t min , see Section 4.1) in the distribution of the particle after 24 h of simulation along with the expected value.Table 2 displays the integration time and RMSE obtained for the various experimental settings.
The time-step profiles are similar, except for the A profile around the region of the maximum of K, where it shows strong variations and, on the average, is longer than the others.Looking at the distribution of particles (lower panel), it can be observed that simulations with A and D interpolation algorithms both satisfy the WMC within the statistical limit, while the simulation with the D algorithm fails to maintain the well-mixed state, in particular near the ground.Additional experiments (not reported) show that in order to obtain a well-mixed solution with D , resolution must be doubled at least.The problem is probably related to the definition of derivatives of K between grid points.In fact, although D computes derivatives at a higher order of approximation than D, they are not consistent with a linear variation of K.Although the use of D can be appropriate for slowly varying and monotone functions like ρ and z, it turns out to be unsuitable for the more complex K profile which, in addition, affects both the Wiener stochastic term and the drift term.For these reasons, the D interpolation scheme is not used in the following experiments.
The second experiment concerns the "peaked" profile.In this case, the K profile is used directly, without the resampling of the fitting function.Simulations with A and D algorithms were performed with both local and non-local sults that simulations with A and D interpolation algorithms both satisfy the WMC within the statistical limit, while the simulation with the D ′ algorithm fails to maintain the well 525 mixed state, in particular near the ground.Additional experiments (not reported) show that in order to obtain a well mixed solution with D ′ , resolution must be doubled, at least.The problem is probably related to the definition of derivatives of K between grid points.In fact, although D ′ computes 530 derivatives at higher order of approximation than D, they are not consistent with a linear variation of K.Although the use of D ′ can be appropriate for slowly varying and monotone functions like ρ and z, it turns out to be unsuitable for the more complex K profile which, in addition, affects both the 535 Wiener stochastic term and the drift term.For these reasons, the D ′ interpolation scheme is not used in the following experiments.
The second experiment concerns the 'peaked' profile.In this case, the K profile is used directly, without the resam-540 pling of the fitting function.Simulations with A and D algorithms were performed with both local and non-local timestep selection algorithm.Figure 8 reports the time-step and concentration profiles, while execution times and RMSEs are shown in Table 3.Although the integration time-step profiles 545 look very similar for the local and non-local algorithms, the small differences have large impact on the results: the local algorithm strongly fails in reproducing the WMC for both  propriate time-step even in presence of strong gradients and isolated maxima.This is reflected on its higher computational cost (see Table 3).

Implementation on the 3-D model
A preliminary test of the algorithms on the 3-D model has 555 been performed.The interpolation algorithm has been implemented in a simplified quasi-1-D form, where the diffusion coefficient has been considered to be horizontally constant between grid points.IL-GLOBO uses the same parallelization of GLOBO, with particle exchanged between processes 560 at each macro time-step.Particles are first advected horizontally for a macro time-step using their deterministic velocity, and then 'diffused' in the vertical according to Equation (12).
After 12 h of spinup, 5 × 10 5 particles are released with a vertical distribution proportional to the average density pro-565 file, and randomly and homogeneously distributed in the horizontal.Particle statistics are computed after 24 h from the release.
A and D interpolation algorithm were tested using the nonlocal time-step selection.It is found that, while interpolation 570 scheme A maintains the WMC reasonably (RMSE=0.024), the time-step selection algorithm for scheme D requires extremely short time-steps ( ≪ ∆t min , see Section 4.1) in the  3.Although the integration timestep profiles look very similar for the local and non-local algorithms, the small differences have a large impact on the results: the local algorithm mostly fails in reproducing the WMC for both interpolation schemes, especially for D. Conversely, the non-local algorithm turns out to be effective in selecting the appropriate time step, even in the presence of strong gradients and isolated maxima.This is reflected in its higher computational cost (see Table 3).

Implementation on the 3-D model
A preliminary test of the algorithms on the 3-D model has been performed.The interpolation algorithm has been implemented in a simplified quasi-1-D form, where the diffusion coefficient has been considered to be horizontally constant between grid points.IL-GLOBO uses the same parallelisation of GLOBO, with particles exchanged between processes at each macro time step.Particles are first advected horizontally for a macro time step using their deterministic velocity and then "diffused" in the vertical according to Eq. ( 12).
After 12 h of spin-up, 5 × 10 5 particles are released with a vertical distribution proportional to the average density profile, and randomly and homogeneously distributed in the horizontal.Particle statistics are computed after 24 h from the release.The initial distribution ∝ ρ (black) and the final distributions obtained using the A interpolation scheme (red) and the D interpolation scheme (blue).Dashed lines show the limit of 3 standard deviation around the initial distribution.

Geosci
region between σ = 0.9 and the lowest boundary (see Figure 9).Figure 10 shows the result of an experiment where 575 the WMC compliance of schemes A and D was tested with the lower limit for ∆t changed to ∆t min = 10 −5 for D. It can be observed that, for the D scheme, strong fluctuations are still present in the same region where the required timestep exceeds the lower limit.This is likely to be caused by the 580 occurrence of strong gradients that can be even larger than in the 'peaked' case, near points with extremely small values of K.In these cases, the A scheme interpolates with a smoother function which reduces the problem.

585
The development of a vertical Lagrangian diffusion model is presented.This constitutes the first step in building IL-GLOBO, a Lagrangian particle model integrated in the Eulerian global circulation model GLOBO.Critical details of the implementation have been analyzed and discussed.

590
The model is developed including the variable density term and the proper coordinate transformation term.The numerical scheme selected to integrate the SDE is the Milstein scheme, which is of order of strong convergence γ = 1.Therefore, it should be regarded as the natural extension of 595 the deterministic Euler scheme, in contrast to the so-called Euler-Maruyama scheme, which is merely the "transcription" of the deterministic Euler scheme, but not its equivalent.
An adaptive time-step scheme is proposed to ensure the 600 consistency of the model implementation with the WMC requirements.The time-step selection algorithm is limited not only by the condition imposed by the spatial scale of gradients, but also takes into account the scale of the width of maxima and minima of the diffusion coefficient, where the 605 former criterium fails.It is shown that this algorithm ensures that the error is within an acceptable range also at the reflecting boundaries.However, in case of isolated maxima, this scheme may fail.The implementation of a non-local algorithm, which evaluates ∆t in 2 additional points, is proposed 610 in order to solve the problem.Two numerical interpolation and derivation schemes are implemented and tested.The first is based on the linear interpolation of K and it is presented in two versions: one (D) keeps a constant first order derivative between two grid 615 points, while the other (D ′ ) uses linearly interpolated derivatives in the same interval.The second scheme (A) is based on a modified Akima (1991) interpolation algorithm with a local algorithm that ensures the positivity of the interpolating function (Fischer et al., 1991).

620
It is found that, although the method D ′ derivatives of higher order of approximation, it creates a local inconsistency between the linearly interpolated function and its derivatives and prevents the model from fulfilling the WMC.The other two schemes (D and A) both satisfy the WMC but 625 extremely peaked profiles of K may require the use of the non-local time-step selection algorithm.
A test with a preliminary implementation of the fully 3D model (IL-GLOBO) shows that, while the A scheme display a correct behavior, the D interpolation scheme requires and 630 extremely strong reduction of the integration time-step that prevents the WMC to be satisfied in reasonable time.A and D interpolation algorithms were tested using the non-local time-step selection.It is found that while interpolation scheme A maintains the WMC reasonably (RMSE = 0.024), the time-step selection algorithm for scheme D requires extremely short time steps ( t min , see Sect.4.1) in the region between σ = 0.9 and the lowest boundary (see Fig. 9). Figure 10 shows the result of an experiment where the WMC compliance of schemes A and D was tested with the lower limit for t changed to t min = 10 −5 for D. It can be observed that, for the D scheme, strong fluctuations are still present in the same region where the required time step exceeds the lower limit.This is likely to be caused by the occurrence of strong gradients that can be even larger than in the "peaked" case, near points with extremely small values of K.In these cases, the A scheme interpolates with a smoother function which reduces the problem.

Conclusions
The development of a vertical Lagrangian diffusion model is presented.This constitutes the first step in building IL-GLOBO, a Lagrangian particle model integrated in the Eulerian global circulation model GLOBO.Critical details of the implementation have been analysed and discussed.
The model is developed with the variable density term and the proper coordinate transformation term.The numerical scheme selected to integrate the SDE is the Milstein scheme, which is of the order of strong convergence γ = 1.Therefore, it should be regarded as the natural extension of the deterministic Euler scheme, in contrast to the so-called Euler-Maruyama scheme, which is merely the "transcription" of the deterministic Euler scheme, but not its equivalent.
An adaptive time-step scheme is proposed to ensure the consistency of the model implementation with the WMC requirements.The time-step selection algorithm is limited not only by the condition imposed by the spatial scale of gradients, but also by taking into account the scale of the width of maxima and minima of the diffusion coefficient, where the former criterium fails.It is shown that this algorithm ensures  The initial distribution ∝ ρ (black) and the final distributions obtained using the A interpolation scheme (red) and the D interpolation scheme (blue).Dashed lines show the limit of 3 standard deviation around the initial distribution.
region between σ = 0.9 and the lowest boundary (see Figure 9).Figure 10 shows the result of an experiment where 575 the WMC compliance of schemes A and D was tested with the lower limit for ∆t changed to ∆t min = 10 −5 for D. It can be observed that, for the D scheme, strong fluctuations are still present in the same region where the required timestep exceeds the lower limit.This is likely to be caused by the 580 occurrence of strong gradients that can be even larger than in the 'peaked' case, near points with extremely small values of K.In these cases, the A scheme interpolates with a smoother function which reduces the problem.that the error is within an acceptable range also at the reflecting boundaries.However, in case of isolated maxima, this scheme may fail.The implementation of a non-local algorithm, which evaluates t in two additional points, is proposed in order to solve the problem.
Two numerical interpolation and derivation schemes are implemented and tested.The first is based on the linear interpolation of K, and it is presented in two versions: one (D) keeps a constant first order derivative between two grid points, while the other (D ) uses linearly interpolated derivatives in the same interval.The second scheme (A) is based on a modified Akima (1991) interpolation algorithm with a local algorithm that ensures the positivity of the interpolating function (Fischer et al., 1991).
It is found that although the method D uses derivatives of higher order of approximation, it creates a local inconsistency between the linearly interpolated function and its derivatives and prevents the model from fulfilling the WMC.The other two schemes (D and A) both satisfy the WMC, but extremely peaked profiles of K may require the use of the non-local time-step selection algorithm.
A test with a preliminary implementation of the fully 3-D model (IL-GLOBO) shows that, while the A scheme displays a correct behaviour, the D interpolation scheme requires an extremely strong reduction of the integration time step that prevents the WMC to be satisfied in reasonable time.and installation on a variety of systems.The code is developed in a modular way, permitting the easy improvement of physical and numerical schemes.
The GLOBO model is available upon agreement with the CNR-ISAC Dynamic Meteorology Group (contact: p.malguzzi@isac.cnr.it).

Fig. 2 .
Fig.2.Values of integration time-step ∆t for the diffusivity profile shown by the red curve.The green line shows the contribution of Eq. (19), the blue line the contribution of Eq. (20), and the black line the combined condition (Eq.21, with ∆T = 432 s and CT = 0.01).

Figure 2 .
Figure 2. Values of integration time step t for the diffusivity profile by the red curve.The green line shows the contribution of Eq. (19), the blue line the contribution of Eq. (20) and the black line the combined condition (Eq.21, with T = 432 s and C T = 0.01).

Fig. 3 .
Fig.3.Average GLOBO profiles of ρ (green symbols) and φ/g (blue symbols) as a function of vertical coordinate σ, and their analytical fits (Eq.23 and Eq. 24, lines of the same colors).

Figure 3 .
Figure 3. Average GLOBO profiles of ρ (green symbols) and φ/g (blue symbols) as a function of vertical coordinate σ and their analytical fits (Eqs.23 and 24).

Fig. 4 .
Fig. 4. Diffusivity profiles used in the experiments.The symbols represents the data from GLOBO and the lines their fitting function.The 'average' profile is shown in red, while the 'peaked' profile is shown in green.The functional form of both profiles is described by Eq. (26).

Figure 4 .
Figure 4. Diffusivity profiles used in the experiments.The symbols represents the data from GLOBO and the lines, their fitting function.The "average" profile is shown in red, while the "peaked" profile is shown in green.The functional form of both profiles is described by Eq. (26).

Fig. 5 .
Fig. 5. Dispersion experiment with different choices of parameter CT .Top panel: diffusivity profile (black line) and ∆t profiles for CT = 0.5 (light blue), CT = 0.1 (green), CT = 0.01 (red) and CT = 0.001 (blue).Bottom panel: normalized concentration profiles for different CT (Line colors as in the top panel).
. In the upper panel, the integration time-step profiles of the three simulations and the Akima interpolated diffusion coefficient profile are displayed.The lower panel shows the normalised www.geosci-model-dev.net/7

Fig. 8 .
Fig. 8. Same as in Fig. 7 for experiments with the 'peaked' diffusivity distribution.Results obtained using the local (left) or non-local (right) ∆t selection algorithm.

Figure 7 .
Figure 7. Experiments with the sampled "average" diffusivity distribution for the interpolation algorithms D (blue), D (green) and A (red).Top panel: diffusivity profile as interpolated by A (black) and t profiles for the different interpolation settings.Bottom panel: normalised final concentration and expected distribution (black).

Fig. 8 .
Fig. 8. Same as in Fig. 7 for experiments with the 'peaked' diffusivity distribution.Results obtained using the local (left) or non-local (right) ∆t selection algorithm.

Figure 8 .
Figure 8. Same as in Fig. 7 for experiments with the "peaked" diffusivity distribution.Results obtained using the local (left) or nonlocal (right) t selection algorithm.

Fig. 9 .Fig. 10 .
Fig. 9. Distribution of ∆t requirement for the condition on the first order derivative (Eq.19) as a function of σ for interpolation algorithms A (left) and D (right).

Figure 9 .
Figure 9. Distribution of t requirement for the conditions of the first-order derivative (Eq.19) as a function of σ for interpolation algorithms A (left) and D (right).

Fig. 9 .
Fig.9.Distribution of ∆t requirement for the condition on the first order derivative (Eq.19) as a function of σ for interpolation algorithms A (left) and D (right).

Fig. 10 .
Fig. 10.Normalized distribution of particles for the 3D experiment.The initial distribution ∝ ρ (black) and the final distributions obtained using the A interpolation scheme (red) and the D interpolation scheme (blue).Dashed lines show the limit of 3 standard deviation around the initial distribution.

Figure 10 .
Figure10.Normalised distribution of particles for the 3-D experiment.The initial distribution ∝ ρ (black) and the final distributions obtained using the A interpolation scheme (red) and the D interpolation scheme (blue).Dashed lines show the limit of 3 standard deviations around the initial distribution.

Table 1 .
RMSE and execution time for different C T .

Table 1 .
RMSE and execution time for different CT .

Table 2 .
Execution time and RMSE for experiments made with the sampled "average" diffusivity distribution and varying interpolation method.
norm σ Fig. 7. Experiments with the sampled 'average' diffusivity distribution for the interpolation algorithms D (blue), D ′ (green) and A (red).Top panel: Diffusivity profile as interpolated by A (black) and ∆t profiles for the different interpolation settings.Bottom panel: Normalized final concentration and expected distribution (black).

Table 2 .
Execution time and RMSE for experiments made with the sampled 'averaged' diffusivity distribution, varying interpolation method.

Table 3 .
Execution time and RMSE for experiments made with the 'peaked' diffusivity distribution, varying interpolation method and ∆t algorithm.

Table 3 .
Execution time and RMSE for experiments made with the "peaked" distribution, varying interpolation method and t selection algorithm.

Table 2 .
Execution time and RMSE for experiments made with the sampled 'averaged' diffusivity distribution, varying interpolation method.

Table 3 .
Execution time and RMSE for experiments made with the 'peaked' diffusivity distribution, varying interpolation method and ∆t algorithm.