Effects of resolution on the relative importance of numerical and physical horizontal diffusion in atmospheric composition modelling

Numerical diffusion induced by advection has a large influence on concentration of substances in atmospheric composition models. At coarse resolution numerical effects dominate, whereas at increasing model resolution a description of physical diffusion is needed. A method to investigate the effects of changing resolution and Courant number is defined here and is applied to the WAF advection scheme (used in BOLCHEM), evidencing a sub-diffusive process. The spread rate from an instantaneous source caused by numerical diffusion is compared to that produced by the physical diffusion necessary to simulate unresolved turbulent motions. The time at which the physical diffusion process overpowers the numerical spread is estimated, and it is shown to reduce as the resolution increases, and to increase with wind velocity.


Introduction
Discretisation of equations used in numerical simulations of atmospheric composition models has two main effects: the numerical solution of the advection equation causes the spread of localised concentrations, which has the effect of smoothing the gradients, referred to as numerical diffusion error (for brevity, numerical diffusion).In addition, discretisation acts as a low-pass filter in the energy spectrum.The unresolved small scales of motion (those smaller that the mesh size) require a parameterization to account for the diffusive effect of the neglected part of the energy spectrum.This effect is referred to as physical diffusion, and is described by the turbulent diffusion coefficient.According Correspondence to: M. D'Isidoro (m.disidoro@isac.cnr.it) to Isichenko (1992), the diffusion coefficient for a cloud of tracer is defined as the limit for large time t of the variance of the position x of the particles underlining that this coefficient is constant with time in a true diffusive process.An exhaustive review of the numerical diffusion issue can be found in Rood (1987), while the problem relating to subgrid turbulence parameterizations is reviewed in Wyngaard (2004).
This work aims to explore the effects of these two aspects on tracer transport and diffusion, described by the equation where U is the non divergent wind field, C is the tracer concentration, and F is the turbulent flux of concentration that needs to be closed in terms of the mean gradient via a turbulent diffusion coefficient.Odman (1997) analysed the numerical diffusion introduced by four different advection schemes, concentrating mainly on numerical aspects.He also compared to few examples of diffusion parameterisation.However, to give a more complete picture, a systematic comparison between numerical and physical diffusion is needed, especially in view of the continually reduced horizontal mesh size in air quality models.As an example, in GEMS (http://gems.ecmwf.int/d/products/raq/) the mesh size of regional models ranges from 0.5 to 0.1 degrees.Orders of 1 km are even attained in other applications (Goncalves et al., 2009;Jimenez-Guerrero et al., 2008).
In the present study, use is made of the BOLCHEM model (Mircea et al., 2008), developed at CNR-ISAC, which is part of the GEMS-RAQ ensemble (Huijnen et al., 2009).BOLCHEM is an on-line coupling of the meteorological mesoscale model BOLAM (Buzzi et al., 2003(Buzzi et al., , 2004) ) with the SAPRC90 gas-phase module (Carter, 1990) and the AERO3 aerosol scheme (Binkowski and Roselle, 2003).The mass conservative Weighted Average Flux (WAF) advection scheme (Billet and Toro, 1997) adopted in the BOLCHEM model is used here to perform a systematic study on the relative importance of numerical and physical diffusion at different spatial resolutions.
The paper is organised as follows.Section 2 discusses a parameterization of diffusion arising from unresolved turbulent motions, as a paradigm of physical diffusion.Section 3 deals with the numerical aspects.Subsequently, Sect. 4 presents some numerical simulations for an idealised case, in order to establish a general method for the evaluation of the numerical vs. physical effects of the diffusion.Finally, some conclusions are drawn.

Sub-grid turbulent diffusion
The term representing the divergence of turbulent fluxes of concentration in Eq. (2), i.e., the diffusive effect of unresolved scale, can be described using a sub-grid turbulent diffusion coefficient D H , by writing ∇ • F = D H ∇ 2 C, where the subscript H emphasises the fact that only horizontal diffusion is considered.The value of D H can be estimated at a given resolution when the properties of turbulence at the scale of the resolution are known.As an example, it can be described in a dynamical fashion by a sub-grid scale model like the one proposed by Smagorinsky (1963).
However, different approaches can be adopted.For example, D H can also be represented statically using known properties of the inertial subrange spectrum and the typical values of turbulent kinetic energy dissipation in the atmosphere.Assuming a Kolmogorov (1941) (K41) spectrum, D H as a function of the wavenumber is given by where C 1 = 0.25C K (C K = 2), C 0 = 6.2 is the Lagrangian structure function constant, and ε can be determined from velocity spectra measured in a variety of flow conditions and experimental arrangements.Following the Tampieri and Maurizi (2007) model, in the boundary layer (z < h), the following are obtained: -unstable conditions (Albertson et al., 1997): -stable conditions (Pahlow et al., 2001): where z−1 = −1 0 +z −1 , in which 0 = 500 m is assumed, and L MO is the Monin-Obukhov length.
Above the boundary layer (z > h) ε = 5 × 10 −5 m 2 s −3 is assumed as being representative of tropospheric data.For model applications, the height h can be determined case by case, using model profiles along with the actual stability.
As an example, in the free-troposphere, for a grid mesh size x = 10 km, the Tampieri and Maurizi (2007) model gives D H = 310 m 2 s −1 .

Discretisation of the advection term
As pointed out in the Introduction, the discretisation of the advection term in Eq. ( 2) produces, in general, the spread of a cloud of tracer advected by the velocity field (see, e.g.Smolarkiewicz, 1984).The nature and magnitude of this spread is a function of the resolution, and depends on the numerical scheme.For simplicity, the spread is referred to as numerical diffusion, regardless of whether the process displays diffusive behaviour or not, i.e. a growth of the cloud size σ proportional to the square root of time.
In order to identify the parameters relevant to the study of numerical diffusion, the horizontal source dimension R and the characteristic wind speed U are selected as scales for length and velocity, respectively, to give to the left-hand side of Eq. ( 2) the non dimensional form: where the prime indicates non-dimensional quantities and operators.The time scale must be defined as T = RU −1 to make Eq.( 6) scale-invariant.This scale represents the time needed for an air parcel to travel across the source.The sodefined non-dimensional time suggests that the solution of the advection equation is the same, if the time is measured in units of the scale T (the time to cross the source) and distance in units of source dimension.To solve Eq. ( 6) numerically, space-time are discretised by x and t, respectively.Combining the non-dimensional grid mesh size x = R −1 x and time step t = U R −1 t, the following set of non-dimensional parameters is defined: -the resolution: -and the Courant number: It is worth noting that ρ ≥ 1 as the actual size of a source is limited by the spatial discretisation.Sources physically smaller than the grid mesh size are represented in numerical simulations by a grid cell, making ρ = 1.Using these parameters, the non-dimensional time is expressed by where N = t t −1 represents the number of integration steps.Any solution of Eq. ( 6) depends on two parameters only, as does the increase in variance σ 2 of the tracer distribution C(x,t) with respect to its initial value, which in nondimensional terms, reads To perform a systematic evaluation of the relative importance of numerical and physical diffusion at different spatial resolutions, in the present study use is made of a mass conservative advection algorithm based on the WAF numerical scheme (Billet and Toro, 1997), which is briefly described below.
The advection part of Eq. ( 2) is solved numerically using an operator splitting approach, namely, by carrying out the computations for each of the space dimensions sequentially.Discretising in i-th direction in space, one obtains where f * is the numerical WAF flux (Hubbard and Nikiforakis, 2003), defined as and φ i are "limiter functions", which generally have the form of an amplification factor applied to the Courant number ν.
In order to eliminate undesired oscillations from the solution, the limiter functions are defined to be functions also of the local flow parameter r = C upwind / C local , which avoids spurious oscillations by adding a numerical dissipation.The limiter functions are defined as where b = max[0,min(2r,1),min(r,2)].

Results and discussion
In order to evaluate the effect of numerical diffusion, several numerical tests were performed.A simplified framework was chosen, with wind field (u,v) = (U,0), in order to single out the effects of sub-grid processes.The initial tracer distribution was chosen with Gaussian-shape and standard deviation σ 0 = R. Aiming to measure the effect of numerical diffusion as a function of the mass distribution resolution, the numerical tests were set up varying the two parameters ρ and ν.The experiments are summarised in Table 1.
Figure 1 shows the increase in normalised puff variance with non-dimensional time for different values of ρ (a) and ν (b).Variations of ρ induce large variations of σ 2 , while varying ν has a slightly weaker impact.Furthermore, in a real numerical simulation, ν is not a fully controllable parameter that can vary in space and time.Only the requirement ν < 1 has to be met.In the following, ν is fixed to 0.6 to simplify the analysis.An example of the experimental setup is given in Fig. 2, showing the time evolution of the source shape for the experiments A 2.5,0.6 and I 0.125,0.6, corresponding to the two extreme resolutions adopted.
For large non-dimensional time, say O(100), it can be assumed that the increase in variance with time can be well approximated by a power law: Fitting Eq. ( 15) to data allows the determination of both slope β and "diffusion coefficient" α.The results are reported in Fig. 3a and b for different ρ.Although β varies, its variation is sufficiently small to justify the direct comparison of α in Fig. 3a.It is worth noting that the representative value of β highlights the sub-diffusive nature of the numerical diffusion process.Moreover, in terms of diffusion coefficient of a normal diffusion process (a process in which variance increases with time), a sub-diffusive process would be represented by a diffusion coefficient that varies with time.The magnitude of this process is mainly driven by resolution, as shown in Fig. 3a, where the sub-diffusion coefficient α varies by orders of magnitude with ρ.This suggests that conditions can be met for physical diffusion (∝ t ) to become dominant, depending on the numerical resolution.To obtain an idea of what such conditions are, a direct comparison with the sub-grid turbulent diffusion described in Sect. 2 was performed.
The increase in non-dimensional variance for the turbulent (Lagrangian) diffusion process is expressed by where D = D(RU ) −1 is the non-dimensional sub-grid diffusion coefficient, which is a function of x.Assuming that the wavenumber is k = π( x) −1 , from Eq. (3) D can be expressed by where γ = (9/2)π −4/3 (0.25C K ) 2 C −1 0 .Note that since the turbulent diffusion process does not depend on U and R, the non-dimensionalisation makes D explicitly dependent on them.Therefore, in order to compare the results of the numerical and turbulent processes, R and U must be selected.Here R, representing the source scale, is given a fixed value, R=12 500 m, representative of the resolution for regional air quality models, where the sub-grid sources are distributed instantaneously over the grid, thus limiting the size of the source itself.The value of U is left to vary in a typical range from 1 to 20 ms −1 .The intersection between σ 2 and σ 2 L ( R,U ) defines the time at which turbulent diffusion starts to dominate over numerical sub-diffusion.Increasing resolution reduces both variances, although, due to the different dependence on ρ, the numerical sub-diffusion coefficient declines more rapidly than the turbulent diffusion coefficient.This makes highresolution simulations more sensitive to turbulence parameterisation.Furthermore, for low wind velocity the time at which turbulent diffusion starts to be dominant arrives earlier.In fact, in the limit U → 0, the numerical effects on diffusion vanish.The non-dimensional time τ at which the size of the puff would be equal for the numerical and turbulent diffusion processes, can be computed combining Eq. ( 15) and Eq. ( 16) by Figure 5a and b reports τ , as a function of ρ (for given R) for fixed ν = 0.6 (varying U ) and fixed U = 5 ms −1 (varying ν), respectively.It is worth noting that wind velocity variations play a major role in determining τ with respect to ν variations, especially for resolutions below O(1).For higher resolutions, however, τ decreases rapidly.
A factor to be taken into account is that for a given grid mesh size x, even when an explicit description of the unresolved energy through a sub-grid scale model is given, the range in which the energy spectrum is not well represented actually extends up to 6 x (Bryan et al., 2003).The energy accounted for in the explicit solution of primitive equations is therefore less than that expected from a K41 in the highwavenumber end of the spectrum.This means that the diffusion coefficient D estimated from Eq. (3), which implies a "perfect" sub-grid model, is probably underestimated, making the results biased towards an underestimation of physical diffusion and, hence, an overestimation of τ .With the present value of the source size R and for grid mesh size between 5 and 20 km (typical of hydrostatic models), the numerical diffusion overpowers the physical one for values of τ approximately between 1 and 10 4 , corresponding to times ranging from few hours to days.Thus the gradients are excessively smoothed and there is no possibility of a more physical description of the diffusion process.
Preliminary experiments conducted for ρ = 3 show that τ drops by order of magnitude.Therefore, for finer grid mesh sizes (∼ 1 km), the role of physical diffusion becomes more and more important.This requires the extension of the present study to values of resolution ρ above 3, moving in to the range attainable by non-hydrostatic models.
Further remark should be made regarding current air quality simulations: the resolution ρ is always about 1, because sources scales are much smaller than the grid mesh size.Therefore in this kind of simulation, numerical diffusion always dominates over the physical one, even for very low wind velocity.

Conclusions
The present article has considered the horizontal spread of a cloud of tracer released instantaneously in a uniform wind field, to study the diffusion induced by the numerical advection scheme WAF in comparison with the sub-grid physical diffusion.The dimensional analysis of the advection equation shows that numerical diffusion depends on two parameters only: resolution (ratio of the source dimension to grid mesh size) and the Courant number.By means of numerical simulations, it is found that the numerical spread induced by WAF is sub-diffusive, i.e., the concentration distribution variance grows with a power of time less than unity.The time at which physical diffusion starts to be larger than numerical diffusion was computed in different configurations and found to decrease as the resolution increases (see Fig. 5).Moreover, it was found to be further reduced in low wind conditions and large Courant number.
Such findings point to the conclusion that for coarse resolution, there is no possibility of correctly describing diffusive processes.The dominating numerical diffusion induced by WAF is larger and displays a "wrong" (sub-diffusive) behaviour.Other schemes may have qualitatively different behaviour, constituting grounds for extending this kind of analysis.

Table 1 .
Summary of the numerical experiments performed, showing parameters ρ and ν.Experiments are organised in groups (A to I), each characterised by a given ρ.The suffix numbers used in the text indicate the values of ρ and ν, respectively, used for each experiment.