A standard test case suite for two-dimensional linear transport on the sphere

,


Introduction
A basic building block in any fluid dynamics solver is the transport operator that approximates the evolution of the bulk motion of a scalar.Despite intense research in transport schemes intended for global modeling on the sphere, only test 1 of the widely used test case suite by (Williamson et al., 1992) seems to be the standard test whereas other (newer) test cases are, in general, only optionally used.Test 1 in (Williamson et al., 1992) is referred to as the solid-body advection test case and the exact solution is simply the translation of the initial condition so that the center of the distribution follows a Figures great-circle.The flow field is non-divergent and does not challenge the transport operator with respect to deformation or divergence.In the last decade other non-divergent global test cases have been proposed such as static (Nair and Machenhauer, 2002) and moving vortices (Nair and Jablonowski, 2008) test cases that include deformation.Also for these tests the analytical solution is known at all times.Scheme developers do, in general, not publish results for all test cases and, perhaps more importantly, they often choose different parameter settings making it more difficult to compare results for different schemes.A purpose of this paper is to provide specific guidelines for test case setup in terms of parameters, resolution, time-step, and diagnostics.Perhaps more challenging analytical wind fields were recently proposed by (Nair and Lauritzen, 2010).The Lagrangian fluid parcels follow complex trajectories (not greatcircles or small circles) making it harder to compute the analytical solution throughout the simulation.Following (LeVeque, 1996) the flow has a "time-reversing" component so that after one period the exact solution equals the initial condition.Half way through the simulation, however, the initial distributions are deformed into thin filaments and an "overlaid" translational flow transports the filaments as they deform.This problem is very challenging.A divergent wind field is proposed in (Nair and Lauritzen, 2010) as well which is in contrast to most idealized wind fields in the literature.The combination of both divergent and deformational flow constitutes a more realistic atmospheric/ocean transport than, for example, solid-body advection flow.
The idealized transport experiments listed above are all based on a single tracer and accuracy is quantified in terms of conventional errors norms, i.e. quantifications of the differences between the analytical (exact) and numerically computed solutions.In some geophysical fluid dynamics problems, such as the transport of long-lived species in the stratosphere and aerosol-cloud interactions (Ovtchinnikov and Easter, 2009), it is not only important that individual species/tracers are transported accurately but also the maintenance of pre-existing functional relations between species/tracers is important.Such models also cannot accept non-physical transport or redistribution of tracer that is not accompanied by resolved motion of air masses.Following Figures

Back Close
Full Lagrangian fluid parcels interrelations between tracers are conserved, however, any non-Lagrangian scheme will almost certainly perturb such relations.Nevertheless, Lagrangian schemes in realistic divergent flows must eventually combine parcels or create new ones, and that process will not likely preserve the relationships.Numerical errors that perturb pre-existing functional relations can resemble "real" mixing similar to what is observed in nature when mixing occurs (hereafter referred to as "real mixing") or the truncation errors can introduce unmixing (that is, spurious mixing) (Thuburn and Mclntyre, 1997).A quantification and classification of mixing between interrelated tracers was recently proposed in (Lauritzen and Thuburn, 2011).For a more extensive overview of test cases for global models and desirable properties for transport schemes intended for atmospheric modeling see, for example, the recent book chapter by (Lauritzen et al., 2011).The purpose of this paper is to propose a minimal and challenging test case suite with specific guidelines on the implementation and diagnostics, thereby facilitating intercomparison of schemes and establishing a benchmark data base for future developers.In the derivation of the test suite we sought to minimize the workload on transport scheme developers while evaluating their schemes in terms of a wide range of quantitative measures of accuracy considered important for geophysical fluid dynamics.Therefore we assume that modelers have already tested their schemes in simpler settings such as with solid-body and static/moving vortices test cases and we do not repeat those tests here.
Almost any test case suite could be extended to include more tests that could provide more insights into specific aspects of accuracy particularly useful for some classes of schemes and applications.For example, (Ullrich et al., 2012) found it insightful (for evaluating higher-order approximations to Lagrangian cell sides) to transport a constant using the initial condition wind wield for the shallow water test case 3 of (Williamson et al., 1992).Similarly one could use the actual observed winds in test case 7 in (Williamson et al., 1992) to generate more complex structures in the tracer field.To produce filaments that eventually become sub-grid-scale in a context where the analytical solution Introduction

Conclusions References
Tables Figures

Back Close
Full is known, one may use the moving vortices test case (Nair and Jablonowski, 2008) and run it out further than the 12 days suggested in the test case description (Pudykiewicz, 2011).The test case suite presented in this paper is not meant to be exhaustive and developers usually have preferred idealized tests specific for their application.It is the intent of this paper to present a minimal test bed based on just two analytic wind fields and four initial conditions that address a wide range of accuracy aspects, challenge both Lagrangian and Eulerian schemes with realistic conditions typical of 3-D flows, and make it straight forward to compare results from different schemes since we provide detailed instructions on test case setup and diagnostics.In doing so we believe this test case suite provides new insights into accuracy beyond the much simpler and most widely used standard solid-body advection test case and associated standard error norms.
The test case proposal is organized as follows.In Sect. 2 the transport equation(s) is introduced followed by formulas for the analytical initial conditions and wind fields.The section is concluded with discretization details such as "definition" of resolution and maximum CFL (Courant-Friedrichs-Lewy) numbers.The actual test case setup is given in Sect. 3 and it is divided into six categories designed to assess numerical convergence rates, "minimal" resolution, filament preservation, transport of discontinuous distributions, maintenance of preexisting non-linear functional relations, and transport under divergent flow conditions.Information important for discussing computational efficiency is proposed in Sect. 4. The paper in concluded with a summary of the test cases.Introduction

Conclusions References
Tables Figures

Back Close
Full 2 General problem formulation

Transport equation(s)
Consider a transport scheme that approximates the solution to the continuity equation for a passive (does not feed back on the flow) and inert (no sources or sinks) tracer, where ρ is the fluid density, V is the two-dimensional flow velocity vector, and φ is the tracer mixing ratio per unit mass.Note that the discretized scheme is not necessarily based on the continuity equation written in flux-form as in Eq. ( 1) but could also be based on the advective form or cell-integrated Lagrangian form where A(t) is a Lagrangian area and D/Dt is the total or material derivative.If tracer density ρφ (flux-form Eq. 1 or cell-integrated Lagrangian Eq. 3) and not mixing ratio φ (advective form Eq. 2) is the prognostic variable, one needs to "extract" the mixing ratio φ from ρφ, which obviously requires the solution to the continuity equation for fluid density ρ (see, e.g., Nair and Lauritzen, 2010 for details).In this test case suite the mixing ratio φ is used for all diagnostics/analysis and not tracer density ρφ.
Define the discrete transport operator T that advances the numerical solution for φ in time where n is the time-level index, k is the index for the grid cell/point, and H is the set of indices defining the "halo" or computational stencil required by T .To compute the numerical solution φ n+1 k an initial condition, a prescribed velocity field and (if applicable) the solution to the continuity equation for fluid density ρ is required.

Initial conditions
For fluid density (if needed) the initial condition is ρ(t = 0) = 1 for all test cases.Four initial conditions for mixing ratio φ are used and defined in sub-sections below.We use one infinitely smooth (C ∞ ) initial condition, one quasi-smooth, one discontinuous, and one non-linearly correlated with the quasi-smooth initial condition.It may be argued that the discrete initial conditions should be as consistent as possible with the numerical method.For example, a finite-volume method is usually based on cell-averaged prognostic variables and the initial condition in cell k, φ k , should be obtained by integrating the continuous initial condition φ over the k-th control volume.Similarly for methods that preserve and transport internal moments of the tracer distribution (e.g., Prather, 1986) should initialize these moments by integrating over the continuous starting distribution.Standard practice, however, is to to use the value of the continuous initial condition evaluated at the centroid of the control volume as representative for the cell averaged value and higher-order moments (if applicable) are zero.It has been shown for finite-volume schemes that standard error norms may vary significantly when using point or cell-averaged values for initial conditions and for computing error norms (e.g., Lauritzen et al., 2010;Zerroukat et al., 2002).However, the conclusions drawn from the results are independent of the choice of exact solution (cell average versus grid-point value) as long as the schemes are compared with the same choice for exact solution in a consistent manner.Therefore the initial condition and exact solution are based on grid-point values at the centroid of the grid cell for finite-volume methods and at quadrature/finite-difference points for basis-function/grid-point methods.All four initial conditions for φ are based on two distributions centered about (λ i ,θ i ), i = 1,2: where λ is longitude and θ is latitude in radians.The distributions are symmetrically placed in the flow field to assess the symmetry of the numerically computed solution.

Gaussian hills
Smooth Gaussian surfaces can be defined as follows: where the height1 and width are determined by h max = 0.95 and b = 5, respectively (Levy et al., 2007).The absolute Cartesian coordinates (X,Y,Z) and spherical (λ,θ) coordinates are related through where radius R is the radius of the sphere.The coordinates for the center of the Gaussian distribution (X i ,Y i ,Z i ) is computed by inserting (λ i ,θ i ) into Eq.( 8) and evaluating the right-hand side.
The Gaussian hills distribution is defined as the sum of the two Gaussian hills h 1 and h 2 (Eq.7), and is graphically shown on Fig. 1a.Note that φ (gh) is infinitely smooth (C ∞ ).Introduction

Conclusions References
Tables Figures

Back Close
Full

Cosine bells
Similarly two symmetrically located cosine bells are defined as follows: where the amplitude h max = 1, base radius r = R/2, and great-circle distance between (λ,θ) and the center (λ i ,θ i ) is r i = r i (λ,θ), with The initial condition φ consists of a background value b and two cosine bells defined above where the background value is b = 0.1 and amplitude c = 0.9 such that φ ∈ [0.1,1.0](see Fig. 1b).

Conclusions References
Tables Figures

Back Close
Full

"Correlated" cosine bells
An initial distribution that is nonlinearly "correlated" with the cosine bells initial condition is defined as The nonlinear functional relation is given by where For a contour plot of the correlated cosine bells see Fig. 1d.

Wind fields
In this test case suite we use two deformational wind fields: one non-divergent and one divergent.The components of the non-divergent velocity vector V (λ,θ,t) and the stream function are given by

Conclusions References
Tables Figures

Back Close
Full respectively, where λ = λ − 2πt/T .In non-dimensional units T = 5 and R = 1.An "Earth"-like dimensionalization of the wind fields may be obtained by setting T = 12 days and R = 6.3172 × 10 6 m.Schemes based on characteristics (typically Lagrangian schemes) may use the algorithm given in (Nair and Lauritzen, 2010) for the computation of parcel trajectories.
When either of the initial conditions given in Sect.2.2 are transported by the nondivergent wind field, they are deformed into thin filaments half way through the simulation and these are simultaneously being transported Eastward by the solid-body component of the flow (see Figs. 2 and 3).At maximum deformation the filaments are approximately 10 • wide when using the cosine bells initial condition.
To challenge schemes under divergent flow conditions we use the following wind field (Nair and Lauritzen, 2010, their case-3 with a "constant background wind field"): where The exact solution for all tests is known at t = T and it is identical to the initial condition We do not have an exact solution throughout the simulation when using either the non-divergent or divergent flow field.Note that the first part of the simulation (t ∈ [0,T/2]) is typical of atmospheric/oceanographic flows in that features collapse to smaller scales whereas the second part (t ∈ [T/2,T ]), in which the reverse occurs, is atypical of atmospheric/oceanographic flows though convenient for getting a problem with an exact solution.The background mean flow ensures that errors, in general, do not cancel when the deformational part of the flow reverses.

Discretization details
We specify resolution in terms of average grid-spacing in degrees at the Equator of the sphere ∆λ.For methods based on quadrature methods the "average resolution" should be specified in terms of mean distance between quadrature points.We define the (maximum) CFL number as where ∆t is the time-step and U max is the maximum zonal wind speed.For the nondivergent flow the non-dimensional and dimensional ("Earth") U max are given by respectively.This definition of CFL number obviously does not emphasize local CFL numbers (in particular for non-isotropic grids); it is defined to facilitate comparison of maximum CFL numbers across discretization grids.
The time-step ∆t should be a "typical/practical" time-step for performing tracer transport with the scheme in question.However, investigating accuracy as a function of 200 Introduction

Conclusions References
Tables Figures

Back Close
Full time-step is also of interest.For example, if the transport scheme permits long timesteps (CFL > 1) it is advised to run the test suite with an "Eulerian" time-step (CFL ≤ 1) as well.
Often limiters/filters are applied to render the numerically computed solution physically realizable.These may be shape-preserving, monotone, and/or non-oscillatory limiters/filters.If schemes have a limiter/filter option, the test suite should be run both without and with limiters/filters.
Accuracy is assessed in terms of several diagnostics.First if all, we use standard error norms that are defined in Appendix A. These require knowledge of the "true" (analytic) solution and are therefore computed at time t = T when the true solution is known.Secondly we use recently proposed mixing diagnostics (Sect.3.5; Appendix B and C) as well as a novel filament preservation diagnostic (Sect.3.3).As these diagnostics do not require an analytical solution we compute them at the time of maximum deformation (t = T/2) before the flow "reversal" which is less physical.
For reference purposes we provide results using the CSLAM (Conservative Semi-Lagrangian Multi-tracer) scheme (Lauritzen et al., 2010) on the cubed-sphere grid.The CSLAM configuration used here is described in detail in (Nair and Lauritzen, 2010).

Test cases
The diagnostics/test cases are designed to assess: 1. numerical order of convergence, 6. ability of transport scheme to operator to deal with divergent flows (Nair and Lauritzen, 2010).
Each category is discussed in separate sections below.

Numerical order of convergence: Gaussian hills
This test is designed to assess the formal (or "optimal") order of convergence of the scheme under quasi-realistic flow conditions on the sphere.This is done as follows.Standard error norms using the Gaussian hills initial condition (Eq.9) and nondivergent wind field (Eqs.18 and 19) at resolutions ranging from approximately ∆λ = 3 • to ∆λ = 0.3 • for fixed CFL number are computed.The choice of resolutions should provide enough data points on a "convergence plot" (e.g., log( 2 ) as a function of log(N)) in the resolution interval of interest, to generate a "credible" estimate of numerical rate of convergence.For example, the following resolutions could be used: The runs should be performed without any limiting/filtering and (if applicable) also with limiters/filters enforcing shape-preservation, monotonicity and/or nonoscillatoriness in the numerically computed solution.
These simulations with infinitely smooth (Gaussian hills) initial conditions should provide a numerical estimate of the "optimal" numerical convergence rate of the scheme.A way to estimate numerical (empirical) convergence rates K 2 and K ∞ , for 2 and ∞ , respectively (see Fig. 4), is to perform a least-squares linear regression of the form (Harris et al., 2010) (26)

"Minimal" resolution ∆λ m : cosine bells
In many geophysical fluid dynamics applications using state-of-the-art physical parameterization packages, increases in horizontal resolution comes at significant computational cost.It is therefore of interest to assess the absolute error in addition to convergence rates.To do that we repeat the experiment described in Sect.3.1 but with cosine

GMDD Introduction Conclusions References
Tables Figures

Back Close
Full bells initial condition (Eq.11) to find the "minimal" resolution.We define the "minimal" resolution ∆λ m as the ∆λ-value for which 2 is approximately 0.033, when using an unlimited scheme and the cosine bells (Eq.11) initial condition (the CFL number used for defining ∆λ m should be one typically used by the scheme).A convergence plot can conveniently be used to find the "minimal" resolution by finding the intersection between the horizontal line 2 = 0.033 and the convergence curve for 2 (see Fig. 5).The quasi but not infinitely smooth initial conditions (Cosine bells instead of Gaussian hills) are used in order to challenge the schemes with respect to weak non-smoothness.
The "minimal" resolution ∆λ m will be used in the remaining test cases.The choice of threshold for ∆λ m is based on results for CSLAM (a resolution for which the thin filaments are marginally resolved).The "minimal" resolution (as defined here) for CSLAM is ∆λ = 1.5 • and ∆λ = 1 • when using a time-step of T/120 (maximum CFL is approximately 5.2) and T/600 (maximum CFL is approximately 1.0).

"Filament" preservation diagnostic f : cosine bells
Realistic flows often deform distributions into thin filaments which, in general, are challenging to represent by Eulerian and semi-Lagrangian transport schemes that use a fixed grid in space (e.g., Behrens et al., 2000).A measure of how well a transport scheme preserves gradients, in particular, thin filaments is relevant for many tracer applications (e.g.transport of long-lived tracers such as chemical species in the stratospheric vortices).Filaments are created when material surfaces stretch and gradients increase.When the thickness of the filaments reach the scale at which molecular diffusion (or some other diffusive process) becomes important the filaments are no longer preserved but gradients are eroded.For the flow and initial conditions considered here the filaments should, for all practical purposes, be preserved by the transport scheme as the physical scale of the filaments is approximately 10 • at maximum deformation.
We do therefore not assess how transport schemes represent the filament erosion process that appears in nature since those "diffusive" processes take place at scales several magnitudes below 10 Full moving vortices test case of (Nair and Jablonowski, 2008) and extend the simulation time so that the filaments are stretched to a level where such processes are important and/or change the parameters in the (Nair and Lauritzen, 2010) flow field to increase the amount of deformation (see, e.g., Kent et al., 2012;Pudykiewicz, 2011).
The "filament" preservation diagnostic is formulated as follows.Define A(τ,t) as the spherical area for which the spatial distribution of the tracer φ(λ,θ) satisfies at time t, where τ is the threshold value.For a non-divergent flow field and a passive and inert tracer φ, the area A(τ,t) is invariant in time.
The discrete definition of A(τ,t) is where ∆A k is the spherical area for which φ k is representative, K is the number of grid cells, and G is the set of indices For Eulerian finite-volume schemes ∆A k is the area of the k-th control volume.For Eulerian grid-point schemes a control volume for which the grid-point value is representative must be defined.Similarly for fully Lagrangian schemes based on point values (parcels) control volumes for which the point values are representative must be defined.Note that the "control volumes" should span the entire domain without overlaps or "cracks" between them.
Define the filament preservation diagnostic For infinite resolution (continuous case) and a non-divergent flow, f (τ,t) is invariant in time: f (τ,t = 0) = f (τ,t) = 100 for all τ.At finite resolution, however, the filament 204 Introduction

Conclusions References
Tables Figures

Back Close
Full diagnostic even for an exact scheme should not necessarily be preserved since the solution must be truncated to the discrete grid.That said, usually the numerical truncation errors are much larger than the grid truncation error at least at moderate resolutions.The experimental setup is as in Sect.3.2, that is, use the non-divergent wind field (Eq.18 and 19) and cosine bells initial condition (Eq.11).At half time, t = T/2, the filament preservation diagnostic f (τ,t = T/2) is computed at 19 equi-distant discrete intervals (τ = 0.10, 0.15, 0.15, 0.20, . . ., 0.95, 1.00) without and (if applicable) with limiters/filters at ∆λ = 1.5 • , ∆λ = 0.75 • as well as at the "minimal" resolution ∆λ = ∆λ m .
The threshold value for which f (t = T/2) is less than,for example 80, is a measure for how well filaments are preserved.Numerical diffusion will tend to decrease f for large τ values (maxima decrease) and increase f for low τ values (gradients are "smeared").An "extreme" situation is shown on Fig. 6a where f is plotted as a function of τ for the highly diffusive 1st-order version of CSLAM.This much improves when using the higher-order version of CSLAM (Fig. 6b).Note that the non-shape-preserving versions of CSLAM produce values of f less than 100.0 for low threshold values (τ < 0.1).This also indicates an error in tracer transport due to undershoots (φ < 0.1), which are not represented in the l f diagnostic.

Transport of "rough" distribution: slotted-cylinders
To challenge shape-preserving filters/limiters (if applicable) we use discontinuous initial conditions, that is, standard error norms for the simulated solution at t = T using the slotted cylinders initial condition and non-divergent winds (Eqs.18 and 19) are computed using the transport scheme without and (if applicable) with limiters/filters at resolutions ∆λ = 1.5 • , ∆λ = 0.75 • as well as at the "minimal" resolution ∆λ m .Contour plots of the solution at t = T/2 and t = T (Fig. 7) using a contour interval of 0. 3.5 Preservation of pre-existing functional relation: cosine bells and correlated cosine bells In the tests described in the previous sections the accuracy is assessed in a singletracer setup.Now we consider two tracers that are both advected by the same nondivergent flow field (Eqs.18 and 19).The initial conditions for the two tracers is the cosine bells initial condition (Eq.11) and correlated cosine bells (Eq.13), respectively (see Fig. 1b,d).The mixing ratio of the two tracers are referred to as χ and ξ.Following Lagrangian parcels any functional relation between tracers should mathematically be preserved at all times and hence any deviation from the pre-existing functional relation between the tracers is essentially numerical errors introduced by the transport scheme.
Note that the "ideal" scheme could be a scheme that does not exactly preserve preexisting functional relations but for which the numerical errors are less than physical diffusive processes in nature.
In any case transport schemes should not disrupt functional relations in unphysical ways.Numerical errors that perturb such relations essentially introduce mixing or unmixing between the tracers.(Lauritzen and Thuburn, 2011) provides a discussion of the physical importance of transport schemes not disrupting tracer interrelationships in unphysical ways with special focus on non-linear chemistry.The numerical errors that perturb pre-existing functional relations between tracers will be referred to as numerical mixing or simply mixing in this paper (one could equally well use terminology such as tracer variance dissipation instead of mixing).In nature such processes that change the correlation between two tracers come about through diffusive processes, and, for reactive tracers, through chemical reactions between tracers.The purpose of this test is to quantify the amount of mixing and the physical realizability of the mixing that a scheme introduces through truncation errors.Note that preserving correlations is, however, no guarantee for accuracy as one may design schemes that satisfy tracer interrelations but are otherwise inaccurate; as formulated by (Thuburn and Mclntyre, 1997): "shaping two tracer fields the same way does not imply shaping them the right way."

Conclusions References
Tables Figures

Back Close
Full Scatter plots, where tracer 1 (χ using cosine bells initial condition) and tracer 2 (ξ using correlated cosine bells initial condition) are plotted against each other, are used to quantify the numerical mixing or unmixing introduced by the scheme (see Fig. 8).As discussed in (Thuburn and Mclntyre, 1997) no Eulerian scheme can exactly preserve pre-existing nonlinear relations between two tracers and hence scatter points (χ k ,ξ k ) will, in general, deviate from the pre-existing functional relation curve ψ as the simulation evolves.The way in which the scatter points deviate from the non-linear ψ-curve has implications for the character of the numerical mixing that the transport scheme introduces.For this test it is crucial that features collapse in scale and we therefore consider scatter plots using prognosed mixing ratios at half time (t = T/2) when the initial condition has deformed into thin filaments and collapse to smaller scales compared to the initial condition.
Following (Lauritzen and Thuburn, 2011) the numerical mixing (deviation of scatter points (χ k ,ξ k ) from ψ-curve) is classified into three categories: -"Real" mixing: Numerical mixing that resembles "real" mixing (e.g., Thuburn and Mclntyre, 1997) when scatter points move to the concave side of ψ.All other deviations from the pre-existing functional curve is spurious unmixing which is accounted for in two separate categories.
-"Range-preserving" umixing: Numerical unmixing within the range of the initial data, that is, scatter points are shifted to the convex side of the pre-existing functional relation or below the convex hull but not outside the range of the initial data.
-Overshooting: Numerical unmixing that is not "range-preserving" unmixing which for this specific test case setup is (χ ,ξ) / ∈ The three diagnostics that quantitatively account for numerical mixing that resembles "real" mixing, "range-preserving" umixing and overshooting are referred to as r , u , and o , respectively, and are formally defined in Appendix C. For more discussion on numerical mixing and the physical reasoning behind the classification of the mixing see (Lauritzen and Thuburn, 2011).Note that knowledge of the exact solution is not needed for the computation of the mixing diagnostics.
Using the non-divergent flow field we compute the mixing diagnostics ( r , u , o ) half way through the simulation t = T/2 using two non-linearly correlated tracer distributions χ = φ (cb) and ξ = φ (ccb) as initial conditions (cosine bells and correlated cosine bells) at resolutions ∆λ = 1.5 • , ∆λ = 0.75 • and ∆λ m using the unlimited and (if applicable) limited/filtered scheme.The scatter plots, that is, the mixing ratio of one tracer (with cosine initial conditions) against the other (with non-linearly correlated cosine bells initial condition) at these resolutions are shown in Fig. 8.
It is noted that transport schemes can be designed to preserve linear pre-existing functional relations.That is, a scheme will preserve linear correlations between species/tracers if the transport operator T satisfies where A and B are constants (Lin and Rood, 1996;Thuburn and Mclntyre, 1997).It is assumed that schemes have already been tested with respect to preservation of linear correlations without and (if applicable) with limiters/filters.

Transport under divergent flow conditions: cosine bells
Most idealized test cases are formulated in terms of non-divergent wind fields.Since realistic flows are divergent it should be demonstrated that the transport operator can handle divergent winds.We repeat the experiment described in Sect.3.4 using the divergent wind field (see Eqs. 21 and 22), cosine bells initial conditions (Eq.11), and the same time-steps.Solutions using CSLAM are shown on Fig. 9.

GMDD Introduction Conclusions References
Tables Figures

Computational cost
Measures that are likely to impact computational efficiency across platforms are also documented.Ultimately computational efficiency is machine dependent, in particular, some computing architectures favor certain types of algorithms over others.Hence it is impossible to define a universal and objective measure of efficiency that is applicable for all computing platforms.Below is a non-exhaustive list of aspects of algorithms that are likely to impact efficiency size of halo/stencil H used to update a cell/grid-point value, number of communications (in parallel setup) that are necessary during one full tracer time-step ∆t, number of integral/functional evaluations (if applicable), maximum CFL for which the transport scheme is stable, amount of information (if any) that can be re-used to transport additional tracers (multi-tracer efficiency).

Summary
Below is a summary of the proposed test case suite.In terms of implementation work using traditional/conventional error norms as well as novel filament-preservation and mixing diagnostics.For convenience the standard error norms i , i = 2,∞, φ min and φ max are computed at the end of the simulation t = T when the exact solution is known (i.e. it equals the initial condition).All mixing diagnostics i , i = r,u,o, and the filament diagnostic f (they do not require knowledge of the analytical solution to the transport equation) are computed half way through the simulation at t = T/2 when the fields are most deformed.
For the non-divergent flow field (Eqs.18 and 19) we: 1. Numerical order of convergence Compute numerical convergence rates K i for i , i = 2,∞, for the resolution range approximately ∆λ = 3 • to ∆λ = 0.3 • using Gaussian initial conditions for the unlimited and (if applicable) limited/filtered scheme (Sect.3.1).
4. "Rough" distribution Compute i , i = 2,∞, φ min and φ max at resolutions ∆λ = 1.5 • , ∆λ = 0.75 • , and ∆λ = ∆λ m for the slotted-cylinder initial conditions (Eq.12) using the unlimited and (if applicable) the shape-preserving scheme (Sect.3.4).where φ T , φ 0 are, respectively, the exact/analytical solution and its initial value, ∆φ 0 is the difference between maximum and minimum value of the initial condition, and the global integral I is defined as follows,  where

Mixing diagnostics
For this particular test case setup R χ = 0.9, R ξ = 0.792, and the "root" χ (root) k is given by where

Appendix C Numerical mixing diagnostics
A For the two-tracer test (Sect.3.5) three mixing diagnostics are used and defined below (Lauritzen and Thuburn, 2011).

Mixing that resembles "real" mixing
"Real" mixing is defined as numerical mixing that resembles "real" mixing in that values are shifted to the concave side of the pre-existing functional relation only (area A on Fig. 1) where K is the total numbers of cells/points in the domain, ∆A k is the spherical area of grid cell k and A is the total area of the domain, A = K k=1 ∆A k .The distance function d k is the shortest normalized distance between the numerically computed scatter point (χ k ,ξ k ) and the preexisting functional curve within the range of the initial conditions.For the quadratic functional relation ψ given in Eq. ( 14) with coefficients (Eq.15), the explicit formula for d k is given in Appendix B. The domain A ("convex hull") is shown on Fig. 1 and is mathematically defined as (min ,χ (max) ] and where F is the straight line that connects (χ (min) ,ξ (max) ) and (χ (max) ,ξ (min) ).Any other mixing (i.e.scatter points not in A) is numerical unmixing that is accounted for in two distinct diagnostics defined next.
"Range-preserving" unmixing "Range-preserving" unmixing is defined as numerical unmixing within the range of the initial data, that is, scatter points are shifted to the convex side of the preexisting functional relation or below the convex hull but not outside the range of the initial data where B are the dark shaded areas in Fig. 1 defined by Note that the shape-preservation constraint is not necessarily enough to guarantee u = 0 since the scheme must be semi-linear and monotone according to (Harten, 1983) to guarantee u = 0 (Thuburn and Mclntyre, 1997).Only first-order schemes will satisfy these constraints (Godunov, 1959).Introduction

Conclusions References
Tables Figures

Back Close
Full ∞ , Gaussian hills Fig. 4. Convergence plots for 2 (first column), and ∞ (second column), respectively, computed with CSLAM with Gaussian hills initial conditions.The keys with "CFL 5.5" and "CFL 1.0" refer to simulations using a time-step of T/120 and T/600, respectively.The keys with the word filter in them refer to simulations using a shape-preserving filter.The upper and lower heavy lines on each plot correspond to the slopes of second-and third-order convergence rates, respectively.Introduction

Conclusions References
Tables Figures

Back Close
Full    Similarly for the area labeled with B (defined in Eq.C4) it is categorized as range-preserving unmixing.The remaining part of the domain is referred to as overshooting.The thick solid line is the preexisting non-linear functional relation curve.See text or (Lauritzen and Thuburn, 2011) for details.
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | and T have the same values as for the non-divergent velocity field.The nondivergent flow field (Eqs.18 and 19) is used for all tests except the test described in Sect.3.6 for which the divergent winds are used (Eqs.21 and 22).Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | 2. "minimal" resolution, 3. ability of the transport scheme to preserve filaments, 4. ability of the transport scheme to transport "rough" distributions, 5. ability of the transport scheme to preserve pre-existing functional relations between tracers, Figures Back Close Full Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | 05 in the range [0.0 : 1.1] are shown.Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | [0.1,1.0] 2 .The deviation of the scatter points from the ψ-curve is quantified in terms of a normalized shortest distance between (χ k ,ξ k ) and the ψ-curve referred to as d k .For the specific parabolic non-linear correlation function used here (Eq.14) the normalized distance function d k is given in Appendix B. Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | only two flows fields and four initial conditions are needed.The accuracy is assessed Discussion Paper | Discussion Paper | Discussion Paper | Compute mixing diagnostics i , i = r,u,o, for the two nonlinearly correlated tracers based on cosine bells for the unlimited and (if applicable) limited/filtered scheme at resolutions ∆λ = 1.5 • , ∆λ = 0.75 • , and ∆λ = ∆λ m (Sect.3.5).Discussion Paper | Discussion Paper | Discussion Paper | φ min = min ∀λ,θ (φ) − min ∀λ,θ (φ T ) ∆φ 0 , minimum" distance function d k is defined as the minimal normalized Euclidean distance between the correlation point (χ k ,ξ k ) and the preexisting functional relation curve (χ ,ψ(χ )) within the range of the initial condition d k = L k distance to the initial condition interval [χ min ,χ max ], and the normalized distance function is given by Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

Fig. 1 .Fig. 2 .Fig. 3 .
Fig. 1.Contour plots for the four initial conditions for mixing ratio φ used in this test suite.(a) depicts the infinitely smooth (C ∞ ) initial condition constructed from Gaussian surfaces, (b) the cosine bells initial condition which is C 1 , (c) the non-smooth slotted cylinders initial condition, and (d) is the initial condition which is nonlinearly correlated with (b).

Fig. 8 .Fig. 9 .
Fig. 8. Scatter plots at t = T/2 for two non-linearly correlated species/tracers based on cosine bells initial conditions using first-order version of CSLAM (a and d), standard CSLAM based on bi-parabolic reconstruction functions (b and e) and standard CSLAM with a shape-preserving filter (c and f).First and second row use ∆λ = 1.5• and ∆λ = 0.75• resolutions, respectively.The solid lines mark the boundaries between the areas used to classify the numerical mixing.On each plot the mixing diagnostics r , u and o are given.

Fig. C1 .
Fig. C1.A schematic of the classification of numerical mixing.If a scatter point is located in the area labeled with A (mathematically defined in Eq.C2) it is categorized as "real" mixing.Similarly for the area labeled with B (defined in Eq.C4) it is categorized as range-preserving unmixing.The remaining part of the domain is referred to as overshooting.The thick solid line is the preexisting non-linear functional relation curve.See text or(Lauritzen and Thuburn, 2011) for details.