On stable and explicit numerical methods for the advection-diffusion equation

Available online at www.sciencedirect.com Mathematics and Computers in Simulation 79 (2008) 561–570 On stable and explicit numerical methods for the advection–diffusion equation Marcin L. Witeka,b,∗, Joao Teixeirac, Piotr J. Flatau d a University of Warsaw, Institute of Geophysics, ul. Pasteura 7, 02-093 Warszawa, Poland b University of Warsaw, Interdisciplinary Centre of Mathematical and Computational Modelling, ul. Pawi´nskiego 5a, 02-106 Warszawa, Poland c NATO Undersea Research Centre, Viale San Bartolomeo 400, 19138 La Spezia, Italy d Scripps Institution of Oceanography, University of California, San Diego, CA 92093, United States Received 7 May 2007; received in revised form 29 February 2008; accepted 1 March 2008 Available online 12 March 2008 Abstract In this paper two stable and explicit numerical methods to integrate the one-dimensional (1D) advection–diffusion equation are presented. These schemes are stable by design and follow the main general concept behind the semi-Lagrangian method by constructing a virtual grid where the explicit method becomes stable. It is shown that the new schemes compare well with analytic solutions and are often more accurate than implicit schemes. In particular, the diffusion-only case is explored in some detail. The error produced by the stable and explicit method is a function of the ratio between the standard deviationσ 0 of the initial Gaussian state and the characteristic virtual grid distance/Delta1S. Larger values of this ratio lead to very accurate results when compared to implicit methods, while lower values lead to less accuracy. It is shown that the σ0//Delta1S ratio is also significant in the advection–diffusion problem: it determines the maximum error generated by new methods, obtained with a certain combination of the advection and diffusion values. In addition, the error becomes smaller when the problem becomes more advective or more diffusive. © 2008 IMACS. Published by Elsevier B.V . All rights reserved. Keywords: Advection–diffusion equation; Numerical stability; Explicit numerical methods 1. Introduction Numerical instability is a major problem for operational weather and climate prediction models which also afflicts several other meteorological and environmental applications [4,6,7]. Stability problems are often associated with the non-linear advection and diffusion components of the equations that form the basis of many of these models and applications. A well-known and important example is the non-linear vertical diffusion equation due to the parameterization of subgrid scale boundary layer turbulence in the atmosphere and in the ocean [1,3,17]. Vertical diffusion algorithms are constrained by the large (in stability terms) time steps imposed by numerical stability considerations related to the large-scale (mostly two-dimensional) transport. This implies that non-linear vertical diffusion algorithms are often utilized with time steps many times larger than required for stable explicit algorithms. ∗ Corresponding author at: University of Warsaw, Institute of Geophysics, ul. Pasteura 7, 02-093, Warszawa, Poland. E-mail address: mwit@igf.fuw.edu.pl (M.L. Witek). 0378-4754/$32.00 © 2008 IMACS. Published by Elsevier B.V . All rights reserved. doi:10.1016/j.matcom.2008.03.001

562 M.L. Witek et al. / Mathematics and Computers in Simulation 79 (2008) 561–570 There are also many situations when both advection and diffusion terms are important and a transport equation needs to be solved as a whole. In the context of the vertical turbulent transport equations in atmospheric and ocean models there are many examples of this type of situation. These include the one-dimensional (1D) vertical advection–diffusion of general atmospheric constituents (e.g. aerosols, chemical species) in the boundary layer as well as recent eddy- diffusivity/mass-flux parameterizations of turbulence and convection[12,14], and are basically the main motivation of the research presented in this manuscript. For the large-scale part of weather and climate prediction models there has been a widespread use of semi-Lagrangian (SL) methods during the last decade or so, in order to tackle the numerical stability problems associated with large-scale advection. The recent use of SL methods, which are explicit methods where a virtual grid is constructed with grid-sizes which are as large as necessary to maintain stability for an explicit algorithm, has been pioneered by Robert [10,11]. Ritchie et al. [9] discuss the way in which the SL scheme was implemented in the European Centre for Medium-Range Weather Forecasts (ECMWF) global numerical weather prediction model. Studies of the SL method in the context of climate models can be found, for example, in Williamson[21], Williamson and Olson [22] and Williamson et al. [23]. There are today several papers discussing the SL technique in atmospheric models and detailed reviews can be found for example in Staniforth and Cˆot´e [15] and Smolarkiewicz and Pudykiewicz [13]. The vertical diffusion equation used for the parameterization of turbulence is typically solved in weather and climate prediction models using an (quasi) implicit method where the variables themselves are considered implicitly in time, while the eddy-diffusivity coefficients are considered explicitly. This method however, does not solve the numerical stability problem due to the strong non-linearities present in the equations and the large time steps that need to be used, often imposed by efficiency constraints on the large-scale part of the model (e.g.[1,3,17]). Due to these stability problems a scheme usually referred to as “more-than-implicit” or “over-implicit”, which corresponds to an implicitness factor larger than 1 has been used in some weather and climate prediction models (e.g. [5]). In Teixeira [16], hereafter T99, a stable and explicit diffusion solver was introduced. This method follows the main general concept behind the SL scheme, by generating a virtual grid with grid-sizes that are as large as necessary in order to guarantee the stability of the explicit diffusion algorithm. In T99 it was shown that for many situations the new method is more stable and more accurate than the implicit scheme. This method has been recently extended successfully to the reaction-diffusion equations [18] and to the vertical turbulent diffusion equation [19]. It is the main objective of the present work to extend the underlying concept of T99 to the advection–diffusion equation. Section 2 describes the two stable explicit solvers for the one-dimensional transport equation. Performance analysis of both methods is presented in Section 3. Section 4 summarizes our findings. 2. Stable explicit advection–diffusion solvers 2.1. Explicit method The one-dimensional advection–diffusion equation for a generic tracer A has the form ∂A ∂t + u∂A ∂x = ∂ ∂x ( k∂A ∂x ) . (1) Assuming that grid spacing is uniform and that the advection velocityu and the diffusion coefficientk are constant and positive, discretizing the equation explicitly in time, centered in space (for diffusion term) and upwind (for advection term) gives A t+/Delta1t x = k/Delta1t /Delta1x2 At x+/Delta1x+ ( 1 − u/Delta1t /Delta1x− 2 k/Delta1t /Delta1x2 ) At x + ( k/Delta1t /Delta1x2 + u/Delta1t /Delta1x ) At x−/Delta1x (2) or At+/Delta1t x = ˜βAt x+/Delta1x+ (1 − ˜α− 2˜β)At x + (˜β+ ˜α)At x−/Delta1x, (3) where the Courant number ˜α = u/Delta1t//Delta1xand the diffusion number ˜β = k/Delta1t//Delta1x2 are two non-dimensional stability coefficients. In this notation /Delta1t denotes the time step and /Delta1x the grid spacing. The algorithm is stable if ˜α+ 2˜β ≤ 1 (e.g. [2]).

M.L. Witek et al. / Mathematics and Computers in Simulation 79 (2008) 561–570 563 2.2. Stable and explicit unified solver The discretized version of the 1D advection–diffusion equation with constant stability coefficients for the new stable and explicit scheme is At+/Delta1t x = βAt x+/Delta1S+ (1 − α− 2β)At x + (β+ α)At x−/Delta1S, (4) where α=( u/Delta1t//Delta1S) and β=( k/Delta1t//Delta1S2). Note that this equation is analogous to (3) but with /Delta1S instead of /Delta1x. To obtain this Eq. (4), the spatial partial derivatives are approximated by finite-differences along the distance /Delta1S, which does not have to be equal the grid length /Delta1x and can be determined by imposing a fixed stability condition (1− α− 2β)= δ(e.g. [2]). As referred above, the von Neumann stability analysis shows that the scheme is stable if δ∈ [0,1). Following this stability condition, an expression for /Delta1S can be obtained: 1 − u/Delta1t /Delta1S− 2 k/Delta1t /Delta1S2 = δ. (5) This second order equation has two solutions for /Delta1S, one of which is negative and therefore neglected. Equations for /Delta1S and corresponding β have then the form: /Delta1S= u/Delta1t+ √ u2/Delta1t2 + 8k/Delta1t(1 − δ) 2(1 − δ) , (6) β = 4k/Delta1t(1 − δ)2 ( u/Delta1t+ √ u2/Delta1t2 + 8k/Delta1t(1 − δ) )2 . (7) Introducing Eq. (5) into Eq. (4) leads to At+/Delta1t x = βAt x+/Delta1S+ δAt x + (1 − β− δ)At x−/Delta1S. (8) The values of A at the points x + /Delta1S and x − /Delta1S can be found by interpolation from the original grid. It is also convenient to describe β and /Delta1S as functions of the initial non-dimensional numbers ˜α and ˜β: /Delta1S= /Delta1x˜α+ √ ˜α2 + 8˜β(1 − δ) 2(1 − δ) , (9) β = 4˜β(1 − δ)2 ( ˜α+ √ ˜α2 + 8˜β(1 − δ) )2 . (10) The new scheme described by Eqs.(6)–(8) (or (8), (9) and (10)) is stable by design, since it imposes a fixed stability condition, which leads to a new virtual grid where the explicit algorithm becomes stable. The values of the generic variable in this virtual grid are obtained by interpolation from the original grid. The performance of the scheme depends on the ˜α and ˜β numbers, the previously mentionedδparameter and the choice of interpolation method, similarly to SL methods and the T99 diffusion scheme. Note that for the case of no advection ( ˜α = 0), the /Delta1S and β parameters reduce to /Delta1S= √ (2k/Delta1t/(1 − δ)); β = (1/2)(1 − δ) and are therefore equivalent to the values presented originally in T99 for the stable and explicit diffusion solver. The full correspondence with Eq.(9) in T99 is achieved by settingδto 2/3. For the case of no diffusion (˜β = 0), β becomes zero and /Delta1S=( u/Delta1t/(1 − δ)). The numerical scheme (Eq. (8)) reduces to At+/Delta1t x = δAt x + (1 − δ)At x−/Delta1S. This is identical to the SL method with the exception of the factor (1− δ). In such a case, δshould be set to zero and a deviation from this will generate spurious errors. This aspect will be further discussed in Section 3.2. 2.3. Fractional-steps version In this section an alternative way for solving the 1D transport equation using the T99 stable and explicit diffusion scheme is presented. The solution to the advection–diffusion problem is obtained using a fractional-steps approach. A semi-Lagrangian method is used to discretize the advective term and the T99 method is used to solve the diffusion

564 M.L. Witek et al. / Mathematics and Computers in Simulation 79 (2008) 561–570 equation. Assuming that grid spacing is uniform and that the advection velocity u and the diffusion coefficient k are constant the solution can be written as At∗ x = At x−u/Delta1t, (11) At+/Delta1t x = βAt∗ x+/Delta1S+ (1 − 2β)At∗ x + βAt∗ x−/Delta1S, (12) whereβ=( k/Delta1t//Delta1S2). The stability of such a fractional-steps method is ensured by the properties of the semi-Lagrangian scheme and the new diffusion algorithm: they are both stable by design. The /Delta1S parameter in this case is given by /Delta1S= √ (2k/Delta1t/1 − δ), and imposing δ= 2/3 corresponds to the best accuracy of the diffusion scheme, /Delta1S= √ 6k/Delta1t. It is important to distinguish /Delta1S obtained with the fractional-steps method from that obtained with the unified solver. In the latter case both ˜α and ˜β numbers affect the distance /Delta1S and consequently the positions of the new virtual grid points. In a mixed advection–diffusion condition, the same/Delta1S value will result in much stronger diffusion in the case of the fractional-steps method than in the new solver described in Section 2.2. 3. Performance analysis As was mentioned in Section 2.2 the accuracy of the new methods depends inter alia on the interpolation method. The situation is similar to that existing in the semi-Lagrangian class of solvers, where the choice of interpolation method plays an important role. This issue was already investigated in many publications (e.g.[8,13,20]) and therefore will not be explored further in this study. In the performance analysis stage a cubic spline interpolation was used during all simulations in order to derive tracer concentrations in virtual grid points. 3.1. The constant diffusion problem In this section a simple situation is being considered in order to analyze in detail the behavior of the solutions of the new method for a case with constant diffusion and a Gaussian initial state. The advantages of such simplified setting include the existence of a simple analytic solution and that it will allow for a detailed study of the dependence of the error on the few free parameters of the problem. In T99 this configuration has already been discussed, but no analysis as the one being carried out was performed. Since this configuration will be extended in the following section in order to include the advection component, a detailed understanding of its properties when only diffusion is being considered is necessary in order to properly analyze the impact of the advection term. In what follows, the sensitivity of the T99 scheme to two fundamental parameters of the scheme will be tested: (i) the standard deviation of the Gaussian initial state σ 0 and (ii) the distance /Delta1S obtained with the new algorithm. The standard deviation of the initial Gaussian is related to the steepness of the horizontal gradient and is a property of the flow: a larger standard deviation implies a smaller gradient. The distance/Delta1S is an important property of the numerical algorithm and is connected to the diffusion coefficient: a larger /Delta1S implies a larger diffusion coefficient. Fig. 1 shows the evolution (up to 50 time steps) of the integrated L2 norm error (in %) of three numerical schemes: the stable and explicit T99 scheme, and two schemes that are both implicit in time but use finite-element or finite- difference discretizations in space. The three panels in Fig. 1 represent situations with different initial Gaussian standard deviations: σ 0 = 10, 5 and 2.5. In each case the maximum initial amplitude is equal to 10 and the background concentration equal to 1. For all the simulations the time step is 1 s and the grid-size is 1 m. This is the default setting used for simulations presented in this manuscript that leads to stability coefficients for advection and diffusion that are equal to their respective velocity and diffusion values. For this particular case the diffusion coefficient is 5 m 2 s−1 which leads to a stability coefficient of 5. For all simulations presented in Fig. 1, it is assumed that β= 1/6. Results in Fig. 1 confirm that the T99 explicit scheme is stable even for large values of the stability coefficient (in this case 10 times larger than the stability limit). For larger values of the initial standard deviation it also produces better results compared to the other schemes. The error associated with the new stable and explicit method is up to two orders of magnitude smaller than with the implicit methods if the initial standard deviation is higher than 5. However, the situation reverses when the initial standard deviation becomes smaller. Similar results are obtained when the diffusion coefficient is changed while the initial standard deviation is kept constant (results not shown). The change of the diffusion coefficient is reflected in the change of the distance/Delta1S.F o r

M.L. Witek et al. / Mathematics and Computers in Simulation 79 (2008) 561–570 565 Fig. 1. The evolution (up to 50 time steps) of the integrated L2 norm error (in %) of three numerical solvers of the diffusion equation: (i) the T99 stable explicit scheme, (ii) the implicit scheme with finite-element (FE) space discretization and (iii) the implicit scheme with finite-difference(FD) space discretization. The initial state is defined as a Gaussian with maximum amplitude of 10 and three different initial standard deviation values: (a) σ 0 = 10, (b)σ0 = 5 and (c)σ0 = 2.5. The background concentration is equal to 1. a weaker diffusion the results are more accurate then for larger values of diffusion. This discussion hints at a relation between the error and a certain function of the ratio between the initial standard deviation and the distance /Delta1S.I n order to investigate if this is the case, the following simulations were performed: for specific values of the diffusion coefficient multiple simulations were performed with a variety of initial standard deviations. The values for the diffusion coefficients ranged from 1 to 7 whereas the standard deviation ranged from 0.5 to 10.5. For each simulation the L2 norm integrated error from the fifth time step was extracted. Fig. 2 shows the error versus the ratio σ 0//Delta1S for three different values of the diffusion coefficient (1, 4 and 7). All the error values seem to collapse along the same curve in a consistent manner. For a value of the ratio σ0//Delta1S equal to 1 the error is less than 0.1%. For decreasing values of the ratio σ0//Delta1S there is a growth of the error asso- ciated with oscillations in the solution. For larger values of the σ0//Delta1S ratio there is a substantial decrease of the error. The solid line represents a functional fit to the results that could allow for the estimation of the error for cer- Fig. 2. The L2 norm integrated error (in %) after five integration time steps as a function of the ratio between the initial standard deviation and /Delta1S. For diffusion coefficients equal to 1, 4 and 7, the distance /Delta1S equals 2.45, 4.90 and 6.47, respectively. The solid line is a fit in the form of f(x) = exp(−20.3x3 + 30.8x2 − 19.4x + 4.4/x3 + 3.4x2 − 4.0x + 1.3), where x = σ0//Delta1S. The dashed line represents a fit to the error obtained from the linear finite-element implicit scheme fimplicit(x) = (0.9x3 − 5.8x2 + 10.9x + 0.16/x3 + 2.1x2 + 0.58x + 2.0), where again x = σ0//Delta1S.

566 M.L. Witek et al. / Mathematics and Computers in Simulation 79 (2008) 561–570 tain values of the diffusion coefficient and the initial standard deviation without having to perform that particular simulation. Similar simulations to those presented for the new scheme were performed with our reference implicit finite-element and finite-difference methods. In these cases the error is also only a function of the ratioσ0//Delta1S. This is a similar behavior to the one found for the new scheme even though the/Delta1S parameter is not used in the implicit methods. This suggests that /Delta1S could serve as a more general scaling parameter for a diffusion problem, suitable for different numerical approaches. A functional fit to the finite-element method’s results is also presented in Fig. 2, with a dashed line. The results from the finite-difference method (not shown) follow closely the line that corresponds to the finite-element method, with only a small shift towards larger errors. Comparing the two curves from Fig. 2 is clearly visible that for values of the ratio σ 0//Delta1S greater than about 0.5 the new scheme is much more accurate than the implicit finite-element (or finite-difference) method. Encouraging is the fact that the new explicit scheme cannot only be stable for large diffusion coefficients but also surpasses in accuracy the implicit scheme for a wide range of σ 0//Delta1S ratios. 3.2. The constant advection–diffusion problem In this section an analysis of the new advection–diffusion schemes and their accuracy is presented. We concentrate mostly on the solver presented in Section2.2, which solves the transport equation in a unified manner. However, results of the fractional-steps method are also presented. Both schemes are additionally compared to a finite-difference implicit scheme with the advection term discretized centrally in space. As an example, Fig. 3 (upper panel) shows the evolution of a Gaussian with an initial standard deviation of 7 with a uniform advection and diffusion specified by the non-dimensional numbers ˜α = 5.5 and ˜β = 3.5. The δparameter is set to zero. The analytic solution (solid lines) and model results (dashed, dotted and dashed–dotted lines) obtained with the three numerical schemes are displayed every 10 integration time steps. The advection direction is from right to left and the initial Gaussian is centered on grid point 500. The integrated L2 norm errors of the three numerical schemes for up to 50 integration time steps are presented in the lower panel (open circles correspond to the results presented in the upper panel). The accuracy of the new method is higher than that of the implicit finite-difference scheme. However, it produces results that are clearly worse than the fractional-steps method, which gives the best accuracy from all the schemes. The fractional-steps results are hardly Fig. 3. (Upper panel) The evolution of an initial Gaussian (every 10 time steps) as simulated with three different numerical schemes solving the 1D advection–diffusion equation: dashed line—the new stable explicit scheme, dotted line—the finite-difference implicit scheme, dashed–dotted line—the fractional-steps version, solid line—analytic solution. The initial state is defined by a Gaussian with a standard deviation 7 and amplitude 10; the background concentration is 1. The advection velocity is 5.5 m/s and the diffusion is 3.5 m 2/s. The results of the fractional-steps version are hardly distinguishable from the analytic solution. (Lower panel) The time evolution of the integrated L2 norm errors (%) produced by each of the solvers. Open circles correspond to results presented in the upper panel.

M.L. Witek et al. / Mathematics and Computers in Simulation 79 (2008) 561–570 567 Fig. 4. The L2 norm error of the stable explicit advection–diffusion solver after five integration time steps obtained with diffusion and advection values varying in such a way that a/Delta1S is always equal to the initial Gaussian standard deviationσ0. Different lines represent different initial standard deviations: 10, 9, 8, 7 and 6, respectively. The initial amplitude is 10. distinguishable from the analytic solution. This is partially due to the fact that the /Delta1S parameter generated with this scheme (only for the diffusion term) is much smaller than the one produced with the more unified solver, since it does not contain the advection term. As will be shown, the error of the scheme depends on the /Delta1S parameter. Previously it was shown that for the diffusion case the error depends only on the ratio σ 0//Delta1S.I nt h e advection–diffusion case, the /Delta1S parameter is a function of both advection and diffusion values. The same value for /Delta1S can be achieved with a variety of combinations of ˜α and ˜β. Fig. 4 shows the error (after five integration time steps) of the new advection–diffusion scheme as a function of the diffusion coefficient ˜β for varying advection and initial standard deviation values. The Courant number for each diffusion value is such that the/Delta1S number derived from Eq. (9) is equal to the Gaussian’s initial standard deviation σ0. Each of the lines in Fig. 4 represents simulations with different initial standard deviations. The error of the new scheme depends on both the diffusion and advection coefficients, even if they correspond to the same /Delta1S. The error reaches the maximum for some specific configuration of ˜α and ˜β and decreases for all other advection and diffusion values leading to the same/Delta1S. This behavior can be understood by investigating two extremes in an advection–diffusion problem: an advection-dominated case and a diffusion-dominated case. In the case of pure advection the new scheme reduces to a semi-Lagrangian scheme. Therefore, an accurate semi-Lagrangian advection solver will dominate even if some diffusion is present. This can explain the error reduction on the left side of the maxima showed in Fig. 4. The right side decrease, on the other hand, is related to the increasing influence of diffusion over advection, in which case the error approaches small values representative for pure diffusion and described in Section3.1. Therefore, the higher error is obtained for a more ‘mixed’ advection–diffusion situation. However, the maximum error attained by any configuration inFig. 4 seems to be independent of the initial standard deviation values. This result suggests a possible relationship between the maximum error andσ 0//Delta1S, somehow similar to the one presented in Section 3.1. To investigate if this is the case, several simulations were performed with the key parameters being in the following intervals: σ0 = (5, 12), /Delta1S =( σ0 − 3, σ0 + 3), ˜α = (/Delta1S− 3.5,/Delta1 S− 0.5) and ˜β = 0.5/Delta1S(/Delta1S− ˜α). For each of the configurations with a fixed σ0 and /Delta1S, and varying ˜α and ˜β, we found the maximum error produced by the scheme. These errors are presented in Fig. 5 as a function of the σ0//Delta1S. The error values collapse along the same line, which clearly shows that the maximum error produced by the new stable advection–diffusion scheme depends only on the σ0//Delta1S parameter. However, this maximum error is obtained only with specific configurations of the ˜α and ˜β parameters, and the more advective or more diffusive is the problem, the smaller the error produced by the new scheme. The accuracy of the fractional-steps method is much higher than that of the unified solver, as illustrated by the black dotted–dashed line inFig. 5. Errors of the fractional-steps scheme correspond to errors generated by the diffusion scheme

568 M.L. Witek et al. / Mathematics and Computers in Simulation 79 (2008) 561–570 Fig. 5. The maximum L2 norm error (after five time steps) produced by the stable advection–diffusion scheme as a function of the ratio between σ0 and /Delta1S (the values of ˜α and ˜β are constrained in such a way as to generate the same /Delta1S value). The black solid line represents a fit to the results: f(x)=( −4.45x2 + 27.2x − 10.0/x2 + 0.497x − 0.412), where x = σ0/Δ S. The dotted–dashed line is the similarly obtained maximum error but produced with the fractional-steps version. alone. Indeed, when the parametric equation from Fig. 2 that corresponds to the new diffusion scheme is applied to the error presented in Fig. 5, almost a perfect match is achieved (results not shown). The error for the advection solver is marginal compared to the diffusion error and therefore does not substantially affect the total accuracy. It was previously mentioned that the accuracy of the new method depends also on the δ parameter. So far we assumed this parameter to be zero, but in a general case it can take values from 0 to less than 1 to assure stability of the scheme. In a pure diffusion problem, for instance, the most accurate results are obtained with δset to 2/3. The δ parameter determines the influence of the central point in deriving concentrations in the following time step. Below an investigation of how the error of the new unified advection–diffusion scheme depends on theδparameter is presented. Fig. 6 presents the integrated L2 norm error (in %) of the new scheme after five integration time steps. Simulations for different diffusion coefficients ranging from 1 up to 10, and different advection coefficients ranging from 0 to 0.4 Fig. 6. The L2 norm error of the new advection–diffusion scheme after five integration time steps as a function of the diffusion coefficient. Different line shades and line markers denote simulations with fixed Courant number: 0, 0.1, 0.2, 0.3 and 0.4. The solid lines represent simulations withδ=0 and the dashed lines with δ= 2/3. The thick black x-marks indicate the crossing points.

M.L. Witek et al. / Mathematics and Computers in Simulation 79 (2008) 561–570 569 were performed. The initial Gaussian standard deviation is 7 in all cases and the maximum amplitude is 10. In Fig. 6 the solid lines represent simulations where the δparameter was set to 0. The dashed lines are for the simulations with the δparameter equal to 2/3. Several other simulations with the δparameter between 0 and 2/3 were also performed, but the error was always in between the values obtained withδset to 0 or 2/3 (results not shown). These results suggest that only with these two choices of the δparameter can the minimum error be achieved. As can be seen in Fig. 6, for each value of the Courant number the solid and the dashed lines cross each other at some specific value of the diffusion coefficient. This crossing point determines whether the scheme is more accurate with the δparameter set to zero, or to 2/3. If there is no advection (plain and dashed black lines), imposing δ= 2/3 gives better results for all diffusion values. But if the Courant number is 0.1 (plain and dashed gray lines), there is a range of ˜β for which the error is smaller for simulations withδ= 0 instead of 2/3, as indicated by the crossing point of the solid and dashed gray lines. For higher values of advection this range shifts towards larger diffusion values. For ˜α = 0.4 the scheme withδset to 2/3 produces smaller errors only when the diffusion coefficient is larger than 10 (black triangled and dashed-triangled lines). These results indicate that the δparameter should be set to 2/3 only if ˜α ≪ ˜β. Such behavior is understandable when we recall that theδparameter specifies how strong is the influence of the central point of the stencil on the following time step (see Eq. (8)). For a pure advection case, when the new scheme reduces to a semi-Lagrangian method, any incursion of the central point would generate spurious errors. Quantitatively this relation can be expressed as: δ= 0.0 if ˜α> −6.8 × 10 −4 ˜β2 + 4.8 × 10−2 ˜β− 3.1 × 10−2. 4. Summary In this paper a stable and explicit method for integrating the diffusion equation, originally proposed in Teixeira[16], is extended to solve the advection–diffusion equation. The main purpose of this work is to highlight the potential of such an approach for advection–diffusion equations. As a preliminary study of the diffusion equation, a simple case with constant diffusion coefficients was explored in some detail. It was shown that for the diffusion of a normally distributed tracer, the error produced by the stable and explicit method is a function of the two free parameters of the problem. In practice it is a function of the ratio between σ 0 and /Delta1S. Larger values of this ratio lead to very accurate results when compared to implicit methods, while lower values can lead to results that are less accurate than with the implicit method. A physical interpretation is that the scheme is more accurate when there is weak diffusion or the tracer’s gradient is small and less accurate in the case of strong diffusion or when the concentration gradient is substantial. Two stable and explicit schemes to solve the transport equation were proposed: (i) a more ‘unified’ scheme where the advection–diffusion equation is solved with a single numerical algorithm; and (ii) a fractional-steps approach that uses the T99 stable and explicit method for the diffusion and a semi-Lagrangian method for the advection. The unified approach was explored in detail and its results were compared to an implicit finite-difference scheme. It was shown that the error produced by this solver varies depending on advection and diffusion values and theδparameter. However, for a combination of advection and diffusion numbers (that lead to the same/Delta1S) the maximum error depends only on the ratio between σ 0 and /Delta1S. It was argued that the more diffusion– or advection-dominated the problem, the smaller the error generated by the explicit advection–diffusion scheme. Better results were obtained with the δparameter set to zero, except for conditions where diffusion highly dominates over advection, in which caseδshould be set to 2/3 to minimize the error. The fractional-steps version gave better results in all analyzed situations, with errors having similar characteristics to the diffusion-only scheme. References [1] A.C.M. Beljaars, Numerical schemes for parameterizations, in: Proceeding Numerical Methods in Atmospheric Model, ECMWF, Reading, United Kingdom, 2, 1991, pp. 1–42. [2] C.A.J. Fletcher, Computational Techniques for Fluid Dynamics V olume 1 Fundamental and General Techniques, Springer-Verlag, New York, Berlin, Heidelberg, 1991. [3] C. Girard, Y . Delage, Stable schemes for the vertical diffusion in atmospheric circulation models, Mon. Weather Rev. 118 (1990) 737–745. [4] G.J. Haltiner, R.T. Williams, Numerical Prediction and Dynamic Meteorology, John Wiley & Sons, 1980. [5] M. Jarraud, A.J. Simmons, M. Kanamitsu, Development of the High Resolution Model, ECMWF Tech. Memo No. 107, ECMWF, Reading, United Kingdom, 1985. [6] G.I. Marchuk, Numerical Methods in Weather Prediction, Academic Press, New York and London, 1974.

570 M.L. Witek et al. / Mathematics and Computers in Simulation 79 (2008) 561–570 [7] N.A. Phillips, An Example of Non-Linear Computational Instability in: The Atmosphere and the Sea in Motion, Rockefeller Institute Press, New York, 1959, pp. 501–504. [8] R.J. Purser, L.M. Leslie, An efficient interpolation procedure for high order three dimensional semi-Lagrangian models, Mon. Weather Rev. 119 (1994) 2492–2498. [9] H. Ritchie, C. Temperton, A. Simmons, M. Hortal, T. Davies, D. Dent, M. Hamrud, Implementation of the semi-lagrangian method in a high-resolution version of the ECMWF forecast model, Mon. Weather Rev. 123 (1995) 489–514. [10] A. Robert, A stable numerical integration scheme for the primitive meteorological equations, Atmos.-Ocean 19 (1981) 35–46. [11] A. Robert, A semi-lagrangian and semi-implicit numerical-integration scheme for the primitive meteorological equations, J. Meteor. Soc. Jpn. 60 (1982) 319–325. [12] A.P. Siebesma, J. Teixeira, An advection-diffusion scheme for the convective boundary layer: description and 1D results, in: Proc. AMS 14th Symposium on Boundary Layers and Turbulence, Aspen, Colorado, 2000. [13] P.K. Smolarkiewicz, J.A. Pudykiewicz, A class of semi-Lagrangian approximations for fluids, J. Atmos. Sci. 49 (1992) 2082–2096. [14] P.M.M. Soares, P.M.A. Miranda, A.P. Siebesma, J. Teixeira, An eddy-diffusivity/mass-flux scheme for dry and shallow convection, Q. J. Roy. Meteor. Soc. 130 (2004) 3365–3384. [15] J. Staniforth, C ˆot´e, Semi-Lagrangian integration schemes for atmospheric models–ar e v i e w ,M o n .Weather Rev. 119 (1991) 2206–2223. [16] J. Teixeira, Stable schemes for partial differential equations: the one-dimensional diffusion equation, J. Comput. Phys. 153 (1999) 403–417. [17] J. Teixeira, Boundary Layer Clouds in Large Scale Atmospheric Models: Cloud Schemes and Numerical Aspects, ECMWF, Reading, United Kingdom, 2000. [18] J. Teixeira, Stable schemes for partial differential equations: the one-dimensional reaction–diffusion equation, Math. Comput. Simulation 64 (2004) 507–520. [19] J. Teixeira, J.P. Ferreira, P.M.A. Miranda, T. Haack, J. Doyle, A.P. Siebesma, R. Salgado, A new mixing-length formulation for the parameter- ization of dry convection: implementation and evaluation in a mesoscale model, Mon. Weather. Rev. 132 (2004) 2698–2704. [20] M. van Loon, Tests on semi-Lagrangian transport and interpolation, Report NM-R9217, Department of Numerical Mathematics, Amsterdam, Netherlands, 1992. [21] D.L. Williamson, Application of semi-lagrangian methods to climate simulation, in: Proc. Semi-Lagrangian Methods, ECMWF, Reading, United Kingdom, 1996, pp. 167–181. [22] D.L. Williamson, J.G. Olson, A comparison of semi-Lagrangian and Eulerian polar climate simulations, Mon. Weather. Rev. 126 (1998) 991–1000. [23] D.L. Williamson, J.G. Olson, B.B. Boville, A comparison of semi-Lagrangian and Eulerian tropical climate simulations, Mon. Weather. Rev. 126 (1998) 1001–1012.