Eulerian Dispersion Model Basd on M. Galperin's Advection Printer-friendly Version Interactive Discussion Construction of an Eulerian Atmospheric Dispersion Model Based on the Advection Algorithm of M. Galperin: Dynamic Cores V.4 and 5 of Silam V.5.5 Gmdd Eulerian Dispersion Model Basd on M. Galperi

The paper presents dynamic cores v.4 and v.5 of the System for Integrated modeLling of Atmospheric coMposition SILAM v.5.5 based on the advection algorithm of Michael Galperin. This advection routine, so far weakly presented in international literature, is non-diffusive, positively defined, stable with regard to Courant number significantly 5 above one, and very efficient computationally. For the first time, we present a rigorous description of its original version, along with several updates that improve its monotonicity and allow applications to long-living species in conditions of complex atmospheric flows. The other extension allows the scheme application to dynamics of aerosol spectra. The scheme is accompanied with the previously developed vertical 10 diffusion algorithm, which encapsulates the dry deposition process as a boundary condition. Connection to chemical transformation modules is outlined, accounting for the specifics of transport scheme. Quality of the advection routine is evaluated using a large set of tests. The original approach has been previously compared with several classic algorithms widely used in 15 operational models. The basic tests were repeated for the updated scheme, along with demanding global 2-D tests recently suggested in literature, which allowed positioning the scheme with regard to sophisticated state-of-the-art approaches. The model performance appeared close to the top of the list with very modest computational costs.


Introduction
One of the key problems in atmospheric composition modelling is the accuracy and reliability of numerical schemes.A less appreciated but important issue is the consistency of the approaches applied in different modules of the modelling system.Usually, model construction follows process-wise split (Yanenko, 1971;Marchuk, 1995;Seinfeld and Pandis, 2006), thus considering separately the advection and diffusion algorithms, physical and chemical transformations, and dry and wet deposition.In practical Introduction

Conclusions References
Tables Figures

Back Close
Full model developments, features of the transport algorithms, first of all, advection scheme, largely shape-up the model and determine its area of application.

Advection schemes
There are numerous types of advection schemes currently employed in atmospheric dispersion models.Two major categories refer to Lagrangian or Eulerian treatment of tracers: as small-size masses (Lagrangian particles) or as the concentration fields discretised in a prescribed grid.The Eulerian schemes, the primary subject of this paper, can be divided to flux-form finite volume, semi-Lagrangian, or expansion-function schemes.The expansion-function schemes approximate the solution with a given set of basis functions and, in turn, can be divided to spectral, pseudospectral and finiteelement approaches.Some classic schemes are also based on finite-difference approximations of the advective term of the transport equation.The basic principles of all these approaches were formulated several decades ago and, with certain modifications, are still in use.Many modern schemes combine several approaches.The large diversity of the advection algorithms is explained by a long list of requirements to such schemes.The most important ones are: positive definition, minimal numerical diffusion, limited non-monotonicity, stability with regard to high Courant number (the number of the model grid cells passed within one advection time step), small phase error, local and global mass conservation, and high numerical efficiency.Some of these requirements contradict to each other.For example, numerical diffusion "blurs" the resulting patterns but also makes them smoother, thus improving the monotonicity.
The finite-difference schemes involve direct discretization of the dispersion equation and various interpolation functions to evaluate derivatives of the concentration fields (see reviews of Richtmyer, 1962;Leith, 1965;Roach, 1980); specific examples are, for instance, (Russell and Lerner, 1981;Van Leer, 1974, 1977, 1979).These schemes, being once popular, usually suffer from large numerical diffusion and limited stability, which sets stringent limitations to the Courant number usually requiring it to be substantially less than one.Therefore, the interest has gradually shifted towards flux, finite-Figures

Back Close
Full element, and semi-Lagrangian schemes, where unconditional stability can be ensured also in explicitly formulated algorithms (Leonard, 2002).The flux schemes represent the transport via mass fluxes at the grid cell borders, which are calculated from concentrations in the neighbouring cells and wind characteristics.They are inherently mass conservative and have become popular in atmospheric models (Kukkonen et al., 2012).Probably the most widely used flux-type scheme is the one made by A. Bott (Bott, 1989(Bott, , 1992(Bott, , 1993)), especially if one would count the numerous Bott-type schemes (see examples in Syrakov, 1996;Syrakov andGalperin, 1997, 2000;Walcek and Aleksic, 1998;Walcek, 2000;Yamartino, 1993), which are based on the same principle but involve different approximation functions, monotonicity and normalization procedures, etc.
The semi-Lagrangian schemes have been among the most-widely used approaches for decades, with numerous algorithms using its basic concept (Crowley, 1968;Egan and Mahoney, 1972;Pedersen and Prahm, 1974;Pepper and Long, 1978;Prather, 1986;Smolarkiewicz, 1982;Staniforth andCote, 1991 andreferences therein, Lowe et al., 2003;Sofiev, 2000, etc).In forward form, these schemes consider the transport of mass starting from the grid mesh points (departure points) and calculate the masses at the grid points closest to the arrival point.Backward algorithms start from arrival grid points and find the grid points near the departure point.The schemes can be based on tracking either grid points or grid cells along their trajectories.The gridpoint-based schemes are not inherently mass-conserving, whereas the volume-based schemes achieve mass conservation by integrating the mass over the departure volume.They can sometimes be described as a combination of finite-volume and semi-Lagrangian methods (Lin andRood, 1996, 1997).Stability of these schemes can be ensured for a wide range of Courant numbers (Leonard, 2002).A general review can be found in Introduction

References
Tables Figures

Back Close
Full but not widely used today.It is based on solving the transport equation in spectral space.Among comparatively recent developments, the schemes for non-Cartesian or unstructured grids can be mentioned (Dendy et al., 2002;Hanert et al., 2004), as well as dynamic-grid algorithms (Jablonowski, 2004;Jablonowski et al., 2006;Lagzi et al., 2004;Staniforth and Cote, 1991 and references therein).The approaches are often combined with semi-Lagrangian or, rarer, with finite-volume or finite-difference advection schemes.Dynamic grids are mainly used for solving the problems in presence of sharp gradients of the computed fields and a wide range of spatial and temporal scales of input forcing.For such problems, advantages of more accurate computations in the sub-domains, which require high resolution, outweigh the extra errors introduced by repeated reprojection of the concentration fields, as well as additional costs of the grid transformation.
One of the main problems of the existing schemes is substantial numerical diffusion originating from the finite-step discretization along space and time.Seemingly inevitable in Eulerian context, this phenomenon, however, does not exist in Lagrangian advection schemes, which do not contain explicit discretization of particle movement.Lagrangian domain is a continuous space rather than a set of pre-defined grid meshes and the position of the particles can be calculated precisely.As a result, numerical diffusion of purely Lagrangian schemes is always zero -at a cost of strongly nonmonotonous concentration fields due to limited spatial representativeness of a single Lagrangian particle and the limited number of particles.
One of ways to reduce the diffusivity of an Eulerian scheme is to store additional prognostic variables to describe the state of each grid cell.It allows to add extra conservation equation for, e.g., first-or higher-order moments (Egan and Mahoney, 1972;Prather, 1986), thus preserving more features of the concentration field.In a series of works, Michael Galperin developed a semi-Lagrangian scheme that was fully nondiffusive, positively defined, and very efficient computationally (Galperin et al., 1994(Galperin et al., , 1995(Galperin et al., , 1997;;Galperin, 1999Galperin, , 2000;;Galperin andSofiev, 1998, 1995).The early ver-Introduction

Conclusions References
Tables Figures

Back Close
Full sion of this scheme applied in the large-scale modelling by (Sofiev, 2000) resembled the "moving-centre" approach widely used in aerosol dynamics models (Kokkola et al., 2008) and shared its characteristic weakness -high non-monotonicity.The later developments substantially reduced it without damaging other features (Galperin, 1999(Galperin, , 2000)).Further development of this scheme will be the subject of the current paper.

Horizontal and vertical diffusion, dry deposition
Diffusion algorithms are less diverse than advection schemes.The physical ground for the diffusion parameterization in numerical models was laid down by (Smagorinsky, 1963), who suggested a formula for grid-scale scalar eddy-diffusivity based on the model resolution and wind speed derivatives, thus connecting the numerical features of the simulations and hydrodynamics.It is widely used by many atmospheric models (Kukkonen et al., 2012).The dry deposition is usually accounted for as a boundary condition to the vertical advection-diffusion equation.Computation of dry deposition for gases practically always follows the electrical analogy, for which one of the first comprehensive parameterizations was suggested by Wesely (1989).Among the extensions of this approach, one was suggested by Sofiev (2002), who combined it with vertical diffusion and connected with the Galperin advection scheme.The algorithm used effective mean diffusion coefficient over thick layers calculated from high-resolution meteorological vertical profiles, the direction also recommended by Venkatram and Pleim (1999).These thick layers were defined using the subgrid information available from the advection scheme, which increased the accuracy of both algorithms.
For aerosol species, the electrical analogy is not correct (Venkatram and Pleim, 1999).Compromising approaches suggested by Slinn (1982), Zhang et al. (2001) and updated by Petroff and Zhang (2010) involve numerous empirical relations, sometimes with thin ground.More rigorous approach unifying the gas and aerosol deposition parameterizations into a single solution was developed by Kouznetsov and Sofiev (2012).Introduction

Conclusions References
Tables Figures

Back Close
Full

Organization of the paper
The current paper describes the Eulerian dynamic core v.4 and v.5 of the System for Integrated modeLling of Atmospheric coMposition SILAM v.5.5, which are both based on the advection scheme of Michael Galperin with several updates.The paper is organised as follows.Section 2 describes the original algorithm of M. Galperin.Section 3 presents the improvements made during its implementation and testing in SILAM, as well as the scheme interconnections to other model parts.
Two sets of standard and advanced model tests are shown in Sect. 5. Finally, discussion in the Sect. 5 includes analysis of the scheme performance in standard and advanced transport tests, as well as comparison with other algorithms.
Following the SILAM standard, all units throughout the paper are the basic SI: [mole] for chemicals, [kg] for aerosols without chemical speciation, [m] for distance and size, [s] for time, etc.The model operates with concentrations, [mass volume −1 ].The dynamic core v.4 is defined in geometric space and uses wind in [m s −1 ], whereas the core v.5 is defined in air-mass space, thus taking wind as an air-mass flux through the corresponding cell border, [kg air m −2 s −1 ].Generic formulations are identical for both core versions.

Basic equations
We consider the forward dispersion equation with the first-order K-theory-based closure for diffusion: where φ is concentration of the pollutant, t is time, E is emission term, x i , i = 1.3 denote the three spatial axes, u i are components of the transport velocity vector along Introduction

Conclusions References
Tables Figures

Back Close
Full these axes, µ i j are components of the turbulent diffusivity tensor, ρ is air density, and ξ represents transformation source and sink processes.The Eq. ( 1) is considered on the time interval t ∈ (t 0 , t 1 ) in the domain {x i } ∈ Ξ = [h 1 , H] × Ω, where the heights h 1 and H are the lower and upper boundaries of the computational domain and Ω is the horizontal computational area with border ∂Ω.The initial conditions are: Boundary conditions depend on type of the simulations.In a general form, they constrain concentration and/or its first derivative: Here the values of α, β, and γ depend on type of the boundary.In particular, dry deposition at the surface x 1 = h 1 is described via α = µ 33 , β = −v d (dry deposition velocity), γ = 0; prescribed concentration ϕ l at the lateral boundaries x 1,2 ∈ ∂Ω implies α = 0, β = 1, γ = ϕ l , etc.A global domain would require periodic longitudinal conditions.

Advection scheme of Michael Galperin
The current section presents, for the first time in international literature, the advection algorithm suggested by Michael Galperin in 2000s as a contribution to Eulerian dynamics of SILAM.The idea of the scheme can be found in a short methodological publication of Galperin (2000, in Russian) and conference materials (Galperin, 1999;Sofiev et al., 2008).
For the 1-D case, let us define the simulation grid as a set of N grid cells i = 1. ..N. Let the coordinate of the centre of the i th grid cell be x i , and coordinates of the cell Introduction

Conclusions References
Tables Figures

Back Close
Full left-and right-hand borders be x i −0.5 and x i +0.5 , respectively.The 1-D cell size is then The advected field is described by a set of rectangular slabs, one for each grid cell i .The slab is represented by the total mass in the cell M i and position of the centre of mass X i , X i ∈ [x i −0.5 , x i +0.5 ]: where k is time step and 5 is distance from the centre of mass X k i to the nearest border of the cell i .The slab therefore has the largest width inside the cell i allowed by position of its mass centre (Fig. 1).At this moment, the non-zero concentration of this slab is fully confined within the cell i .The first moment of this slab is: Advection of the slab does not change its shape within the time step δt but can move it anywhere over the domain or bring outside.In explicit form, it reads: The concentration distribution over the domain at k + 1 step is considered as a sum of the advected slabs over al grid cells 1. ..N: Note that now the slabs can be located anywhere irrespective to their original grid cells.
Projection of the distribution Eq. ( 6) back to the grid concludes the advection step (Fig. 1).For each grid cell i , the new mass and its position are computed from the mass Figures

Back Close
Full and the first moment conservation requirements: x i −0.5 ϕ k+1 (x)dx ( 7) and are stored and used to construct the new slabs at the next advection time step following the definition (Eq.4).
Boundary conditions are defined separately for in-and out-flow.Since Eqs. ( 5) and ( 6) do not limit the slab position by the grid borders, the out-flow boundary is transparent and the mass leaving the domain during the kth step is computed the same way as Eq. ( 7): The original scheme was formulated for zero inflow boundary conditions.Generalization of the above 1-D algorithm to 2-and 3-D spaces is straightforward: slabs become rectangles or cuboids and 1-D integrals are replaced with 2-and 3-D formulations, respectively.An alternative is to perform the 1-D advection along each dimension sequentially.It simplifies the formulations but requires additional measures to ensure conservation of the cross-transport-dimension moments.Namely, the mass M xyz in the cell (x, y, z) is accompanied with its three moment components P x , P y , P z .Introduction

Conclusions References
Tables Figures

Back Close
Full For the time step k, they read: When the advection step is performed along some axis, the moment along it is recomputed as stated in Eq. ( 8).Two other moments are transported in exactly the same manner as the mass itself, thus requiring the repetition of the advection step three times.In practice, the intermediate results of the mass advection can be used, so that only the last step of the grid projection has to be repeated for each moment component.

Construction of the SILAM v.5.5 dynamic cores v.4 and v.5
Within this section, we present the updates and extensions to the Galperin advection scheme, organization of the vertical advection-diffusion interactions, and connections to emission and chemical transformation modules.Particular attention is paid to utilising the strengths of the advection routine, i.e. the sub-grid information on the slab position, to increase the overall model accuracy.

Lateral and top boundary conditions
The open boundary for the outgoing masses is kept in all SILAM regional simulations.
The inflow into a limited-area domain is defined via prescribed concentration at the boundary (Eq.3), α = 0, β = 1, γ = ϕ l .The mass coming into the domain during a single time step is equal to:

GMDD Introduction Conclusions References
Tables Figures

Back Close
Full Here ℵ(u) is Heaviside function (= 1 if u > 0, = 0 if u ≤ 0).Assuming locally constant wind we obtain that M in is distributed uniformly inside the slab similar to that of Eq. ( 4).For, e.g., the left-hand-border, the continuous form will read: It is then projected to the calculation grid following Eqs.( 7) and ( 8).
The top boundary follows the same rules as the lateral boundaries.At the surface, the vertical wind component becomes zero, which is equivalent to closure of the domain with regard to advection.
Global-domain calculations require certain care: SILAM operates in longitudelatitude grids, (possibly, rotated with repositioned southern pole and zeroth meridian), i.e. it has singularity points at the poles and a cut along the 180 Vertical motion in the pole cylinders is calculated separately using vertical wind component diagnosed from global non-divergence requirement.Introduction

Conclusions References
Tables Figures

Back Close
Full

Extension of the scheme for complex wind pattern
The original Galperin's scheme Eqs. ( 4)-( 10) has a drawback that was discovered during its tests in real-life complex-terrain applications: it utilises the constant-wind assumption at a spatial scale of a single grid cell (or at a distance of uδt, whatever is larger).This is not correct for regions with complex terrain, in the vicinity of atmospheric fronts, etc.As a result, the scheme tends to accumulate mass at stagnation points where one of the wind components is small.Similar problems were reported by Ghods et al. (2000) for the scheme of Egan and Mahoney (1972).Ghods et al. also suggested the generic principle for solving the problem: increasing the number of points at which the wind is evaluated inside the grid cell.In application to Galperin's scheme, it can be achieved by separate advection of each slab edge instead of advecting the slab as a whole.This allows for shrinking and stretching the slab following the gradient of the velocity field.Formally, this can be written as follows.
Let's again consider the 1-D slab that has been formed according to Eq. ( 4).Its edges are: Let the wind be defined at the borders of the grid cells.The advection step takes the wind velocity at the left and right edges of the slab and transports these edges accordingly.The transport of each edge accounts for the changing velocity within and between the grid cells and integrates it along the path through the whole time step.
If wind is assumed to be linearly interpolated inside the cell, this integration can be made analytically.As a result, the time step modifies the slab width and the new one is formed as a uniform distribution between the new positions of the edges: are the new positions of the slab edges at the end of the time step.This equation replaces the step Eq.(5).

Adaptation of the scheme to aerosol size spectrum computations
The scheme Eqs. ( 4)-( 15) can be adapted to compute the evolution of the particle size spectrum due to aerosol-dynamics processes.This extension requires certain care because not all processes can be presented in a form of advection Eq. ( 1).There is a direct analogy only between the advection in the particle size space and the condensation process (Seinfeld and Pandis, 2006).Nucleation and coagulation, however, cannot be described in such a form as the number of particles is not conserved in such processes.
Let us represent the considered range of particle sizes via a set of bins, which form the grid in the discrete size space.For each bin i , let us denote the total volume of particles belonging to it as V i , volume of individual particle as v i , total number of particles in the bin as N i and a volume-weighted mean radius as R i .Note that: Let us consider two steps: modification of the particle size in each bin, and projection of the new distribution back onto the bins.Also, let's consider the process in the singleparticle-volume space instead of the size space.Then the processes changing the particle size in the i th bin will result in a change of v i : which can be considered as a shift of the whole bin with δv i without changing its shape.Such transformation is analogous to the slab movement Eq. ( 5).Replacing Eq. ( 5) with ( 17), one can apply Galperin's advection scheme in the discrete aerosol particle volume space, no matter what aerosol dynamic processes are taken into consideration.Introduction

Conclusions References
Tables Figures

Back Close
Full Aerosol spectra sometimes exhibit very sharp gradients due to, e.g., evolution of the spectrum due to condensation or nucleation bursts (Seinfeld and Pandis, 2006).In such cases, sectional approach often lacks resolution resulting in too crude representation of the spectrum even accounting for the sub-grid information on the slab positions.A marginal improvement can be obtained by allowing non-homogeneous number distributions inside the bins, thus refining the utilization of the sub-grid information by accounting for values in the neighbouring cells.The whole spectrum is then split to ranges with dN/dR > 0 (range 1) and dN/dR < 0 (range 2).For the first range, the number distribution with regard to volume is approximated by N(v) = const.For the second range, N(v) ∼ v −1/2 .These two distributions follow the monotonicity of N(R), thus smoothing the distribution shape and allowing analytical solution of the advection Eqs. ( 5)-( 8), if written in notations Eqs. ( 16) and ( 17).

Shape preservation for "puff over background" patterns
Along with strong scores of the scheme in most tests (discussed below), one can point out an undesirable feature visible in conditions of strong sharp plumes over a low background level (Fig. 2).The scheme tends to artificially sharpen the plume edges and to gradually redistribute the background mass in the vicinity of the plume towards it (Fig. 2, upper row).This feature, albeit comparatively small (the error in Fig. 2 took 300 steps to evolve), leads to noticeable distortions in the vicinity of strong point sources or around dense plumes.
The reason for the feature follows from Eqs. ( 7) and ( 8): if the mass is concentrated at one of the grid cell sides, the resulting centre of mass location becomes insensitive to the low-mass part of the cell, i.e. the small mass that appears there from the neighbouring cell is just added to the big slab with little effect on its position.
The cheapest, albeit not rigorous, way to confront the effect is to compensate it via additional filtering.In SILAM, the filter is made as a three-point 1-D symmetric moving average, which redistributes a fraction ψ of the central cell mass to each of Introduction

Conclusions References
Tables Figures

Back Close
Full The optimal value for ψ max was empirically found to be 0.01 for real-life applications where it additionally compensates the absence of horizontal diffusion presently not included in SILAM dynamics.For tests it is decreased down to 0.001 to (almost) preserve the zero diffusivity of the scheme.

Vertical axis: combined advection, diffusion, and dry deposition
The vertical diffusion algorithm is implemented according to the extended resistance analogy (Sofiev, 2002), which considers air column as a sequence of thick layers.Vertical slabs within these layers are controlled by the same 1-D advection, which is performed in absolute coordinates -either z-or p-depending on the vertical (height above the surface or hybrid).
Diffusional exchange coefficient between the neighbouring layers is taken as the inverse of aerodynamic resistance between centres of mass within the layers.The effective thickness of the layers is taken proportional to pressure drop across the layer, which assures equilibration of mixing ratios due to diffusion.Settling of particles is included into advection for all layers except for the first one, where the exchange with the surface is treated by the dry deposition scheme.It takes the dry deposition velocity at the height of centre of mass in the lowest layer and follows the approach of Kouznetsov and Sofiev (2012) for incorporating the sedimentation velocity.Since the diffusion is calculated between the vertical positions of the mass centres, the subgrid information provided by the vertical advection is taken into account.Introduction

Conclusions References
Tables Figures

Back Close
Full

Connection with other model components
Interfaces between the dynamic core and other SILAM v.5.5 modules utilise the positions of mass centres to increase the overall system accuracy.Emission-to-dispersion interface distinguishes between the point sources and area sources.For point sources, the mass is added to the corresponding grid cell accounting for its centre-mass position and affecting both mass and moment in the cell following Eqs.( 7) and ( 8).
For area sources, the approach depends on the source grid.If it is the same as the computational one, the mass centre is put to the middle of the cell (no extra information can be obtained).If the grids are different, the source is reprojected to the computation grid using the mass-and the first moment-conserving algorithm.
Meteo-to-dispersion interface has to meet the stringent requirements of the advection procedure regarding the non-divergent air flux.Indeed, as the scheme has zero numerical diffusivity, even slight inconsistency between the wind fields and advection will create local mass conservation problems, which are not smeared out by numerical diffusion.Therefore, the wind fields are pre-processed to ensure exact satisfaction of the continuity equation.
The methodology for boundary layer mixing parameterization implemented in SILAM since v.4.x is described in (Sofiev et al., 2010).For the free troposphere and the stratosphere, eddy diffusivity is diagnosed following (ECMWF 2008).The vertical exchange coefficient is expressed as a function of gradient Richardson number and wind shear calculated between neighbouring meteorological levels.The characteristic length scale is prescribed as a universal function of height that equals to metric height near the surface and asymptotically approaches 150 m above.Despite having little physics behind, the approach provides realistic wind fields in the stratosphere with the ECMWF IFS model, and results in realistic patterns of age of air in the stratosphere calculated by SILAM.Introduction

Conclusions References
Tables Figures

Back Close
Full Chemical module interface is implemented in a very simple manner: the mass centres are not affected by the transformations.Chemical module deals exclusively with concetrations in the grid cells, which are used to calculate new moments after the transformation routine finishes.

Standard tests
A set of basic tests and comparison with some classical approaches has been presented by Galperin (1999) for the original scheme.It was also tested by Petrova et al. (2008), along with Bott, Holmgren, and TRAP (Bott-type) schemes.The tests demonstrated the scheme quality and positioned it along with several classic approaches.
More sophisticated tests added during SILAM v.4.x and v.5.x implementation (dynamic core v.4) included puff-over-background (Fig. 2), high-Courant number wind (Fig. 3), constant-level background field in eight vortices with stagnation points (Fig. 4), and rotation tests for various shapes (Fig. 5).The scheme stays stable at arbitrarily high Courant numbers and maintains very low distortions of the shapes (the L 2 norm of the error varies from 0.1 up to 3.8 % of the initial-shape norm -for the most challenging task in Fig. 5).The test with eight vortices was essentially failed by the original scheme (Fig. 4a) but the improvements Eqs. ( 14) and ( 15) resolved the problem (Fig. 4b).conditions and a prescribed non-divergent velocity field, with the analytical solution at the end of integration being identical to the initial pattern.

Test setup
For the sake of completeness, below we provide a brief overview of the tests, whose detailed description can be found in Lauritzen et al. (2012).The velocity field is defined on a sphere via the stream function dependent on longitude λ, latitude θ and time t: where R is the radius of the sphere, and T is the simulated time interval -the spatial and temporal scales of the exercise, respectively (e.g., Earth radius and 12 h).The velocity components u and v are then given by: The velocity field Eqs. ( 19) and ( 20) consists of deformation and rotation and reconstructs the initial concentration field at t = T , providing the exact solution φ(t = 0) = φ(t = T ).

GMDD Introduction Conclusions References
Tables Figures

Back Close
Full correspond to the model time step of T/12 = 1 h, T/24 = 30 min, and T/72 = 5 min), total 18 runs for each initial pattern.Examples of the most challenging runs with slotted cylinders in the middle (t = T/2) and at the end of the simulations t = T are shown in Figs.7 and 8, respectively.The corresponding error fields are collected in Fig. 9 as decimal logarithms of the absolute difference between the corresponding field in Fig. 8 and the slotted-cylinder initial shape of Fig. 6.

Performance metrics
Deviation of the resulting field ϕ T = ϕ(λ, θ, t = T ) from the initial shape, which is also the exact solution ϕ 0 = ϕ(λ, θ, t = 0), was considered in three spaces: L 2 , L ∞ , L 1 .The corresponding distance metrics are defined as follows: where S[•] is an area-weighted sum over latitude and longitude.The values of these three metrics for all model runs are presented in Fig. 10.Advection should also keep the local ratio of the tracer's concentrations.Such ratio between "cosine bells" and "correlated cosine bells" was calculated at t = T/2 and t = T .Since these initial patterns are related by Eq. ( 21), the concentration fields during the tests should maintain the same relation.The scatter plots of the concentrations in these two tests give an indication on how the ratio is kept.Ideal advection would keep all points on a line given by Eq. ( 21).The results of the tests for t = T/2 are shown in Fig. 11.

Efficiency of Galperin advection scheme
Evaluation of the scheme efficiency is always very difficult as it depends on computer, parallelization, compiler options, etc.Nevertheless, some basic characteristics

Conclusions References
Tables Figures

Back Close
Full of the scheme have been deduced from comparison of the simple cases for classical schemes (Galperin, 2000).It was shown to be about twice faster than, e.g., Bott scheme.
The run with 0.75 • resolution and 120 time steps (took 47 s) can be related to performance of HEL and CSLAM schemes, which were tested against the same collection by (Kaas et al., 2013).With all ambiguity of the runtime parameter, this run took about 200 s for HEL and 300 s for CSLAM, i.e. about 4 and 6 times longer than in the SILAM v.5.5, dynamic core v.5.

Discussion
The presented SILAM dynamic cores v.4 and v.5 is based on semi-Lagrangian advection scheme of M. Galperin with subgrid information are available through the centre of mass position.It poses certain challenges in implementation but demonstrates strong performance with low computational costs.In particular, by attributing the release from a point source to its actual location one can reduce the impact of the common problem of Eulerian models: point release is immediately diluted over the model grid cell.This substantially improves the transport but does not solve the problem completely: (i) the chemical module still receives the diluted plume concentration, (ii) the slab size in case of the source near the centre of the grid cell will still be as large as the grid cell itself.
An accurate solution would be the plume-in-grid or similar approaches, which is being built in SILAM.Another example of the sub-grid information usage is utilisation of full

Conclusions References
Tables Figures

Back Close
Full meteorological vertical resolution to calculate effective values of meteo variables for thick dispersion layers (Sofiev, 2002).The model dynamic core can operate at any Courant number (Fig. 3).Its time step is limited not by grid cell size but by a spatial scale of the wind-shear field, i.e. has to satisfy much less restrictive Lipchitz criterion, which relates spatial and temporal truncation errors (Pudykiewicz et al., 1985).It follows from the advection Eq. ( 5) and new-concentration computation Eq. ( 6), which do not restrict new positions of the slabs: they can find themselves anywhere in the grid or outside it after the time step is made.
SILAM heavily relies on such features of Galperin's scheme as mass conservation and accountability: the scheme provides complete mass budget including transport across the domain boundaries.In particular, nesting of the calculations is straightforward and does not need the relaxation buffer at the edges of the inner domain: the inflow through the boundaries is described by the same slabs as the main advection.The scheme is also shape-preserving -in the sense this term is used by L14, -i.e. it does not result in unphysical solutions, such as negative mixing ratio.Some distortions are still possible (Fig. 2), which can be reduced by filtering described in Sect.3.4, Eq. (18).

Standard advection tests
Evaluating the Galperin's scheme with the simple tests (Figs.2-5), one can point out the known issues of the classical schemes resolved in Galperin's approach: high-order algorithms suffer from numerical diffusion, oscillations at sharp gradients (require special efforts for limiting their amplitude), high computational costs and stringent limits to Courant number.None of these affect the Galperin scheme.
The main issue noticed during the implementation of the original scheme Eqs. ( 4)-( 8) was the unrealistically high concentrations near the wind stagnation points.Thus, the concentration pattern at the test Fig. 4a resembles the situation of divergent wind field.However, it is not the case: the 2-D wind pattern is strictly solenoidal.The actual reason is insufficient resolution of the advection grid: one centre of mass point is not enough Introduction

Conclusions References
Tables Figures

Back Close
Full if spatial scale of the wind variation is comparable with the grid cell size.Tracking the edges of the slab rather than its centre resolves the problem (Fig. 4b).
The other challenging tasks for Galperin's algorithm were those with smooth background and soft gradients, a frequent issue for semi-Lagrangian schemes, which is easily handled by more diffusive approaches.This feature was also visible in the tests of (Petrova et al., 2008), where the scheme noticeably distorts the Gaussian and conical plumes.
For the puff-over-background pattern, the scheme makes a single low-mass dip in the vicinity of the puff, which receives this mass (Fig. 2).From formal point of view, the scheme does not conserve the higher moments inside the grid cell, which becomes a problem when the pattern changes at a spatial scale shorter than a grid cell size.The shape-preserving filtering step Eq. ( 18), Fig. 2, may be advised despite it has no rigorous ground and, as in L14 evaluation of other schemes, somewhat damages the formal quality scores.In general, such test is among the harshest for advection routines, which develop oscillations, smooth the gradient, and/or produce inhomogeneous solution for both the plain and the puff.

Advanced advection tests
The application of the scheme to highly challenging tests of Lauritzen et al. (2012) allowed its evaluation in a tough 2-D case and comparison with state-of-the-art schemes evaluated by L14 and Kaas et al. (2013).
Performing these tests with different spatial and temporal resolutions, as well as Courant numbers, suggested that the scheme has an "optimal" Courant number for each spatial resolution where the error metrics reach their minimum, so that the increase of temporal resolution is not beneficial.Indeed, in Fig. 10 the low-Courant further runs are by no means the most accurate.This is not surprising: for an ideal scheme, increasing the grid resolution and reducing the time step should both lead to gradual convergence of the algorithm, i.e. the error metrics should reduce.For real schemes, higher temporal resolution competes with accumulation of the scheme errors Introduction

Conclusions References
Tables Figures

Back Close
Full with increasing number of steps.Convergence in L14 tests was still solid for all fixed-Courant-number series (Fig. 10) but excessive temporal resolution (specific for each particular grid cell size) was penalised by higher errors.Similarly, the most-accurate representation of the correlated patterns is obtained from the runs with the intermediate Courant numbers (Fig. 11).This seems to be a common feature: the same behaviour was noticed by L14 for several other schemes.High optimal Courant numbers, however, should be taken with care.For L14, the smooth wind fields reduced the dimension-split error and made long time steps particularly beneficial.For more realistic wind fields, it is likely that the splitting problems will more significantly contribute to the error and require shorter time steps.
It is also seen (Fig. 9) that the best performance, in case of near-optimal Courant, is demonstrated by the high-resolution simulations, which have reproduced both the sharp edges of the slotted cylinders, the flat background and the cylinder's top planes.
The scheme demonstrated convergence rate higher than one for all metrics and all tests with smooth initial patterns.Even for the stringiest test with the slotted cylinders, the scheme reached the linear convergence rate in metric L 1 (Fig. 10).

5.3
Other parts of the SILAM v.5 system SILAM v.5.5 does not include horizontal diffusion, largely for historical reasons.It has been estimated that the contribution of resolved wind fluctuations and vertical wind shear are more important (Pekar, 1989).Conversely, some other old works suggested that for global stratospheric transport even the non-diagonal terms of the diffusivity tensor might be of importance (Karol, 1970).This constitutes one of near-future extensions of the model.
The better performance of the advection at Courant number greater than 1 poses significant challenges to implementation of other modules, first of all, chemistry and emission, which should be programmed accordingly.Indeed, introduction of emitted mass once per long time step would result in a broken plume unless the mass is spread downwind over the corresponding distance.Similar problems show up in chemical Introduction

Conclusions References
Tables Figures

Back Close
Full transformation calculations, especially for fast reactions.At present, the actual SILAM applications are performed with Courant close to but smaller than one to avoid such problems.
The original scheme was formulated for the bulk mass of all transported tracers, thus performing the advection step for all species at once: the tracer's mass in the slab definition Eq. ( 4) was the sum of masses of all species.They were transported within the slab like within a single "container".This algorithm is faster than the specieswise advection used in SILAM since v.4.3 (dynamic core v.2) and reduces the number of the moments per dimension down to one regardless the number of tracers.It can also be useful in case of strong chemical binds between the species in coarse-grid and sub-optimal Courant number: as seen from Fig. 11, such runs can have noticeable non-linearity between the tracer concentrations.The bulk advection does not have such problem but instead loses much of its quality if the species have substantially different life times in the atmosphere, are emitted from substantially different sources or otherwise decorrelated in space.Use of the bulk algorithm in such cases can lead to significant errors.

Summary
Current paper presents the dynamic cores v.4 and v.5 of System for Integrated mod-eLling of Atmospheric coMposition SILAM v.5.5, which are both based on the improved advection routine of Michael Galperin combined with separate developments for vertical diffusion and dry deposition.
The corner stone of the advection scheme is the subgrid information on distribution of masses inside the grid cells, which is generated at the emission calculation stage and maintained in a consistent way throughout the whole model, including chemical transformation, deposition, and transport itself.This information, albeit requiring substantial storage for handling, allows for accurate representation of transport.Introduction

Conclusions References
Tables Figures

Back Close
Full The scheme is shown to be particularly efficient for point sources and sharp gradients of the concentration fields, still showing solid performance for smooth patterns.The most challenging task was found to be the puff-over-plain test, where the scheme showed noticeable distortions of the concentration pattern.Shape-preserving filtration efficiently reduces the problem at a cost of non-zero diffusivity of the resulting scheme.
Advanced tests and comparison with state-of-art algorithms confirmed the compromise between the efficiency and accuracy.SILAM performance was fully comparable with much more computationally demanding algorithms, outperforming some of them.
Among the future developments, introduction of physically grounded horizontal diffusion procedure and replacement of the shape-preserving compensation filter with extensions of the core advection algorithm, are probably the most-pressing ones.

Code availability
SILAM is a publicly available model.Our experience shows however that its successful application critically depends on the user's modelling skills and understanding of the model concepts.Therefore, SILAM is available on-request basis from the authors Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |its neighbours.The filter strength ψ is dynamic and depends on significance of the distortion due to mass centre recalculation: in case of Courant number close to any integer value, no compensation is needed.Therefore, we use the following relation: ψ reaches maximum ψ max at C = 0.5 and decreases to zero towards C = 0 and C = 1: Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Performance of Galperin's advection scheme in global spherical domain was assessed with the collection of demanding tests of Lauritzen et al. (2012) -for SILAM v.5.5 dynamic core v.5.The cases are designed to evaluate the accuracy of transport schemes at a wide range of resolutions and Courant numbers.The tests use four sets of initial Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |