Articles | Volume 16, issue 24
https://doi.org/10.5194/gmd-16-7509-2023
https://doi.org/10.5194/gmd-16-7509-2023
Development and technical paper
 | 
22 Dec 2023
Development and technical paper |  | 22 Dec 2023

An improved version of the piecewise parabolic method advection scheme: description and performance assessment in a bidimensional test case with stiff chemistry in toyCTM v1.0.1

Sylvain Mailler, Romain Pennel, Laurent Menut, and Arineh Cholakian
Abstract

This study presents a novel method to estimate the performance of advection schemes in numerical experiments along with a semi-realistic, non-linear, stiff chemical system. This method is based on the examination of the “signature function”, an invariant of the advection equation. Apart from exposing this concept in a particular numerical test case, we show that a new numerical scheme based on a combination of the piecewise parabolic method (PPM) with the flux adjustments of Walcek outperforms both the PPM and the Walcek schemes for inert tracer advection as well as for advection of chemically active species. From a fundamental point of view, we think that our evaluation method, based on the invariance of the signature function under the effect of advection, offers a new way to evaluate objectively the performance of advection schemes in the presence of active chemistry. More immediately, we show that the new PPM + W (“piecewise parabolic method + Walcek”) advection scheme offers chemistry-transport modellers an alternative, high-performance scheme designed for Cartesian-grid Eulerian chemistry-transport models, with improved performance over the classical PPM scheme. The computational cost of PPM + W is not higher than that of PPM. With improved accuracy and controlled computational cost, this new scheme may find applications in other fields such as ocean models or atmospheric circulation models.

1 Introduction

Chemistry-transport models are models that aim at representing the concentration of trace gases and particles in the atmosphere. Many such tools exist and are used for several purposes, including research and operational forecast. The core of such models consists of a chemical solver adapted to stiff ordinary differential equation (ODE) systems along with a framework for solving the advection equation for all the chemical species.

Among the possible strategies to solve the advection equation in chemistry-transport models are the flux-based advection schemes, based on the ideas of Godunov (1959), including the Van Leer (1977) and Colella and Woodward (1984) schemes. The Walcek (2000) scheme, an improvement of Van Leer (1977), has also been of common use in chemistry-transport models. For example, GEOS-Chem provides an advection framework based on the FV3 module implementing the Putman and Lin (2007) method (Martin et al.2022), based on the Colella and Woodward (1984) scheme for 1D advection. The Colella and Woodward (1984) PPM scheme is also implemented in the CMAQ model (Byun and Schere2006; Zhao et al.2020). The CHIMERE model (Menut et al.2021) also provides the Van Leer (1977) and Colella and Woodward (1984) schemes for horizontal advection, while vertical advection can be treated with either the Van Leer (1977) scheme or the Després and Lagoutière (1999) antidiffusive advection scheme (Lachatre et al.2020). The Walcek (2000) advection scheme, an improved version of Van Leer (1977) with reduced numerical diffusion, is also used in chemistry-transport models including CCATT-BRAMS (Freitas et al.2012) and LOTOS-EUROS (Timmermans et al.2022). Therefore, while not the only ones used in chemistry-transport models, the schemes we study here are the base of the numerical resolution of advection in some of the most common chemistry-transport models. Other popular schemes include versions of the Bott (1989) scheme, which is less diffusive than PPM but has the inconvenience of being non-monotonic and therefore tends to generate extreme values or oscillation in the presence of large concentration gradients (Byun and Schere2006), while, by construction, PPM, Van Leer and Walcek schemes enforce mass conservation and monotonicity (Van Leer1977; Colella and Woodward1984; Walcek2000). Despite the availability of higher-order and less diffusive schemes such as the Prather scheme (Prather1986) or MPDATA (e.g. Waruszewski et al.2018), PPM is still considered “highly accurate and efficient enough to be useful” (Harris et al.2021). A new implementation of the PPM scheme in the CAMx chemistry-transport model has been designed recently to take advantage of GPU computational facilities (Cao et al.2023). Since the PPM scheme and, to a lesser extent, the Van Leer (1977) and Walcek (2000) schemes are widely used at least in chemistry-transport modelling, it is important to look for ways to improve these schemes while maintaining their desirable properties of robustness and numerical efficiency. In this direction, the goal of the present study is to compare and assess the performance of these schemes in a bidimensional, academic framework including active chemistry; to build an improved version of the PPM scheme, the PPM + W scheme; and to compare the performance of this new scheme to the above-cited classical schemes.

In the past, many studies have focused on developing, improving and evaluating advection schemes (e.g. LeVeque1996; Nair and Lauritzen2010; Lauritzen et al.2012, 2014), but very few studies have tackled the evaluation of numerical systems combining advection and chemistry in an academic framework. The most significant step in this direction is the study of Lauritzen et al. (2015), who introduced a toy chemistry scheme mimicking the photolysis and recombination of a virtual stratospheric species. With this simplified non-linear chemical system, they have tested chemistry-advection combinations, and the errors they generate as diagnostics of mass-conservation issues. They have noted that combining such a chemical system with the advection solver may reveal problems that are not generated by inert tracer advection. More recently, Lachatre et al. (2022) have shown that changes in the advection formulation may have significant effects on the behaviour of non-linear chemical processes in the troposphere (in their case, the oxidation pathways of SO2 in a mid-tropospheric volcanic plume). Therefore, like Lauritzen et al. (2015), we feel that it is important to test advection schemes not only with inert tracers but also with active chemistry.

Our goal is to provide such a test case for conditions more representative of tropospheric chemistry at the scale of an urban area and to deploy new tools to evaluate advection schemes in the presence of active chemistry. To meet this objective, apart from classical methods and metrics, we introduce a novel idea, the “signature function”, that permits giving a lower bound of model error compared to the exact solution for problems with inert tracer advection and isolating the error due to advection itself in problems including active chemistry. Even though this method is related to the area coordinate introduced by Nakamura (1996), we introduce a new formulation of this idea along with a way to use it to construct a new error estimate which can be used in problems of pure advection as well as in advection with active chemistry. Apart from this novel way of evaluating advection error in cases with active chemistry, we also propose a new “hybrid” advection scheme, the piecewise parabolic method + Walcek (PPM + W) scheme, made of the PPM scheme with the Walcek (2000) flux adjustments in the vicinity of the extrema.

The PPM + W advection scheme and the concept of the signature function are tested within the toyCTM academic chemistry-transport model, already used in Mailler et al. (2021) to test the use of the antidiffusive scheme of Després and Lagoutière (1999) for vertical advection in the atmosphere. The model version used for this study, toyCTM v1.0 (Mailler and Pennel2023), includes horizontal advection with the following schemes (according to user's choice): Godunov (1959), Van Leer (1977), Walcek (2000), Colella and Woodward (1984), and PPM + W (present study); the Després and Lagoutière (1999) scheme is also available for the vertical direction. Chemical processes are solved using an Euler backward iterative (EBI) method (Hertel et al.1993). As reviewed in Cariolle et al. (2017), this EBI scheme or closely related schemes are used in the MOZART model (Emmons et al.2010), the ECHAM5-HAMMOZ model (Pozzoli et al.2008), the TM5 model (Huijnen et al.2010) and the UKCA climate-composition model (O'Connor et al.2014; Esentürk et al.2018).

In Sect. 2, we describe the flux and chemistry of the numerical experiment we have implemented. In Sect. 3, we present the set of simulations we have performed and analysed as well as a description of the advection schemes (Sect. 3.1), the chemical solver (Sect. 3.2) and the time-stepping strategy (Sect. 3.3) implemented in toyCTM v1.0. In Sect. 4, we present the concept of signature function that we introduce in this study for the analysis of simulation results. Section 5 compares and discusses the results obtained with the various numerical schemes, and our conclusions are presented in Sect. 6.

2 Numerical-experiment description

2.1 Chemical mechanism

The chemical mechanism used here (Reactions R1R12) includes a subset of the main reactions of tropospheric gas-phase chemistry. Reactions (R1)–(R3) are the three reactions that constitute the Leighton system, and Reactions (R4) to (R7) account for the formation of the hydroxyl radical OH through the photolysis of ozone in the presence of water vapour. Reactions (R8) and (R9) account for the production of the hydroperoxyl radical through oxidation of CO and the oxidation of NO into NO2 by HO2. Reactions (R10)–(R11) are “termination reactions” that consume the radical species, and Reaction (R12) describes the final consumption of the NOx species by formation of nitric acid.

(R1)NO2+hνNO+O(R2)O+O2+MO3+M(R3)NO+O3NO2+O2(R4)O3+hνO(1D)+O2(R5)O(1D)+H2OOH+OH(R6)O(1D)+N2O+N2(R7)O(1D)+O2O+O2(R8)CO+OHCO2+HO2(R9)NO+HO2NO2+OH(R10)HO2+HO2H2O2+O2(R11)OH+HO2H2O+O2(R12)NO2+OHHNO3

The reaction constants of Reactions (R1)–(R12) have been taken mostly from Seinfeld and Pandis (1997), with a temperature of 298 K and pressure of 101 325 Pa. The photolysis rates have been set to typical midday values (e.g. Mailler et al.2016):

  • jR1=8×10-3s-1

  • jR4=2.5×10-5s-1.

Apart from the chemically active species defined in Reactions (R1)–(R12), we introduce an inert tracer denoted TRC, which undergoes no chemical reaction and is passively advected by the flow. In chemistry-transport models, Reactions (R4)–(R7) are typically lumped in one single reaction, O3+hνOH+OH, with a pseudo-reaction rate constant that depends on the concentration of air molecules and of water vapour molecules and on applying the quasi-steady-state approximation to O(1D). For this study, we have chosen to treat O(1D) as a prognostic species to preserve the full chemical stiffness of the problem, with lifetimes ranging from 4×10-9 s for O(1D) to several days for CO.

Of course key processes like oxidation of methane and of other volatile organic compounds are not taken into account in the above mechanism, but it retains some key features of tropospheric chemistry, which we think important:

  • extreme stiffness;

  • OH production, which depends on the presence of ozone, water vapour and sunlight;

  • non-linear behaviour of ozone production (in this simplified system, ozone production depends on the simultaneous presence of nitrogen oxides, OH and available CO for oxidation).

2.2 Definition of test case

2.2.1 Simulation domain

The simulations are performed on a domain D=0,L×0,L×0,H, where L=105 m and H=1000 m. Since we will use only barotropic winds, the problem is in fact bidimensional in xy, with no z dependance. However, the choice has been made to formally treat the problem as tridimensional in order to be able to use quantities such as air density and reaction rate constants with their usual magnitudes and units. Due to the barotropic nature of the problem, discretization in the vertical direction is in one single cell, while the x and y dimensions are split evenly into n=25 subintervals each. This corresponds to a resolution δx=4×103 m, rather typical of regional-scale chemistry-transport modelling. Domain 𝒟 is therefore discretized into n2 cells, each cell with thickness H and horizontal section δx2=L2n2.

2.2.2 Wind field

The flow we use in this study is the swirling deformational flow introduced by LeVeque (1996) (their Eqs. 9.5–9.6):

(1) u = L T sin 2 π x L sin 2 π y g t , v = - L T sin 2 π y L sin 2 π x g t ,

with

(2) g ( t ) = cos π t T .

T is the half period of the experiment, and the design of the flow is such that all fluid particles are back at their original location after time T, but in between they have undergone a deformation, which is maximal at time T2. Here, while the LeVeque (1996) study formulates the problem with non-dimensional scales for time and space, we set a dimensional scale length L=105 m and half period T=86 400 s. The velocity field corresponding to these values is depicted in Fig. 1. Equation (1) ensures that the wind is zero at domain boundaries (x0,L or y0,L) so that no mass enters or leaves domain 𝒟. Therefore, no boundary conditions for concentrations are needed.

The time-dependent streamfunction for this flow is

(3) ψ x , y , t = - L 2 π T sin 2 π x L sin 2 π y L g t .
https://gmd.copernicus.org/articles/16/7509/2023/gmd-16-7509-2023-f01

Figure 1Streamlines (black contours), wind vectors and wind module in m s−1 (colour shades) at t=0 for the swirling deformational wind field defined in Eq. (1) with L=105 m and T=86 400 s.

Download

2.3 Initial conditions

The numerical experiments will be conducted in the domain 𝒟 defined above, with the chemical scheme described above. To define the initial conditions, we introduce a concentration profile (between 0 and 1) as follows:

(4)φx,y=sin22πxLsin22πyLifx<L2andx<L2,(5)φx,y=0otherwise.

The initial conditions are defined as follows (in terms of the mixing ratio):

  • αTRC=100ppb×φx,y (see Fig. 2)

  • αO3=30 ppb

  • αCO=500 ppb

  • αNO=100ppb×φx,y

  • αNO2=10ppb×φx,y

  • αH2O=8.044×10-3.

For water vapour, the mixing ratio αH2O=8.044×10-3 corresponds to a specific humidity of 5 g kg−1. The initial concentrations of the other active species (O, O(1D), OH, HO2, CO2, H2O2, HNO3) are initialized to zero. They will be produced by Reactions (R1), (R4), (R5), (R8) and (R12).

https://gmd.copernicus.org/articles/16/7509/2023/gmd-16-7509-2023-f02

Figure 2Initial mixing ratio of TRC, proportional to φ(x,y) (defined in Eqs. 45), discretized on domain 𝒟 with n=25 subintervals (δx=δy=4000 m). Note that, at t=0, αNO=αTRC and αNO2=αTRC10 have the same spatial distribution.

Download

Finally, another species of inert tracer, TRCb, is introduced so that

(6) α TRCb + α NO + α NO 2 = 110 ppb .

This species is designed so that, at the initial time and later along the run, in the exact solution (but not necessarily in numerical solutions), all along the run, the sum αTRCb+αNO+αNO2+αHNO3 is uniform, constant and equal to 110 ppb.

3 Numerical methods

3.1 Advection schemes

3.1.1 Existing advection schemes

The following existing advection schemes have been tested in the study:

  1. Godunov (1959)

  2. Van Leer (1977)

  3. Walcek (2000)

  4. PPM (Colella and Woodward1984).

These schemes are flux-based, upwind-biased, semi-Lagrangian schemes based on polynomial reconstructions of the average concentrations. These polynomial reconstructions are piecewise-constant for Godunov (1959), piecewise-linear for Van Leer (1977) and Walcek (2000) (Fig. 3a, b), and piecewise parabolic except in the vicinity of the extrema for PPM (Fig. 3c). The Van Leer (1977) scheme exhibits a discontinuity in the vicinity of the maximum, with the concentrations having a positive jump towards the maximum (Fig. 3a). As a consequence of this discontinuity, due to the upwind-biased strategy, the fluxes going out of the maximum (from the high side of the discontinuity) will be systematically overestimated compared to the fluxes going into the maximum (from the low side of the discontinuity), thereby tending to advect too much mass out of the maximum and not enough mass into the maximum. Walcek (2000) presents his scheme as a way to counteract this bias, by adjusting the flux estimates in the cells next to the maximum, in order to intentionally overestimate the fluxes going into the maximum to counteract the excessive estimation of the fluxes out of the maxima (Fig. 3b).

3.1.2 The PPM + W scheme

The PPM scheme presents the same caveat as Van Leer in the vicinity of extrema, with a strong discontinuity on each side of the extremum (Fig. 4a), with the effect of underestimating the mass flux into the maximum and overestimating the mass flux out of the maximum. Therefore, since Walcek (2000) has proven that his flux adjustments in the vicinity of the extrema are successful in improving the Van Leer scheme, it makes sense to try applying the same flux adjustments to PPM, which seems to have a behaviour similar to Van Leer (1977) in the vicinity of mixing-ratio extrema (Fig. 3a, c). To test this idea, we design a new scheme based on PPM but applying the Walcek flux adjustments in the vicinity of the extrema (Figs. 3d and 4b). We call this scheme PPM + W, standing for the “piecewise parabolic method + Walcek flux adjustments”. The PPM + W has the same behaviour as Walcek in the extrema and the neighbouring cells and the same behaviour as PPM in all other cells.

We detail here the procedure applied for this scheme. Let (αi) be the values of the mixing ratio in the model cells numbered by the 1D index i, with a Courant number ν. For the sake of simplicity we assume that δx=1. The objective of this procedure is to calculate the average mixing ratio between xi and xi+12-ν, α̃i+12 (the rest of the implementation of advection from this estimate is detailed in Lachatre et al.2020).

The procedure is as follows.

  • If αi+1-αiαi+2-αi+1>0 and αi-1-αi-2αi-αi-1>0,
    the current cell is not a neighbour of a maximum. We estimate α̃i+12 following the piecewise parabolic method procedure described in Colella and Woodward (1984).

  • Otherwise we estimate the Walcek-adjusted flux as follows:
    s=signαi+1-αimin12αi+1-αi-1,2αi+1-αi,2αi-αi-1, the Van Leer slope;
    if αi+1-αiαi+2-αi+10, β=1.75-0.45ν,
    else β=max1.5,1.2+0.6ν,
    α̃i+12=αi+121-ν×β×s.

The β coefficient (β>1 by construction) is introduced by Walcek (2000) to steepen the Van Leer slopes in the vicinity of the maxima to obtain the desired overestimation of tracer fluxes into the maxima. These steepened slopes are visible in Fig. 3b for the Walcek (2000) scheme and in Figs. 3d and 4b for the PPM + W scheme.

https://gmd.copernicus.org/articles/16/7509/2023/gmd-16-7509-2023-f03

Figure 3Reconstruction of a Gaussian mixing-ratio profile by (a) the Van Leer (1977) scheme, (b) the Walcek (2000) scheme, (c) the PPM scheme and (d) the PPM + W scheme. The x axis is a non-dimensional space coordinate. The reconstruction has been performed for a Courant number ν=0.4. The reconstructed fields are presented with alternating red and blue colours to enhance the discontinuities between neighbouring cells.

Download

https://gmd.copernicus.org/articles/16/7509/2023/gmd-16-7509-2023-f04

Figure 4(a) Same as Fig. 3c but zoomed in to within the vicinity of the mixing-ratio maximum and (b) same as Fig. 3d but zoomed in to within the vicinity of the mixing-ratio maximum.

Download

3.1.3 Computational cost of the advection schemes

To evaluate the computational cost of these advection schemes, advection of a 1D vector composed of 2×105 cells has been performed over 520 time steps, corresponding to 1.04×108 calls to the reconstruction routine, plus the update of the mixing-ratio values at each time step. The calculation time for all these advection schemes is presented in Table 1, showing that the schemes using linear reconstruction (Van Leer and Walcek) are less costly than the schemes using parabolic reconstruction (PPM and PPM + W) due to the simpler calculation. Interestingly, the computation cost of PPM + W is slightly smaller than the cost of PPM, possibly because in the cells neighbouring a concentration extremum the reconstruction is linear in PPM + W instead of parabolic in PPM (Fig. 4).

Table 1Mean calculation time per cell and per time step for the five advection schemes retained for the present study. The calculation has been performed in Fortran, a programming language frequently used for operational chemistry-transport models, on a laptop with an Intel Core i7-1165G7 CPU.

Download Print Version | Download XLSX

3.1.4 Convergence properties

A first comparison between the PPM + W advection scheme and the other four tested schemes has been performed in terms of numerical convergence. For this purpose, an inert tracer with a squared cosine bell distribution (initially) has been advected over a 1D periodic domain of unit length, with a constant and uniform speed, with a unit duration. This convergence test has been performed dividing the domain into 10, 20, 40, 80, 160 and 320 grid cells and with a Courant number of 0.5. The results of this convergence test are shown in Fig. 5. Several features can be observed in Fig. 5. The convergence rates, defined as the opposite of the log–log slope between the two last data points (nx=160 and nx=320), are given in Table 2. From these results, several observations can be made.

First, the PPM + W advection scheme performs better on this simple convergence test than the other tested schemes. In particular, throughout all the resolution range, the error obtained in both ∥⋅∥1 (Fig. 5a) and ∥⋅∥2 (Fig. 5b) with the PPM + W scheme is 30 % to 50 % lower than with the classical PPM scheme. This difference persists even for high resolutions. On the contrary, the Walcek scheme strongly outperforms the Van Leer scheme for coarse resolutions, but this difference tends to diminish for higher resolution. In other terms, it looks like the Walcek flux corrections when applied to the Van Leer scheme permit the improvement of accuracy when model resolution is coarse but that the same flux corrections when applied to the PPM scheme improve accuracy even for fine resolutions. Finally, the convergence rate for PPM + W is similar to that of PPM, in our case around 2.5 for the ∥⋅∥1 convergence rate and around 2 for ∥⋅∥2 (Fig. 2).

https://gmd.copernicus.org/articles/16/7509/2023/gmd-16-7509-2023-f05

Figure 5(a) The ∥⋅∥1 error as a function of the number of nx for the five advection schemes used in the present study; (b) same as (a) for the ∥⋅∥2 error.

Download

Table 2Convergence rate of the five advection schemes used in the present study for ∥⋅∥1 and ∥⋅∥2.

Download Print Version | Download XLSX

3.2 Chemistry solver: the Euler backward iterative method

The stiff chemical system is integrated using an EBI method. As described in Hertel et al. (1993) and Cariolle et al. (2017), we obtain the concentration vector c(t+δtchem) as the solution of

(7) c t + δ t chem = c t + δ t chem P t + δ t chem 1 + δ t chem L t + δ t chem .

For the present study, the focus is to test the performance of the advection scheme in articulation with active chemistry. Due to this focus, we limit errors in the resolution of the chemical system using a short time step for chemistry (δtchem=20 s).

Equation (7) is a non-linear, fixed-point equation and can be solved only numerically, usually with an iterative method. Formally, Eq. (7) guarantees exact mass conservation. However, this is true only if a very good convergence of the solution is reached (Cariolle et al.2017). To limit violation of mass conservation in our study, we have set a very strict convergence criterion for the iterative resolution of Eq. (7), stopping iteration when the estimate of c(t+δtchem) yields a relative difference less than ε for each species between c(t+δtchem) and ct+δtchemPt+δtchem1+δtchemLt+δtchem. The convergence parameter is set to ε=10-4 in the UKCA chemistry-transport model (Esentürk et al.2018) and ε=10-6 in the present study. This very strict convergence criterion is in line with the short chemical time step to obtain the best-possible numerical solution of the chemical evolution of the system, even at the cost of slow computations. Again, this choice is due to the purpose of this study to test the performance of the advection scheme, limiting as much as possible the errors in the chemical solver.

3.3 Time stepping

The advection time step δtadv has been set to δtadv=1800 s. With maximal wind module U≃1.8 m s−1 (Fig. 1) and δx=4×103 m, this yields a maximal Courant number νmax≃0.8. The time-stepping strategy follows Strang-style time stepping (Strang1968), with the steps as follows:

  1. Integrate chemistry over δtadv2 (45 chemical time steps).

  2. Integrate zonal advection over δtadv2.

  3. Integrate meridional advection over δtadv.

  4. Integrate zonal advection over δtadv2.

  5. Integrate chemistry over δtadv2 (45 chemical time steps).

Table 3 summarizes the six simulations that have been performed. The above-described case with Reactions (R1)–(R12) has initial conditions as described in Sect. 2.3, a domain shape and discretization as described in Sect. 2.2.1, advection schemes as described in Sect. 3.1, and a chemical solver as described in Sect. 3.2. Along with the simulations performed with each of the five advection schemes, a “base” simulation has been performed with the same setup as the other simulations but without advection. As discussed later (Sect. 4), this base simulation will serve as a benchmark to estimate advection errors in the other five simulations.

Godunov (1959)Van Leer (1977)Walcek (2000)Colella and Woodward (1984)

Table 3Summary of the main characteristics of the simulations that have been performed.

Download Print Version | Download XLSX

3.4 Conservation properties

It is worth noting that, by construction, flux-based advection integration is mass-conservative since the mass flux out of a cell through a facet is compensated exactly by the mass flux into the neighbouring cell through the same facet. Equation (7) also guarantees mass conservation in the Euler backward scheme as soon as the chemical reactions themselves are balanced (Hertel et al.1993) (which is the case of Reactions R1R12 except for the imbalance of dioxygen in Reaction R8 due to integrating reaction H+O2HO2 into the kinetically limiting time step CO+OHCO2+H).

However, due to the finite number of iterations in the iterative resolution of Eq. (7), mass conservation is only enforced with a finite precision of ε=10-6 (see Sect. 3.2). Therefore, the relative mass imbalance in the outputs for C, active N (without taking into account N2) and H can be expected to be of the order of ε, since the Euler backward iterative scheme (Eq. 7) is in principle mass-conservative, but mass conservation is obtained “only if a good convergence of the solution is reached” (Cariolle et al.2017).

Mass calculations for C, active N, H and TRC have been performed between the beginning and the end of the simulations. The results of this calculation for the PPM + W and base simulations are given in Table 4, showing that the relative mass imbalance at the end of simulation is 10-6 for active N, 10-12 for C and H, and 10-15 for TRC (which has no chemistry and is therefore affected only by advection which, as discussed above, ensures mass conservation up to numerical accuracy for the flux-based schemes we have implemented).

The imbalance results are similar for all simulations, including the base simulation that has no advection, which shows that the small mass imbalance for chemically active species (up to 10-6) is essentially due to the finite precision in the Euler backward iterative chemistry solver.

Table 4Mass-conservation diagnostic for C, active N, H and TRC in the PPM + W and base simulations.

Download Print Version | Download XLSX

In summary, the integration strategy we introduce above permits the conservation of mass (up to numerical accuracy for inert species and up to an arbitrary numerical tolerance defined by the user, in our case ε=10-6, for chemically active species) and conserves initially uniform mixing ratios up to numerical accuracy for the species with no active chemistry. All the advection schemes presented above are also built to respect monotonicity: they do not create new mixing-ratio extrema.

4 The signature function

4.1 Accuracy of inert tracer concentrations

As the LeVeque (1996) flow is designed so that at t=T every Lagrangian particle is back at its original location, it is possible to estimate the accuracy of numerical simulation by comparing the final simulated tracer concentration field to its initial value, therefore giving access to the magnitude of numerical error.

In the present study, we will estimate model error in ∥⋅∥1, introducing E1 as the normalized 1 error in the mixing ratio:

(8) E 1 = i = 1 N ρ i V i α i t - α i 0 i = 1 N ρ i V i α i 0 ,

where index i spans the entire domain. In the present case where density ρi and cell volume do not vary across cells, Eq. (8) boils down to

(9) E 1 = i = 1 N α i t - α i 0 i = 1 N α i 0 .

4.2 The signature function for inert tracer advection

Here we introduce a new idea to evaluate advection schemes. As far as we know, this idea has not been tested in the past literature but resembles the area-coordinate formulation used by Nakamura (1996).

Let us imagine a fluid with density ρx;y;z;t in a three-dimensional domain 𝒟, advecting a tracer having initially a mixing ratio αx;y;z;t=00,1, following equations

(10)ρt+ρu=0,(11)αt+uα=0.

For any given time t and any mixing ratio 0X1, we can define 𝒮t(X) as the mass of fluid in the volume 𝒟t(X) defined as the set of all x;y;z where the tracer mixing ratio αx;y;z;t<X, divided by the entire mass of fluid in 𝒟:

(12) S t X = D H X - α x ; y ; z ; t ρ d V D ρ d V ,

where H is the Heaviside step function (H(u)=1 if u>0; H(u)=0 if u≤0). The 𝒮t function can be, in some sense, interpreted as a mass-weighted cumulative probability density function of the tracer mixing ratio. If we reduce this definition to 2D flows with uniform density, 𝒮t is related to the reciprocal function of the area-coordinate formulation of Nakamura (1996).

With this definition, we always have St(0)=0 (for all t values) and St(X)→1 when X+ (more precisely, St(X)=1 as soon as X is larger than the maximum value of αx;y;z;t=0 over domain 𝒟).

Equation (12) makes clear that function 𝒮t is invariant during the motion: for any given value of X, 𝒮0(X) is the (normalized) mass of the fluid parcel 𝒟0(X) that has a tracer mixing ratio α<X at t=0. We can observe that, since mixing ratio in Lagrangian parcels is preserved by pure advection (Eq. 11), 𝒟t(X) and 𝒟0(X) represent the same Lagrangian fluid parcel at a different time. We also know that the mass of fluid in Lagrangian parcels is constant in time due to mass conservation for the carrier fluid (Eq. 10). Since the total mass of fluid in 𝒟 is also constant in time for the same reason, this implies that, for all t and X values, 𝒮t(X)=𝒮0(X). In other words, the signature function 𝒮t is an invariant in time of the advection equation. The invariance of this function holds for both divergent and non-divergent flows.

Therefore, since we know that, for the exact solution, 𝒮t=𝒮0, the departure of the numerical evaluation of function 𝒮t from 𝒮0 in numerical simulations can be used as an objective measurement of discretization errors.

In practice, in an Eulerian model discretized into N cells, each cell has an evaluation of the tracer mixing ratio αit. The evaluation of 𝒮t is straightforward:

(13) S t X i = 1 N H X - α i ρ i V i i = 1 N ρ i V i .

Generally, the numerical evaluation of 𝒮t in an Eulerian model will differ from 𝒮0: the signature of the initial tracer distribution will evolve under the effect of the errors of the advection scheme, and the magnitude of the signature modification can serve as a measure of the extent of advection error. In the particular case in which the flow is non-divergent and the carrier fluid mass ρiVi is the same in all model cells, the norm-1 difference between 𝒮t and 𝒮0 can be calculated as

(14) 0 S t X - S 0 X d X = 1 N i = 1 N α t ̃ i - α 0 ̃ i ,

where αt̃i is the vector of all mixing ratios in model cells sorted in increasing order. This simplified formula holds only in the specific case in which the flow is invariant and the mass of fluid is distributed evenly between all the model cells. In this case, and only in this case, it is also convenient to introduce the normalized norm-1 difference 1 between 𝒮t and 𝒮0 as

(15) E 1 = i = 1 N α t ̃ i - α 0 ̃ i i = 1 N α 0 ̃ i .

Figure 6 shows an example of how comparing the 𝒮t signature function with 𝒮0 permits comparison of the accuracy of two simulations at a time when no analytic solution is easily accessible. Figure 6a–b show the mixing ratio for TRC at T2 in the Godunov and Van Leer simulations, respectively. Without access to the exact solution, it is hard to compare quantitatively the quality of these simulations at that stage, even though indicators such as the maximal tracer mixing ratio can give partial information. Figure 6c–d show the ST2 function compared with 𝒮0 for the Godunov and Van Leer simulations, respectively. This graphical representation permits giving an intuitive meaning to the normalized signature error 1 as the total area between the representative curves of 𝒮t and 𝒮0 (shaded in Fig. 6c–d), divided by the total area left of the representative curve of 𝒮0 (hatched in Fig. 6c–d). This measure gives an indication of model error based on not only one particular point representative of a part of the tracer distribution (e.g. the maximal value, a very partial indicator), but also an integrated error indicator taking into account the maximal and minimal values as well as the entire tracer distribution in between these values. In the case presented in Fig. 6c–d, the area between 𝒮0 and 𝒮t is smaller in the Van Leer simulation (Fig. 6d, with 1=0.218) than in the Godunov simulation (Fig. 6c, with 1=0.571).

https://gmd.copernicus.org/articles/16/7509/2023/gmd-16-7509-2023-f06

Figure 6(a) The TRC mixing ratio as simulated in the Godunov simulation at time T2; (b) same as (a) for the Van Leer simulation. (c) The signature function for the TRC mixing ratio in the Godunov simulation compared to the base simulation (𝒮0, green line); (d) same as (c) for the Van Leer simulation. The signature error 1 is equal to the ratio of the shaded area between the representative curves of 𝒮0 and 𝒮t to the hatched area left of 𝒮0.

Download

4.3 The signature function for advection of active species

The concentrations of active species evolve not only under the effect of advection but also due to chemical reactions. Therefore, the time invariance of the signature function does not hold for these species. However, the signature function can still be used to compare and evaluate simulations if one remarks that, in the case we study here with no variations in air density and air temperature, and with no emissions, the chemistry that takes place in each Lagrangian air parcel is independent of its position. Therefore, for all species, the signature function at time t, 𝒮t, should theoretically be the same as in the base simulation with no advection. If we note St, the signature function at time t in the simulation without advection, for all t values we have St=St: at any time, the signature of the distribution of every chemical species should be the same in the case with advection as in the case without. In other words, advection should only deform the map of all chemical species and not change the chemistry within each Lagrangian air parcel. With a non-linear chemical system such as Reactions (R1)–(R12), there is no easy access to the exact solution of the system even without advection: St is not known exactly for the chemically active species. However, the base simulation resolves the chemical reactions with exactly the same chemical solver as the other simulations but without motion. Therefore, comparing the distribution of a chemical species at time t in a simulation with advection to its distribution in the base simulation without advection will permit having an estimate of how the numerical errors in advection affect the distribution of chemically active species.

As an illustration of this, Fig. 7 shows the ozone mixing ratio at T2, as simulated in three simulations with advection (Fig. 7a–c) and the base simulation without advection (Fig. 7d). This figure alone does not make it easy to discriminate between the three numerical simulations. Figure 8 shows the signature function for the same three simulations (Van Leer, PPM and PPM + W): this time, both the visualization of the agreement between 𝒮t and St and the objective calculation of 1 clearly show that the PPM + W simulation agrees better with the base simulation in terms of the signature function (1=0.0622), followed by the PPM simulation (1=0.0679) and the Van Leer simulation (1=0.0787). Apart from this quantitative agreement, some qualitative and local conclusions can be drawn from the graphical comparisons between 𝒮t and St, such as the fact that the representation of the ozone minimum in the PPM + W simulation is more accurate than in the PPM simulation, which is visible in the initial part of the signature function (αO3<30 ppb) where the shaded area is smaller in PPM + W (Fig. 7c) than in PPM and Van Leer.

https://gmd.copernicus.org/articles/16/7509/2023/gmd-16-7509-2023-f07

Figure 7Ozone mixing ratio at T2 as simulated in (a) the Van Leer simulation, (b) the PPM simulation, (c) the PPM + W simulation and (d) the base simulation.

Download

https://gmd.copernicus.org/articles/16/7509/2023/gmd-16-7509-2023-f08

Figure 8The signature function for the O3 mixing ratio at T2, ST2, as simulated in (a) the Van Leer simulation; (b) the PPM simulation and (c) the PPM + W simulation. As in Fig. 6, the signature error 1 is equal to the ratio of the shaded area between the representative curves of 𝒮0 and 𝒮t to the hatched area left of 𝒮0.

Download

5 Results and discussion

5.1 Results for classical performance indicators

Figure 9 shows the evolution of some performance indicators along the course of the experiments. For the same reasons as for the signature function, these metrics should be the same in all simulations in the absence of numerical errors in the representation of advection, so the differences between the time series obtained in the simulations with advection and the base simulation reveal the effects of numerical errors in the representation of advection. It is interesting to note that all the metrics represented in Fig. 9 can be derived directly from the signature function.

Regarding ozone extrema throughout the simulation (Fig. 9a), we observe smaller differences between the different simulations. The PPM + W, PPM and Walcek simulations produce very comparable values of ozone maxima. Surprisingly, the simple Godunov schemes represents the ozone maximum slightly better than all other schemes except Walcek towards the end of the simulation. However, the representation of the ozone minimum by the Godunov simulation is very bad, failing to represent a minimum lower than the background value towards the end of the simulation. The PPM + W and Walcek simulations represent the smallest (and closest to the base simulation) minimum value for O3, followed by PPM and Van Leer.

Regarding the preservation of tracer maxima (Fig. 9b), the PPM + W simulation performs best, with a clear edge over the Walcek and PPM simulations (Fig. 9b). For this metric, the Walcek simulation gives better results than the PPM simulation, with a smaller half plume, closer to the base simulation, Walcek and PPM. The PPM and Walcek simulations perform similarly for this metric.

Figure 9c shows evolution of the non-dimensional half volume of the tracer plume during the simulation (defined as the smallest volume containing half of the mass of TRC divided by the total volume of domain 𝒟). As discussed in Lachatre et al. (2020), this is a measure of the diffusivity of the advection schemes. With this metric, we see that the PPM + W simulation performs best, followed by Walcek and PPM.

https://gmd.copernicus.org/articles/16/7509/2023/gmd-16-7509-2023-f09

Figure 9Time series for all the simulations for the (a) minimum and maximum of the O3 mixing ratio, (b) maximum TRC mixing ratio, and (c) half volume of the TRC plume (relative to the volume of the entire domain).

Download

5.2 Accuracy and signature error

Unlike the partial metrics presented in Fig. 9, the normalized 1 signature error in the mixing ratios 1 is a performance diagnostic for the simulations that takes into account the entire distribution of mixing ratios and not just particular values such as the maximum or minimum. Figure 10 shows 1 for TRC, O3, NO and NO2 for all the simulations, with the smallest 1 indicating the best-performing simulation.

The first information we get from Fig. 10 is that, for all these compounds, the PPM + W simulation performs best in this regard. For the case of ozone (Fig. 10b), the 1 time series discriminates much more between the different simulations, with a clear superiority of the PPM + W simulation over the Walcek and PPM simulations, while the differences between these three simulations in the ozone min–max plot (Fig. 9b) appeared small. Analysing the evolution of 1 throughout the simulation shows that, for TRC, NO and NO2, the best-performing simulations are PPM + W, Walcek and PPM, in this order, but that for ozone, PPM performs better than Walcek, a conclusion that could not be drawn from the min–max plot (Fig. 9b), which indicated a better performance of Walcek regarding both the ozone maximum and the ozone minimum.

https://gmd.copernicus.org/articles/16/7509/2023/gmd-16-7509-2023-f10

Figure 10Time series for 1 (1 normalized signature error) for (a) TRC, (b) O3, (c) NO and (d) NO2.

Download

Table 5Normalized 1 error E1 (Eq. 9) and normalized 1 signature error 1 (Eq. 15) at the end of the simulations for O3, NO, NO2 and TRC, compared to the base simulation with no advection. In each column, the lowest error value is in bold font the second-lowest is underlined.

Download Print Version | Download XLSX

While we have shown so far that the analysis of 1 permits the drawing of clearer conclusions regarding the performance of the various simulations, in the present case we can confirm these results by comparing 1 to a more classical metric, the normalized 1 error E1 (Eq. 9) of the simulations with advection relative to the base simulation without advection. Unlike 1, E1 can only be calculated at the final time step (when all the Lagrangian particles are back at their initial location), since the exact solution for t<T is not accessible. Several observations can be made from the values of E1 and 1 in Table 5.

First, in all cases we always have 1<E1. Qualitatively, this can be interpreted as 1 being a weakened form of the 1 error, retaining the differences in the distribution of mixing ratios but eliminating the differences in the location of the tracer plume.

Interestingly, the performance ranking between the five simulations obtained by analysing the signature error 1 is the same for all variables as with E1: for TRC, NO and NO2, PPM + W performs best, followed in this order by Walcek, PPM, Van Leer and Godunov, but for O3, PPM performs better than Walcek. In our study case, analysing 1 permits comparison of the performance of the various simulations without access to the exact solution and gives the same results as the analysis of E1 (which requires access to the exact solution). Even for simulations with comparable performance to the Walcek and PPM simulations, the signature error 1 permits diagnosing which of the two simulations performs better for each variable. Even though the differences between these two simulations are not drastic and depend on the species of interest, the conclusions drawn from the analysis of the signature error 1 are the same and correspond to E1. This gives confidence in the usability of the signature error 1 as a proxy for simulation accuracy when the exact solution is not available.

Having verified this, it is useful to go back to Fig. 10 and interpret the evolution of 1 in time as a hint of when error appears along the simulation. In this regard, we see two very different behaviours between the analysed variables. For TRC, NO and NO2, substantial errors appear almost instantly after 1800 s of simulation (one single time step). This is due to the action of the wind field on the initially very sharp peak of these species (Fig. 2). On the contrary, for O3, having an initially uniform distribution, errors due to advection appear only when sufficient heterogeneity is introduced into the O3 map by chemical processes, since all the advection schemes are built to advect exactly a uniform mixing ratio, maintaining its uniformity. Therefore, the onset of advection errors in O3 is much slower than in the three species that have initially heterogeneous distributions.

5.3 Results at 1 km resolution

To test the impact of higher resolution on our results, we have performed the same test case as above but refining the resolution from 4 to 1 km and accordingly reducing the time step from 1800 to 450 s. The results for this higher-resolution simulation are shown in Table 6. This table shows that, also at this resolution, the PPM + W offers the best performance of all the tested schemes, but, unlike in the 4 km resolution case, this does not hold for all species: at 1 km resolution, the PPM scheme performs comparably to PPM + W for the inert tracer TRC and slightly better for NO2. However, PPM + W clearly performs better for O3 and NO. These results suggest that the improvement brought by PPM + W over PPM may tend to become smaller for higher resolutions (contrary to what we observe in Fig. 5).

Table 6Normalized 1 error E1 and normalized 1 signature error 1 in the 1 km simulations at the end of the simulations for O3, NO, NO2 and TRC, compared to the base simulation at 1 km resolution with no advection. In each column, the lowest error value is in bold font and the second-lowest is underlined.

Download Print Version | Download XLSX

5.4 Conservation of sum species

As shown by Godunov (1959), “among [linear] schemes of second order accuracy [for the advection equation], there is none which satisfies the monotonicity condition”. As a consequence of this result, known as the Godunov theorem, second-order, monotonic schemes for the advection equation are all non-linear. In particular, among all the schemes we study here, all are non-linear due to their specific treatment of the maxima to preserve monotonicity. However, these non-linear effects should be kept as small as possible. To test this property, we have tested the conservation of two sum species, which should stay uniform and constant along the run:

  • The first is αTRCb+αNO+αNO2+αHNO3, where TRCb is an inert tracer whose initial distribution is defined by Eq. (6). The mixing ratio of this sum species should stay constant, uniform and equal to 110 ppb.

  • The second is αCO+αCO2, which should stay constant, uniform and equal to 500 ppb due to the conservation of carbon.

Figure 11 shows that, as imposed by the Godunov theorem, only the Godunov (1959) scheme preserves exactly the uniformity of sum species. Other schemes are higher-order and monotonic, but therefore non-linear, and do not preserve the uniformity of these sum species. The magnitude of the departure of the mixing ratio for these species from its theoretical value is a measure of the non-linearity of advection schemes, which is an undesirable property since the advection equation itself is linear. Figure 11 shows that non-linearities of up to 1 %–2 % appear for all schemes in the representation of CO+CO2 and that they tend to be strongest in the Walcek and the PPM + W schemes. Regarding TRCb+NO+NO2+HNO3, the non-linearities reach up to 10 % of the expected value and are strongest in the PPM + W scheme. These results suggest that non-conservation of sum species in the vicinity of mixing-ratio extrema, which is observed in all second-order monotonic schemes, tends to be stronger when the Walcek (2000) flux corrections are applied, which might represent a drawback of PPM + W compared to PPM. In spite of (or due to) this increasingly non-linear behaviour in the vicinity of the extrema, PPM + W preserves the values of extrema better than PPM (Fig. 9b).

https://gmd.copernicus.org/articles/16/7509/2023/gmd-16-7509-2023-f11

Figure 11Minimum and maximum values in the 4 km simulations (presented in the paper) for αCO+αCO2 (panel a) and αNO+αNO2+αHNO3+αTRCb (panel b).

Download

6 Conclusions

We have introduced the signature function 𝒮t as a sort of mass-weighted cumulative probability distribution function of the tracer mixing ratio (Eq. 12) and shown that 𝒮t is an invariant of the advection equation. This invariant is not a scalar like other classical invariants (tracer mass, minimum and maximum mixing ratio, etc.) but a function and therefore contains much more information than the above-cited, more classical invariants. In fact, these invariants (tracer mass, minimum and maximum mixing ratios, etc.) can be derived directly from 𝒮t so that 𝒮t can be considered a “stronger” invariant. Like other invariants such as the minimum and maximum values of the mixing ratio, it is usually not preserved perfectly by the advection schemes, and the degree of non-conservation of 𝒮t gives a proxy for model error. Based on this idea, we propose the normalized 1 signature error 1 (Eq. 15) as a measure of model error. Graphically, 1 is the normalized area between the simulated curve of 𝒮t and its theoretical curve (e.g. Fig. 6c–d).

In this context, we have shown that the signature function and its normalized 1 error 1 can also be used as an error estimate for chemically active species, even though 𝒮t is not an invariant of the system in this case. This is achieved by comparing a simulation including chemistry and advection to a “companion” simulation with chemistry only but no advection (the base simulation in the present study). For each chemical species and at any time t, these two simulations should in theory have the same 𝒮t, but in practice they do not due to errors in advection. Therefore, errors due to advection can be estimated even for chemically active species by calculating the 1 error between the two functions (Fig. 8).

We have used this new invariant in order to evaluate a new advection scheme that we have designed for this study, based on the PPM scheme (Colella and Woodward1984) with flux corrections in the vicinity of the extrema as in Walcek (2000) (Figs. 34). This new advection scheme, which we propose to name PPM + W (for piecewise parabolic method + Walcek), has been tested for both an inert tracer and chemically active species, with a velocity flux from LeVeque (1996) (Eq. 1). A simplified chemical scheme with 12 reactions has been designed, representing daytime photochemistry of nitrogen oxides in the presence of CO, including short-lived species such as O(1D) or O (Reactions R1R12). With an initial peak of NO and NO2 concentrated over a small area (Fig. 2) and initially uniform concentrations of CO and O3, this case generates a sharp ozone minimum colocalized with the NOx peak, a large area with background O3 concentration and a belt of high ozone concentrations in between (Fig. 7). With this test case, we have evaluated the PPM + W scheme along with the PPM scheme, the Walcek (2000) scheme, the Van Leer (1977) scheme and the Godunov (1959) scheme. In this case, we have shown that, for all species and all the metrics we have tested, the PPM + W scheme performs better than all the other tested schemes (Table 5), with a normalized 1 error E1 of 16.4 % in O3 (19.3 % for PPM), 20.7 % in TRC (24.3 % for Walcek), 28.3 % in NO2 (31.2 % for Walcek) and 37.2 % in NO (42.5 % for Walcek). Table 5 also shows the 1 signature error, which is always smaller than E1. Interestingly, examining the 1 errors for the same variables and the same simulations yields exactly the same conclusions as examining the normalized 1 error E1. This shows that, even without access to an exact solution, the 1 signature error permits comparison of the simulations against each other for inert and reactive species, giving the same conclusions as an accuracy analysis with E1. This being shown, Fig. 10 permits visualization of the evolution of error along the course of the simulation, while E1 can be calculated only at t=T because the LeVeque (1996) flux field is designed to guarantee that at that time all the Lagrangian particles should be back at their initial locations.

Therefore, the conclusion of this study is twofold. First, regarding the signature function as an invariant of the advection equation, we feel that this invariant contains much more information than other invariants that have been typically used to check advection schemes, such as the minimum or maximum values of mixing ratios, while not requiring access to the exact solution. In the case of chemistry-transport models, generalizing this concept to more dimensions (by studying the mass-weighted probability distribution function of all species simultaneously rather than one signature function per species) could be promising. For the same reasons as exposed above, this multidimensional probability distribution function should be an invariant of the advection equation. This concept could be explored to quantify model error in a synthetic way across all variables, instead of separately for each variable. The approach introduced here with the 1 signature error could possibly be generalized using statistical tools such as the Wasserstein distance to compare these multidimensional probability distribution functions with each other. The main limit to the use of the signature function as a tool to evaluate the advection framework in real-world geophysical models is that it relies on mixing-ratio conservation, which holds for pure advection but does not hold in the presence of mixing or diffusion. However, some geophysical compartments like the deep ocean or the stratosphere are substantially affected by mixing only for very long timescales, so the signature function could be a useful tool to verify the behaviour of advection frameworks in such compartments.

From a more applied point of view, the PPM + W advection scheme introduced here performs better than both the Walcek (2000) scheme and the PPM scheme, with all the metrics we have tested and for both inert and active species. It not only preserves the tracer maxima better than the Walcek (2000) scheme (Fig. 9), but also is more accurate than the PPM scheme for the representation of ozone and other active species (Table 5). Of course, these results are proven in the present study only on one 2D test case with active chemistry and on a 1D convergence test. While the results of both the test case and the convergence test consistently indicate a better performance of PPM + W compared to the other advection schemes tested here, they do not include the full range of Courant numbers and tracer patterns that occur in realistic models. For example, an additional numerical experiment presented here (Sect. 5.3 and Table 6) shows that, when the resolution is refined four times compared to our main test case, the advantage of PPM + W over PPM seems to be reduced and, for one species (NO2), PPM + W is even outperformed by PPM in terms of accuracy. However, even in this high-resolution test case, the performance of PPM + W is better than that of PPM for the three other tested species. The convergence test performed in this study (Fig. 5) suggests that there is no systematic reduction in the performance of PPM + W at higher resolution. A possible drawback of the PPM + W scheme when compared to the PPM scheme is its stronger non-linearity in the vicinity of maxima (Fig. 11), which permits better conservation of the maxima themselves but induces more numerical artefacts in the representation of sum species in the vicinity of extrema.

We have also observed (Table 1) that the computation cost of PPM + W is slightly lower than the cost of the PPM scheme, which is used in some of the most popular chemistry-transport models. The improved performance in terms of accuracy and of preservation of tracer extrema without increasing the computational load makes this scheme a very interesting option for chemistry-transport models, in an effort to reduce numerical diffusion, which is important in particular in the presence of non-linear chemistry as discussed in Lachatre et al. (2022).

Code and data availability

This study has been performed using toyCTM v1.0.1 (https://doi.org/10.5281/zenodo.10018706, Mailler and Pennel2023). All the Python scripts used to launch the model and to perform the post-processing of model outputs to obtain the figures in this paper and the numbers in Table 5 are available from https://doi.org/10.5281/zenodo.10018761 (Mailler et al.2023).

The Fortran code AdvBench v1.0.0 used to evaluate advection performance (Table 1) is available from https://doi.org/10.5281/zenodo.7937121 (Mailler2023).

The model toyctm v1.0.1, AdvBench v1.0.0 and all the scripts used to launch the model and post-process its outputs for the present study are distributed under the GNU General Public License v2.0. No datasets were used for this study.

Author contributions

All the authors contributed to the design of the simulated cases; SM performed and analysed the simulations; SM developed the software with RP.

Competing interests

The contact author has declared that none of the authors has any competing interests.

Disclaimer

Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.

Acknowledgements

The numerical calculations in this study benefited from the IPSL Data and Computing Center ESPRI, which is supported by CNRS, SU, CNES and École Polytechnique.

Financial support

This research has been supported by the Agence de l'Environnement et de la Maîtrise de l'Énergie (ESCALAIR).

Review statement

This paper was edited by Sylwester Arabas and reviewed by Christopher Walcek, Hilary Weller, and one anonymous referee.

References

Bott, A.: A Positive Definite Advection Scheme Obtained by Nonlinear Renormalization of the Advective Fluxes, Mon. Weather Rev., 117, 1006–1016, https://doi.org/10.1175/1520-0493(1989)117<1006:APDASO>2.0.CO;2, 1989. a

Byun, D. and Schere, K. L.: Review of the Governing Equations, Computational Algorithms, and Other Components of the Models-3 Community Multiscale Air Quality (CMAQ) Modeling System, Appl. Mech. Rev., 59, 51–77, https://doi.org/10.1115/1.2128636, 2006. a, b

Cao, K., Wu, Q., Wang, L., Wang, N., Cheng, H., Tang, X., Li, D., and Wang, L.: GPU-HADVPPM V1.0: a high-efficiency parallel GPU design of the piecewise parabolic method (PPM) for horizontal advection in an air quality model (CAMx V6.10), Geosci. Model Dev., 16, 4367–4383, https://doi.org/10.5194/gmd-16-4367-2023, 2023. a

Cariolle, D., Moinat, P., Teyssèdre, H., Giraud, L., Josse, B., and Lefèvre, F.: ASIS v1.0: an adaptive solver for the simulation of atmospheric chemistry, Geosci. Model Dev., 10, 1467–1485, https://doi.org/10.5194/gmd-10-1467-2017, 2017. a, b, c, d

Colella, P. and Woodward, P. R.: The piecewise parabolic method (PPM) for gas-dynamical simulations, J. Comput. Phys., 11, 38–39, https://doi.org/10.1016/0021-9991(84)90143-8, 1984. a, b, c, d, e, f, g, h, i, j

Després, B. and Lagoutière, F.: Un schéma non linéaire anti-dissipatif pour l'équation d'advection linéaire, CR Acad. Sci. I – Math., 328, 939–943, https://doi.org/10.1016/S0764-4442(99)80301-2, 1999. a, b, c

Emmons, L. K., Walters, S., Hess, P. G., Lamarque, J.-F., Pfister, G. G., Fillmore, D., Granier, C., Guenther, A., Kinnison, D., Laepple, T., Orlando, J., Tie, X., Tyndall, G., Wiedinmyer, C., Baughcum, S. L., and Kloster, S.: Description and evaluation of the Model for Ozone and Related chemical Tracers, version 4 (MOZART-4), Geosci. Model Dev., 3, 43–67, https://doi.org/10.5194/gmd-3-43-2010, 2010. a

Esentürk, E., Abraham, N. L., Archer-Nicholls, S., Mitsakou, C., Griffiths, P., Archibald, A., and Pyle, J.: Quasi-Newton methods for atmospheric chemistry simulations: implementation in UKCA UM vn10.8, Geosci. Model Dev., 11, 3089–3108, https://doi.org/10.5194/gmd-11-3089-2018, 2018. a, b

Freitas, S. R., Rodrigues, L. F., Longo, K. M., and Panetta, J.: Impact of a monotonic advection scheme with low numerical diffusion on transport modeling of emissions from biomass burning, J. Adv. Model. Earth Sy., 4, M01001, https://doi.org/10.1029/2011MS000084, 2012. a

Godunov, S. K.: Finite difference method for numerical computation of discontinuous solutions of the equations of fluid dynamics, Matematičeskij Sbornik, 47, 271–306, https://hal.archives-ouvertes.fr/hal-01620642 (last access: 18 December 2023), 1959. a, b, c, d, e, f, g, h

Harris, L., Chen, X., Putman, W., and Zhou, Linjiong andChen, J.-H.: A Scientific Description of the GFDL Finite-Volume Cubed-Sphere Dynamical Core, Tech. Rep., National Oceanic and Atmospheric Administration, Office of Oceanic and Atmospheric Research, Geophysical Fluid Dynamics Laboratory (U.S.), https://doi.org/10.25923/6nhs-5897, 2021. a

Hertel, O., Berkowicz, R., Christensen, J., and Hoc, O.: Test of two numerical schemes for use in atmospheric transport-chemistry models, Atmos. Environ., 27A, 2591–2611, https://doi.org/10.1016/0960-1686(93)90032-T, 1993. a, b, c

Huijnen, V., Williams, J., van Weele, M., van Noije, T., Krol, M., Dentener, F., Segers, A., Houweling, S., Peters, W., de Laat, J., Boersma, F., Bergamaschi, P., van Velthoven, P., Le Sager, P., Eskes, H., Alkemade, F., Scheele, R., Nédélec, P., and Pätz, H.-W.: The global chemistry transport model TM5: description and evaluation of the tropospheric chemistry version 3.0, Geosci. Model Dev., 3, 445–473, https://doi.org/10.5194/gmd-3-445-2010, 2010. a

Lachatre, M., Mailler, S., Menut, L., Turquety, S., Sellitto, P., Guermazi, H., Salerno, G., Caltabiano, T., and Carboni, E.: New strategies for vertical transport in chemistry transport models: application to the case of the Mount Etna eruption on 18 March 2012 with CHIMERE v2017r4, Geosci. Model Dev., 13, 5707–5723, https://doi.org/10.5194/gmd-13-5707-2020, 2020. a, b, c

Lachatre, M., Mailler, S., Menut, L., Cholakian, A., Sellitto, P., Siour, G., Guermazi, H., Salerno, G., and Giammanco, S.: Modelling SO2 conversion into sulfates in the mid-troposphere with a 3D chemistry transport model: the case of Mount Etna's eruption on 12 April 2012, Atmos. Chem. Phys., 22, 13861–13879, https://doi.org/10.5194/acp-22-13861-2022, 2022. a, b

Lauritzen, P. H., Skamarock, W. C., Prather, M. J., and Taylor, M. A.: A standard test case suite for two-dimensional linear transport on the sphere, Geosci. Model Dev., 5, 887–901, https://doi.org/10.5194/gmd-5-887-2012, 2012. a

Lauritzen, P. H., Ullrich, P. A., Jablonowski, C., Bosler, P. A., Calhoun, D., Conley, A. J., Enomoto, T., Dong, L., Dubey, S., Guba, O., Hansen, A. B., Kaas, E., Kent, J., Lamarque, J.-F., Prather, M. J., Reinert, D., Shashkin, V. V., Skamarock, W. C., Sørensen, B., Taylor, M. A., and Tolstykh, M. A.: A standard test case suite for two-dimensional linear transport on the sphere: results from a collection of state-of-the-art schemes, Geosci. Model Dev., 7, 105–145, https://doi.org/10.5194/gmd-7-105-2014, 2014. a

Lauritzen, P. H., Conley, A. J., Lamarque, J.-F., Vitt, F., and Taylor, M. A.: The terminator “toy” chemistry test: a simple tool to assess errors in transport schemes, Geosci. Model Dev., 8, 1299–1313, https://doi.org/10.5194/gmd-8-1299-2015, 2015. a, b

LeVeque, R. J.: High-Resolution Conservative Algorithms for Advection in Incompressible Flow, SIAM J. Numer. Anal., 33, 627–665, https://doi.org/10.1137/0733033, 1996. a, b, c, d, e, f

Mailler, S.: AdvBench (v1.0.0), Zenodo [code], https://doi.org/10.5281/zenodo.7937121, 2023. a

Mailler, S. and Pennel, R.: ToyCTM v1.0.1, Zenodo [code], https://doi.org/10.5281/zenodo.10018706, 2023. a, b

Mailler, S., Menut, L., di Sarra, A. G., Becagli, S., Di Iorio, T., Bessagnet, B., Briant, R., Formenti, P., Doussin, J.-F., Gómez-Amo, J. L., Mallet, M., Rea, G., Siour, G., Sferlazzo, D. M., Traversi, R., Udisti, R., and Turquety, S.: On the radiative impact of aerosols on photolysis rates: comparison of simulations and observations in the Lampedusa island during the ChArMEx/ADRIMED campaign, Atmos. Chem. Phys., 16, 1219–1244, https://doi.org/10.5194/acp-16-1219-2016, 2016. a

Mailler, S., Pennel, R., Menut, L., and Lachâtre, M.: Using the Després and Lagoutière (1999) antidiffusive transport scheme: a promising and novel method against excessive vertical diffusion in chemistry-transport models, Geosci. Model Dev., 14, 2221–2233, https://doi.org/10.5194/gmd-14-2221-2021, 2021. a

Mailler, S., Pennel, R., Menut, L., and Cholakian, A.: Scripts to reproduce the results of gmd-2023-78. In Geoscientific Model Development (v1.0.1), Zenodo [code], https://doi.org/10.5281/zenodo.10018761, 2023. a

Martin, R. V., Eastham, S. D., Bindle, L., Lundgren, E. W., Clune, T. L., Keller, C. A., Downs, W., Zhang, D., Lucchesi, R. A., Sulprizio, M. P., Yantosca, R. M., Li, Y., Estrada, L., Putman, W. M., Auer, B. M., Trayanov, A. L., Pawson, S., and Jacob, D. J.: Improved advection, resolution, performance, and community access in the new generation (version 13) of the high-performance GEOS-Chem global atmospheric chemistry model (GCHP), Geosci. Model Dev., 15, 8731–8748, https://doi.org/10.5194/gmd-15-8731-2022, 2022. a

Menut, L., Bessagnet, B., Briant, R., Cholakian, A., Couvidat, F., Mailler, S., Pennel, R., Siour, G., Tuccella, P., Turquety, S., and Valari, M.: The CHIMERE v2020r1 online chemistry-transport model, Geosci. Model Dev., 14, 6781–6811, https://doi.org/10.5194/gmd-14-6781-2021, 2021. a

Nair, R. D. and Lauritzen, P. H.: A class of deformational flow test cases for linear transport problems on the sphere, J. Comput. Phys., 229, 8868–8887, https://doi.org/10.1016/j.jcp.2010.08.014, 2010. a

Nakamura, N.: Two-dimensional mixing, edge formation, and permeability diagnosed in area-coordinate, J. Atmos. Sci., 53, 1524–1537, https://doi.org/10.1175/1520-0469(1996)053<1524:TDMEFA>2.0.CO;2, 1996. a, b, c

O'Connor, F. M., Johnson, C. E., Morgenstern, O., Abraham, N. L., Braesicke, P., Dalvi, M., Folberth, G. A., Sanderson, M. G., Telford, P. J., Voulgarakis, A., Young, P. J., Zeng, G., Collins, W. J., and Pyle, J. A.: Evaluation of the new UKCA climate-composition model – Part 2: The Troposphere, Geosci. Model Dev., 7, 41–91, https://doi.org/10.5194/gmd-7-41-2014, 2014. a

Pozzoli, L., Bey, I., Rast, S., Schultz, M. G., Stier, P., and Feichter, J.: Trace gas and aerosol interactions in the fully coupled model of aerosol-chemistry-climate ECHAM5-HAMMOZ: 1. Model description and insights from the spring 2001 TRACE-P experiment, J. Geophys. Res.-Atmos., 113, D07308, https://doi.org/10.1029/2007JD009007, 2008.  a

Prather, M. J.: Numerical advection by conservation of second-order moments, J. Geophys. Res.-Atmos., 91, 6671–6681, https://doi.org/10.1029/JD091iD06p06671, 1986. a

Putman, W. M. and Lin, S.-J.: Finite-volume transport on various cubed-sphere grids, J. Comput. Phys., 227, 55–78, https://doi.org/10.1016/j.jcp.2007.07.022, 2007. a

Seinfeld, J. H. and Pandis, S. N.: Atmospheric chemistry and physics: From air pollution to climate change, Wiley-Interscience, ISBN 0-471-17815-2, 1997. a

Strang, G.: On the construction and comparison of difference schemes, SIAM J. Numer. Anal., 5, 506–517, https://doi.org/10.1137/0705041, 1968. a

Timmermans, R., van Pinxteren, D., Kranenburg, R., Hendriks, C., Fomba, K., Herrmann, H., and Schaap, M.: Evaluation of modelled LOTOS-EUROS with observational based PM10 source attribution, Atmospheric Environment: X, 14, 100173, https://doi.org/10.1016/j.aeaoa.2022.100173, 2022. a

Van Leer, B.: Towards the ultimate conservative difference scheme. IV. A new approach to numerical convection, J. Comput. Phys., 23, 276–299, https://doi.org/10.1016/0021-9991(77)90095-X, 1977. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o

Walcek, C. J.: Minor flux adjustment near mixing ratio extremes for simplified yet highly accurate monotonic calculation of tracer advection, J. Geophys. Res., 105, 9335–9348, https://doi.org/10.1029/1999JD901142, 2000. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s

Waruszewski, M., Kühnlein, C., Pawlowska, H., and Smolarkiewicz, P. K.: MPDATA: Third-order accuracy for variable flows, J. Comput. Phys., 359, 361–379, https://doi.org/10.1016/j.jcp.2018.01.005, 2018. a

Zhao, S., Russell, M. G., Hakami, A., Capps, S. L., Turner, M. D., Henze, D. K., Percell, P. B., Resler, J., Shen, H., Russell, A. G., Nenes, A., Pappin, A. J., Napelenok, S. L., Bash, J. O., Fahey, K. M., Carmichael, G. R., Stanier, C. O., and Chai, T.: A multiphase CMAQ version 5.0 adjoint, Geosci. Model Dev., 13, 2925–2944, https://doi.org/10.5194/gmd-13-2925-2020, 2020. a

Download
Short summary
We show that a new advection scheme named PPM + W (piecewise parabolic method + Walcek) offers geoscientific modellers an alternative, high-performance scheme designed for Cartesian-grid advection, with improved performance over the classical PPM scheme. The computational cost of PPM + W is not higher than that of PPM. With improved accuracy and controlled computational cost, this new scheme may find applications in chemistry-transport models, ocean models or atmospheric circulation models.