The Lagrangian chemistry and transport model ATLAS: validation of advective transport and mixing

We present a new global Chemical Transport Model (CTM) with full stratospheric chemistry and Lagrangian transport and mixing called ATLAS (Alfred Wegener InsTitute LAgrangian Chemistry/Transport System). Lagrangian (trajectory-based) models have several important advantages over conventional Eulerian (grid-based) models, including the absence of spurious numerical diffusion, efficient code parallelization and no limitation of the largest time step by the Courant-Friedrichs-Lewy criterion. The basic concept of transport and mixing is similar to the approach in the commonly used CLaMS model. Several aspects of the model are different from CLaMS and are introduced and validated here, including a different mixing algorithm for lower resolutions which is less diffusive and agrees better with observations with the same mixing parameters. In addition, values for the vertical and horizontal stratospheric bulk diffusion coefficients are inferred and compared to other studies. This work focusses on the description of the dynamical part of the model and the validation of the mixing algorithm. The chemistry module, which contains 49 species, 170 reactions and a detailed treatment of heterogeneous chemistry, will be presented in a separate paper.


Introduction
Chemical Transport Models are commonly based on a fixed spatial grid and include an Eulerian advection scheme for the transport of chemical species (e.g.SLIMCAT, Chipperfield, 2006).Although much effort has been put into the development of advection schemes that minimize numerical diffusion (e.g.Prather, 1986), numerical diffusion is in principle unavoidable in Eulerian models.
Correspondence to: I. Wohltmann (ingo.wohltmann@awi.de) In the context of this paper, numerical diffusion is defined as any process that is triggered by the model calculations and behaves like a diffusive process acting on the chemical species in the model.Typically, numerical diffusion will be spurious and unrealistic, e.g.present in every grid cell with the same magnitude.We do not consider processes as numerical diffusion here which are deliberately introduced into the model, e.g. to mimic the observed diffusion.Eulerian approaches suffer from numerical diffusion because they contain an advection step that transfers information between different grid cells.Usually, this step is connected to averaging over several grid points, which introduces diffusion on the spatial scale of the grid resolution at the temporal scale of the model time step.An example is the interpolation of species to trajectory end points in Semi-Lagrangian approaches.
Since diffusion is very small in the stratosphere, numerical diffusion in Eulerian models is always much larger than the real atmospheric diffusion (e.g.Konopka et al., 2005), leading to problems like excessive mixing across atmospheric mixing barriers.Particularly at the low resolutions typically required for long runs (e.g.where the CTM module is coupled to a General Circulation Model, GCM, as part of a Chemistry Climate Model, CCM), numerical diffusion can lead to problems in the representation of the polar vortices and hence the Antarctic ozone hole.With an Eulerian model it is difficult to maintain the steep gradients of species at the boundary of the ozone hole, which are not only relevant for a proper representation of ozone depletion, but also for the radiative and dynamical feedback from the ozone hole on climate processes.We have therefore developed a Lagrangian CTM that can be used as a stand alone model or as module in a CCM system.
Lagrangian models have several advantages over Eulerian models: -There is no numerical diffusion.Actually, some explicit form of diffusion has to be included in Lagrangian models to obtain realistic results, which gives Published by Copernicus Publications on behalf of the European Geosciences Union.
the opportunity to model the real atmospheric diffusion on a physical basis.
-The model architecture allows for an efficient parallelization of code and runs fast on low-cost computers.
-The time step of the integration of the transport equations can be chosen freely.There is no limitation of the largest time step by the Courant-Friedrichs-Lewy criterion.
-Mixing ratios of inert tracers are conserved by design (no negative mixing ratios by transport).
-Transport of additional species does not add to computing time.
In addition, our specific model is free of any grid in the vertical or horizontal.The absence of numerical diffusion, the conservation of mixing ratios and the fact that transport of additional species does not add to computing time are basically due to the fact that the Lagrangian time derivative of mixing ratios is zero in the absence of local sources and sinks, while the Eulerian derivative contains an advection term.Despite these obvious advantages, only very few global Lagrangian chemical transport models with explicit mixing exist so far (e.g.Collins et al., 1997;Fairlie et al., 1999).The most prominent examples are the CLaMS model (McKenna et al., 2002;Konopka et al., 2004Konopka et al., , 2007) ) and the ATTILA model (Reithmeier and Sausen, 2002), which has recently been coupled to the General Circulation Model ECHAM (Stenke et al., 2009).
The main focus of this model version is on a proper representation of the chemistry and transport of the stratosphere, although a basic troposphere is included.Hence, the validation focusses on the stratosphere.It is planned to extend the model to the troposphere in the future.
Both resolution and trajectory time step can be chosen freely, so the model can be used both for highly resolved process studies and long-term runs over several decades.One of the specific advantages of Lagrangian models is that they allow to study the fine filamentary structures that are usually below the resolution of conventional Eulerian models and that disappear rapidly due to numerical diffusion.This is possible since the information over the small scale behaviour of tracers is contained in the large-scale wind fields (Batchelor regime in fluid dynamics, see Haynes and Vanneste, 2004).
The main concept of the model is similar to the CLaMS 3-D model (Konopka et al., 2004), which has a more sophisticated mixing parameterization compared to models like e.g.ATILLA or Collins et al. (1997).However, ATLAS is a completely new model developed from scratch, which has no code in common with CLaMS.Some changes and improvements are introduced to the transport and mixing approach of CLaMS here, and their impact on the results is elucidated.The method of validation presented in Konopka et al. (2004) is also analyzed, and some changes are proposed for a more robust calculation of the validation parameters.Finally, the model is used to infer values for the effective horizontal and vertical diffusion coefficients for the stratosphere, and results are compared with other studies.
Section 2 contains the basic model description.Section 2.1 gives an overview of the model concept.Sections 2.2 and 2.3 describe the model implementation.While Sect.2.3 gives a detailed treatment of the mixing algorithm, Sect.2.2 describes the other relevant features of the model.The validation of the mixing algorithm is presented in Sect.3. Results for diffusion coefficients are discussed in Sect. 4. Conclusions are given in Sect. 5.

Model concept
The model design is based upon a modular approach.This allows a flexible configuration and easy exchange of components.The basic layout of the model is shown in Fig. 1.
Transport and chemistry in the model are driven by external meteorological fields.These can be provided either by a General Circulation Model (GCM module in the sketch, red box) or by existing analysis data like ECMWF ERA Interim (Simmons et al., 2006).
With this input, a large number of trajectories, each symbolizing an air parcel, is initialized and advected for some hours (filling the domain of the complete atmosphere) (trajectory module, green box).Chemistry is calculated on every trajectory like in a box model (chemistry module, yellow box).
Since there is no numerical diffusion in a Lagrangian model, some form of mixing has to be introduced to obtain realistic trace gas distributions.This has the advantage that mixing can be modeled from physical considerations and that the strength of mixing can be adjusted to observed values, while numerical diffusion cannot be controlled in strength and location.Every few hours, a mixing step is introduced (mixing module, orange box).In this step, new air parcels are added and existing air parcels are merged in regions of large flow deformation (shear and strain) to keep the density of air parcels constant.The actual mixing is accomplished by averaging the chemical species on the new or merged air parcels.As a side effect, adding and merging parcels also ensures that voids and crowded regions are avoided.After the mixing step, the model continues with the advection step again.A detailed justification of the validity of the mixing approach can be found in Sect. 2 of McKenna et al. (2002) and, shortly recapitulated, in Sect.2.3.2.

Basic model features
The initial air parcels are distributed randomly over the model domain by drawing uniformly distributed random numbers for longitude, the cosine of latitude and the vertical coordinate.The initial value of the random generator can be set to a fixed value to obtain reproducible results.Every air parcel contains variables for a configurable number of chemically active species and additionally, for a configurable number of chemically inert tracers.
The vertical coordinate is limited by an upper and lower boundary, which can be chosen freely.The number of air parcels is then chosen such that a predefined initial resolution r 0 is obtained (defined by the mean horizontal parcel distance in layers of the vertical depth z introduced in Sect.2.3.3).While it is possible to set the lower boundary to the surface, the model does only contain a basic representation of the troposphere, mainly used to provide the stratosphere with tropospheric source gases.Several components that would be necessary for a realistic representation of the troposphere are currently not implemented, like tropospheric chemistry, clouds or a treatment of sub-scale convection.
The trajectory module is based on the parallel trajectory code used in Wohltmann and Rex (2008).Time step and integration method (e.g.4th order Runge-Kutta) are freely configurable.The vertical coordinate system and the associated vertical winds can be set to four different options: -Pressure as vertical coordinate and vertical winds taken from the input files of the analysis data or GCM, which usually means winds that satisfy mass conservation by the continuity equation.
-Pressure and vertical winds from the thermodynamic equation (satisfying conservation of thermal energy) as in Wohltmann and Rex (2008).
-Potential temperature and heating rates (satisfying conservation of thermal energy).
-A hybrid coordinate that transforms from a pressure coordinate at the surface to a potential temperature coordinate in the stratosphere, as defined in Konopka et al. (2007).The hybrid winds transform from analysis winds in pressure coordinates to heating rates.
The model contains a lower and an upper boundary layer.

Mixing algorithm
The basic idea of the mixing algorithm is identical to the approach used in the CLaMS 3-D model (Konopka et al., 2004).The algorithm is as follows: -For every air parcel, determine all air parcels that are, in a certain sense (see Sect. 2.3.1 below), next neighbors of that air parcel.
-Advect air parcels with trajectory model and optionally calculate chemistry.
-If the distance between a particular air parcel and one of its former next neighbors (prior to advection) exceeds a critical distance r + after advection, insert a new air parcel in the middle between the original air parcel and the neighbor, the new mixing ratio is the mean of the mixing ratios of the air parcel and its neighbor.Do that for every neighbor.
-Determine next neighbors again.
-If the distance between a particular air parcel and several of its neighbors is lower than a critical distance r − , merge all of these neighbors and the air parcel itself to a new air parcel, the new mixing ratio is the mean of the old mixing ratios.All merged air parcels are removed immediately and can only contribute to one new air parcel.
-Advance to next mixing step.

Next neighbors
Next neighbors are determined by Delaunay triangulation: the air parcels are connected by a mesh of triangles (2-D case) or tetrahedrons (3-D case), such that the circumcircles of any triangle (or circumspheres of any tetrahedron) are empty (e.g.O'Rourke, 1998).Air parcels that are directly connected in the mesh are considered next neighbors.Depending on the model resolution, we use two different methods, which are outlined in Fig. 2: -Global layer approach (CLaMS): for high model resolutions (r 0 <150 km), the algorithm also applied in CLaMS 3-D (Konopka et al., 2004) is used.The atmosphere is divided into a set of layers.In each layer, the vertical coordinate of the air parcels is ignored and a global 2-D triangulation on the surface of a sphere is performed to obtain the next neighbors.The vertical depth z of the layers is determined with one of the methods described in Sect.2.3.3 to obtain realistic vertical diffusion coefficients and mixing (with z as the log-pressure altitude for the pressure coordinate and appropriate definitions for the potential temperature and hybrid coordinates).A staggered set of layers is used in every second mixing step to avoid a clustering of the air parcels in the vertical direction by the averaging of the vertical coordinate in new and merged air parcels (see Sect. 2.3.4).
-Local layer approach (ATLAS): for low model resolutions (r 0 ≥150 km), a new method is introduced.The main difference to the global layer approach is that we use a local layer with a limited horizontal extent for every air parcel, which is vertically centered on every single air parcel.The advantage of this approach is that it is less diffusive and gives superior results with the same mixing parameters (see Sect. 3.5).Conceptually, it is also a more realistic approximation to atmospheric mixing, since air parcels that are vertically close, but in adjacent layers do not mix in the CLaMS approach in contrast to the ATLAS approach.Unfortunately, the approach is computationally too expensive for very high resolutions.
As in the global layer approach, the air parcels tend to cluster vertically in the local layer approach.We add a random component to the vertical coordinate of newly inserted or merged air parcels to avoid the generation of artificial layers in the local layer approach.The addition of the random component is explained in more detail in Sect.2.3.4.
In the first step of the method, a global 3-D triangulation is carried out to get a preliminary set of next neighbors for every air parcel (for some details, see Appendix A).The list of next neighbors is then extended by the second neighbors (next neighbors of next neighbors) to get a cloud of air parcels around every original air parcel (the obvious approach to include all air parcels below a given distance to the original parcel in the cloud is computationally to expensive).Now, all air parcels outside an interval z− z low ,z+ z upp are deleted from the list of neighbors to center the air parcels in the local layer.The vertical mixing depth z= z low + z upp is determined as described in Sect.2.3.3.The remaining parcels and the original air parcel are now treated with a local 2-D triangulation ignoring the vertical coordinate to get a final list of next neighbors.In contrast to the global layer approach, a new 2-D triangulation is necessary for every single air parcel, but the number of points in each triangulation is much lower.
In principle, the final next neighbors could also be obtained directly from the 3-D triangulation.However, there is no constraint on the minimum vertical distance in this approach: newly inserted air parcels will tend to introduce new layers of air parcels with a lower mean vertical distance between parcels than in the last mixing step.These parcels would shield the old parcels that are a little bit further away in the 3-D triangulation and cause a finer mesh in the vertical.This Fig. 2. Sketch of the two methods for finding next neighbors.Top: in the CLaMS approach the atmosphere is divided into layers (left, blue and cyan dots denote the air parcels).In every layer, the z coordinate is ignored and a global 2-D triangulation is performed to find the next neighbors (right, view from above, in the example the orange parcels are the next neighbors of the red parcel, clipped lines symbolize the global triangulation without boundaries).Bottom: in the ATLAS approach, a global 3-D triangulation is performed (left, clipped lines again symbolize the global triangulation) and in a first step, next and second neighbors are determined for each air parcel (middle, in the example all parcels are next or second neighbors of the red parcel).All parcels outside a vertical depth z are sorted out (green parcels).The remaining blue parcels are triangulated locally in 2-D ignoring the z coordinate (right).
problem will get worse in every mixing step.In the 2-D triangulation, all these air parcels will tend to end up in the same layer and increase the density of parcels and hence, the number of parcels below the critical distance r − for merging.
Note that there is no fixed model grid at all in the ATLAS triangulation approach.There are neither distinct layers nor a horizontal grid in the model.The (virtual) model layers in the following discussion just mean all points within a certain interval [z− z low ,z+ z upp ] with z chosen freely and the vertical extent chosen to get an appropriate density of points in the horizontal.This is done in the post-processing of the model results and completely independent from the model parameterization.

Horizontal parameterization of mixing
Since the stratosphere is stably stratified, vertical turbulence is small and it is mainly the two-dimensional horizontal chaotic advection that will produce fine horizontal filaments and fractal structures in tracers on different scales.In conjunction with the vertical wind shear, sloping structures with large horizontal-to-vertical aspect ratios will be created (Haynes and Anglade, 1997).These structures will get smaller by the rapid scale collapse induced by the large-scale flow (but maintain their aspect ratio) and finally be dissipated by molecular diffusion and with the help of small-scale turbulence (e.g. by breaking gravity waves), which mixes the structures down to the molecular scale (e.g.McKenna et al., 2002;Legras et al., 2005;Balluch and Haynes, 1997;Wilson, 2004).This turbulence is intermittent in space and time due to the stable nature of the stratosphere.
The mixing concept is based on the assumption that mixing by small-scale turbulence predominantly occurs in regions of large flow deformation, e.g. in regions of large vertical shear or horizontal strain.The basic idea is to use the Lyapunov exponent to determine the critical distances r + and r − .It can be shown that large values of the Lyapunov exponent are directly related to large shear or strain rates (McKenna et al., 2002).
The Lyapunov exponent λ is a measure for the rate of separation of trajectories that are initially in close neighborship.In a simplified sense, it can be defined for trajectories by dr = λr dt (1) where r is the initial separation of the trajectories at t=0 and r+dr is the separation of the trajectories at t=dt.If the flow exceeds a critical value λ c of the Lyapunov exponent, we will insert a point.If integrating over a time t with constant λ c and initial separation r 0 , this corresponds to using a critical distance r + =r 0 exp(λ c t). λ c is a free parameter of the model and has to be optimized by comparison of model results with observations.The same is true for the length of the mixing time step t. r 0 determines the model resolution and can be chosen freely (it is identical to the r 0 used in the initialization).To relate r − to r + , we note that for incompressible flow, a small circle of radius r 0 will be deformed into an ellipse with the same area after a sufficiently small time t.If the minor and major axes are r − and r + , that means r 2 0 =r + r − and we get avoiding another free parameter in the model.In a typical shear flow, r + will apply to diffusion along the trajectory and r − to diffusion orthogonal to it.
The mixing concept ensures that mixing in the model predominantly takes place where the flow deformation is large.This is in broad agreement with observed mixing barriers like the polar vortex edge and observed well-mixed regions like the surf zone, which are indeed situated where deformation is small or large, respectively (e.g.Haynes and Shuckburgh, 2000).
λ c controls the intensity of mixing (in terms of the percentage of mixed parcels per time interval).The lower λ c is set, the more air parcels are affected by mixing.In conjunction with r 0 , which specifies the horizontal distance affected by one mixing event, λ c determines the horizontal diffusion coefficient K h (for an explicit equation, see Sect.2.3.5).Hence, λ c is a crucial parameter to adapt the strength of mixing in the model to observed values.
The parameter r 0 does not need to be the actual resolution r eff 0 of the model, since there is no constraint on the total number of parcels in the model.While the initialization starts with a mean parcel distance of r 0 , the number of parcels will change as long as the number of inserted parcels and the number of deleted parcels do not agree.Eventually, an equilibrium will be reached where this is the case.The equilibrium value for the effective resolution (defined by the mean parcel distance) is close to r 0 both for small values of λ c (high mixing intensity leads to high grid adaption frequency) and high values of λ c (no mixing cannot alter the number of parcels), but lower inbetween (see Fig. 3).
Reasonable values of t are limited to a range from approximately 6 to 24 h.For too short mixing times, the mixing step gets computationally too expensive.For too long mixing times, neighbor relationships are not conserved and the non-linearity of the flow prevents the application of the mixing algorithm.t is no crucial parameter in the parameterization, since validation results only weakly depend on the product λ c t.

Vertical diffusion coefficient
The actual magnitude of vertical diffusion enters into the model through the depth z of the layer mixed in one mixing time step, which is directly related to the vertical diffusion coefficient: the greater the diffusion coefficient, the broader the layer mixed in a given time period.Turbulent diffusion in the stratosphere is dominated by vertical diffusion processes, since the horizontal-to-vertical aspect ratio α of tracer structures in the stratosphere is large (about 250) (Haynes and Anglade, 1997).Tracer filaments will be mixed to background values in vertical direction over the whole horizontal extent of the filament quickly, while in the same time period, horizontal diffusion will only affect the outer edges of the filament.For this reason, the apparent (effective) horizontal diffusion is much larger than the actual horizontal diffusion and directly dependent on the vertical diffusion and aspect ratio by α 2 =K h /K z .The standard deviation σ of the vertical positions of air parcels that start at the same vertical coordinate z after a time t is related to the vertical diffusion coefficient by This and a normal distribution of the air parcels follows from Fick's second law with a delta function as initial condition.
In this sense, if K z = z 2 /(8 t), a layer of approximately z=2σ (containing 69% of the air parcels) is mixed after t, which gives the possibility to estimate the mixing layer depth from the diffusion coefficient.
Note that the log-pressure altitude z=H log(p 0 /p), H scale height, is used in all derivations in this section, which do directly apply to the pressure coordinate system.For the potential temperature and hybrid coordinates, z is transformed into appropriate quantities in a last step.
A first guess of the vertical diffusion coefficient is obtained from a climatology K clim z that is a function of log-pressure height (Fig. 4).We denote the diffusion coefficient as K clim z here to discriminate it from the effective vertical diffusion coefficient K eff z that is actually applied in the model (see the following paragraphs).The climatology below 50 km is taken from Massie and Hunten (1981), while the part between 50 km and 80 km is following an exponential increase (e.g.Shimazaki and Wuebbles, 1973).In the troposphere, the coefficient is set to a constant value of 10 m 2 /s.The boundary between the troposphere and the stratosphere is determined by the altitude of the tropopause.In the stratosphere below 15 km, the coefficient is set to a constant value of 0.58 m 2 /s.There is considerable uncertainty in the values of the vertical diffusion coefficient up to the present day (e.g.Wilson, 2004;Legras et al., 2005), and our parameterization, which is based on some rather old, but easily available data, fits well into the published range of values.Note also that the climatology was derived from tracer data by neglecting advective transport, which may introduce uncertainty.Since the climatology is scaled in the following, only the shape and not the absolute values are of importance anyway.
Alternatively, the vertical diffusion coefficient can be set proportional to entropy (the approach taken by CLaMS, Konopka et al., 2007) or to a constant value.Actually, the shape of the diffusion profile based on entropy is quite similar to our approach (Fig. 4, red line), which makes the choice of the method a matter of philosophy rather than being of practical implications.
The vertical mixing depth z is obtained individually for each air parcel in the ATLAS triangulation approach.Ideally, z would be obtained from where z is the log-pressure altitude of the air parcel.However, the equation actually used is α( 20km) is another free parameter of the model and defines the aspect ratio between the vertical scale z and the horizontal scale r 0 of the model in the lower stratosphere.This follows from the above equation with z=20 km: The factor 2 is due to the fact that in a layer of depth z with uniformly distributed points, the mean distance of any two points will be z/2.The reason for skipping the (physically correct) dependency on t and introducing the dependency on r 0 is that it is important to match the horizontal and the vertical scale of the model.The mean horizontal distance between air parcels in a layer of depth z must be of the order r 0 (z)=α(z) z(z)/2 to ensure that mixing events affecting a certain distance in the vertical do affect an area of the correct horizontal extent and vice versa, as given by α.A simple example may illustrate that: a typical filament in the stratosphere could have a thickness of 1 km and a horizontal extent of 250 km, corresponding to an aspect ratio of 250.If the vertical diffusion coefficient would mix a layer of 1 km in 12 h to background values, it would also mix the structure of 250 km extent horizontally to background values in that time.A filament of thickness 2 km needs to have a width of 500 km and will need correspondingly longer to be mixed.If the horizontal and vertical resolution of the model would be decoupled, the mixing in the model would not work anymore in a realistic way.
If Eq. ( 4) would be used for the vertical resolution and the aspect ratio would still be respected, the horizontal resolution would be no free parameter and e.g.r 0 (20 km)=α(20 km) z(20 km)/2.With z as small as inferred from the vertical diffusion coefficient climatology (e.g.z≈450 m and r 0 ≈55 km for K z =0.58 m 2 /s and t=12 h), the number of air parcels and model layers would increase to values computationally too expensive.That is, the method used in the model overestimates the effective vertical diffusion coefficient due to computational constraints.Note that CLaMS effectively implements the layer thickness in an equivalent way and suffers from the same problem.This can be partly compensated by setting λ c to high values, which reduces the number of mixing events.However, while this will lower the effective bulk diffusivity of the model, every single mixing event will still take place over a vertical distance that is too large.
Finally, z is split into an upper part z upp , which we use to search for all points between z and z+ z upp and a lower part z low to search for points between z and z− z low .The simplest approach would be to set both values to z/2.In our approach, z low and z upp can take different values.First, if the interval [z− z(z)/2,z] crosses the tropopause, z low is set such that the interval ends at the tropopause.The same is done with z upp and the interval [z,z+ z(z)/2].To avoid problems in other regions where K z changes rapidly with altitude, z low is calculated as half of the mean value of z(z ) between z and z− z(z)/2, An analoguous aproach is taken with z upp .
www.geosci-model-dev.net/2/153/2009/Geosci.Model Dev., 2, 153-173, 2009 Equation ( 5) can be interpreted as introducing an effective α(z) by z(z)=2r 0 /α(z) and This implies a higher effective aspect ratio at altitudes with a lower z and vice versa.Very probably, this ad hoc assumption of the variation of α(z) with altitude is not compatible with the real atmospheric variation, but this variation follows from the horizontal resolution of the model held constant with altitude.In future versions, that may be changed to allow a varying horizontal resolution.For the following validation the effect is small, since the model is compared to measurements around 20 km altitude.

Inserting new parcels
If a new air parcel is inserted, its coordinates are calculated as follows: the longitude and latitude of the new air parcel are calculated by transforming all longitudes and latitudes of the old air parcels to cartesian coordinates with a local orthographic projection, and averaging the x and y coordinates.The result is transformed back.
In the CLaMS approach, the vertical coordinate is the average of the old vertical coordinates.This would lead to a clustering of the air parcels in the middle of each layer with fixed layers: the average of some random numbers from an interval will typically be close to the middle of the interval and is bounded by the original numbers.This increases the number of air parcels near the middle of the layer in every mixing step and leads to a positive feedback which results in the generation of layers in the vertical distribution of parcels.For this reason, a staggered set of vertical layers is used in every second mixing step in the CLaMS global layer approach.However, this approach still cannot avoid a certain clustering of the parcels into layers (Fig. 5).
Somewhat surprisingly, there is the same problem in the ATLAS local layer approach of triangulation, if the CLaMS method of averaging is used.Although there are no predefined layers, the same process can be triggered once there is a small random initial deviation from a uniform vertical distribution of the air parcels: the deviation will tend to be amplified in the next mixing steps by the feedback described above.This results in a process of self-organisation of the vertical coordinate values, where neighboring layers of enhanced parcel density are in concurrence for the parcels inbetween.
It is necessary to work against this effect to avoid the formation of pronounced vertical layers and associated problems with vertical mixing.Since there is no possibility of a staggered set of layers as in the CLaMS approach, a random number is added to the average of the vertical coordinates in every newly inserted or merged air parcel to achieve a more uniform vertical distribution of air parcels (Fig. 5).The random number is chosen from an interval − z low , z upp , with these values again calculated as in Sect.2.3.3.In most cases, the addition of the random component only slightly increases the diffusivity.
Since there is still a certain vertical clustering of the air parcels in the CLaMS global layer triangulation approach with simple vertical averaging, we also add a random number to the vertical coordinate in the CLaMS triangulation approach by default.This modified CLaMS approach is the only computationally feasible method for higher resolutions, if we want to achieve a distribution as uniform as possible also here.A discussion of the differences of the approaches in the results is given in Sect.3.5.

The effective large-scale diffusion in the model
It is possible to estimate the diffusion coefficient from the mean distance r of any two points in a Gaussian distribution of air parcels as the one used for Eq. ( 3).It turns out that r 2 =2σ 2 and hence the estimate is According to this equation, we estimate the diffusion coefficent in the model for every newly inserted or merged parcel.For comparison with other published values of the diffusion coefficient, the diffusion coefficient averaged over a certain geographical area and not for an individual air parcel is needed.Hence, a mean diffusion coefficient is calculated for every mixing step by averaging over the diffusion coefficients of the air parcels in a given area.Parcels not affected by mixing are assigned a coefficient of zero.It is possible to derive equations that provide some insight on how the diffusion behaves if we change the free model parameters r 0 , λ c , α and t, see also Fig. 5 of Konopka et al. (2003).If two air parcels are mixed in the r + or r − step, the horizontal diffusion that the newly inserted air parcel experiences can be roughly estimated to and the vertical diffusion to since the typical distances between the air parcels that contribute to the new parcel are r ± and z/2.It is obvious that the diffusivity of the model increases for decreasing resolution, just because every single mixing event occurs over a larger distance.For decreasing λ c (more mixing), K eff h+ (diffusion in flow direction) decreases, while K eff h− (diffusion orthogonal to flow direction) increases.For the spatially averaged diffusion and also for the average diffusion that a single parcel experiences over a large number of mixing steps, the percentages of added and merged parcels also matters (see Fig. 11).These are also a function of the mixing parameters, particularly of λ c , where lower values correspond to higher percentages of mixed parcels and a higher average diffusion coefficient (and vice versa).Supplementary to the analysis in Konopka et al. (2003), we find that the increase of diffusivity with lower λ c in CLaMS and ATLAS is not only due to the increase in K eff h− , but that the higher percentage of mixed parcels contributes significantly.Note that the vertical and horizontal diffusion can be varied independently from each other within certain limits, since λ c appears only in the horizontal diffusion coefficients and α only in the vertical.

Additional patchyness from noise in the wind fields
In principle the strength of mixing between air parcels can be tuned to the real atmospheric mixing.But in practice random motions from noise in the driving wind fields and their Eulerian discretization result in spurious random transport of air masses.This contributes to unavoidable excessive dispersion even in a purely Lagrangian transport approach, which generates spurious transport across mixing barriers and steep tracer gradients.This results in the artificial generation of small scale structure and patchyness in the tracer fields.The mixing in the model has to be tuned to keep this patchyness at realistic levels.Hence, the mixing in the Lagrangian model will always tend to overestimate the real diffusion in the stratosphere somewhat, though to a much lesser degree than an Eulerian model.Eulerian models suffer from noise and discretization in the driving wind fields in a similar way, but additional numerical diffusion largely dominates these effects such that they do not become visible.

Introduction
The validation of ATLAS is based on the same aircraft data set that is used in the validation of CLaMS 3-D in Konopka et al. (2004)  We will first describe the setup of the model runs and show some examples for model results.A quantitative validation will then be done by defining parameters for the agreement of the modeled and observed tracer-tracer relationships and the agreement of the spatial roughness of model and observations.

Setup
The model runs use the hybrid coordinate and are driven by meteorological data from the ECMWF ERA Interim reanalysis (Simmons et al., 2006(Simmons et al., , 2007) ) on 60 model levels (6 h temporal resolution, 2 • ×2 • horizontal resolution).Heating rates (clear sky) are also taken from ERA Interim.ERA Interim data in the stratosphere has improved considerably compared to the ERA-40 reanalysis (Simmons et al., 2006(Simmons et al., , 2007)), and, without any modifications, produces good agreement with observations in the following.Vertical diffusion coefficients are taken from the built-in climatology.The trajectory model uses a 4th order Runge-Kutta method for integration and a time step of 30 min.A change of the time step to 10 min does not lead to any significant improvements in the results.
Model runs are performed for all combinations of the resolutions r 0 =50, 100, 150, 200 and 300 km and critical Lyapunov exponents λ c =1, 2, 2.5, 3, 3.5, 4, 4.5, 5, and 10 d −1 .The number of air parcels ranges between 10 4 (300 km) and 10 7 (50 km).All runs at and below 150 km resolution are started both with the CLaMS and the ATLAS method of finding the next neighbors.For higher resolutions, only the CLaMS method can be used.t is set to 12 h and α to 250.Additionally, runs with fixed r 0 (100 km) and λ c (3 d −1 ) are started for different values of t or α.The effect of switching the random vertical coordinate on or off is also studied.
The lower model boundary is at the hybrid level 350 K and the upper boundary is around 1500 K.The number of virtual levels of the model (in the sense of the number of mixing depths z one needs to stack on top of each other to span the model domain) is 8 for 300 km resolution, 12 for 200 km, 16 for 150 km, 24 for 100 km, and 48 for 50 km.
The model is started on 1 November 1999, while the tracers are initialized on 1 December.This is done to avoid any spin-up effects in the mixing: usually, the equilibrium number of air parcels is lower than the number of air parcels the model is initialized with.This leads to a spin-up phase where more points are deleted in every model cycle than inserted, which in turn alters the mixing properties of the air parcels in an unrealistic way.
Outside the polar vortex, the methane tracer is initialized with monthly mean HALOE data from November and December 1999 (weighted equally) as a function of pressure and equivalent latitude ϕ E (Grooß and Russell III, 2005).Inside the lower vortex (ϕ E >64 • N, θ <712 K), an OMS balloon profile from the LACE instrument measured at 19 November is used (Ray et al., 2002), since HALOE is not able to measure in high latitudes and darkness.The profile is corrected for the descent between 19 November and 1 December, as inferred from a model run driven by the heating rates described above.The profile above 712 K in the vortex is a profile from the highest equivalent latitude bin of HALOE scaled to match the LACE values below.The Halon-1211 tracer is initialized by a methane to Halon-1211 tracer-tracer relationship measured at the OMS balloon flight.
The model results for the methane tracer are compared to measurements of the ACATS-IV (Romashkin et al., 2001), ALIAS (Webster et al., 1994) and Argus (Loewenstein et al., 2002) instruments on 11 flights of the ER-2 aircraft between 20 January and 12 March (Newman et al., 2002).The Halon-1211 model results are compared to ACATS measurements on the same flights.Flights started from Kiruna (Sweden) and probed the stratospheric polar vortex and surf zone.Halon and methane measurements from ACATS are available every 70 s, corresponding to a resolution of about 20 km, while ALIAS and Argus measurements are available every 2 s (about 1 km resolution).Wherever possible, we use the high resolution data, tracer-tracer relationships and Halon plots use the low resolution.

Examples
Figure 6 shows model results and measurements for the 27 January flight as an example.Results are shown for 50 km resolution and optimal mixing parameters (λ c =3 d −1 , α=250, t=12 h, see Sect.3.5).Only the highest tested resolution of 50 km is able resolve the fine filamentary structures at the edge of the polar vortex.The flight crosses the vortex boundary twice at a mean flight level of 50-60 hPa and contains a dive to 100 hPa around 12:00 UTC.For the comparison, the flight path is transformed to the position of the probed air parcels at 12:00 UTC by calculating short forward or backward trajectories starting at every measurement.Model results at 12:00 UTC are interpolated to the transformed measurement positions by averaging the mixing ratios over the next model neighbors of every measurement location, with a 2-D triangulation performed in the same way as in the mixing step of the model.The model values in the average are weighted according to their distance to the measurement: the weight is 1 for distances smaller than 2r 0 and decreases linearly to zero at 3. quite good, both for the gradient across the vortex edge and the vertical structure.Inner and outer vortex air masses can clearly be distinguished both in the model and the observations below 70 hPa and the gradual transition from inner to outer vortex values above 70 hPa is also matched quite well by the model.The modeling of Halon is challenging, since the mixing ratios of Halon very sharply drop to zero above 50 hPa.Hence, the location of the drop is a good test for the accuracy of the heating rates, the accuracy of the winds and the mixing parameterization.Panel e shows the modeled and observed tracer-tracer relationships, which we will use later to validate the model.Panel f shows the methane tracer field at the flight level, with the polar vortex characterized by low methane values.
The interpolation to the probed air parcels of the flight will inevitably lead to additional diffusion in the model results, and different approaches for finding corresponding model values show considerable differences in the results, see also Konopka et al. (2003).Figure 7 shows the remarkable differences for the flight on 11 March, which crossed a filament  of extra-vortex air inside the vortex twice during the flight.While the filament is clearly visible if using only the nearest model neighbor to each measurement as the corresponding model value (free of any interpolation), it is more difficult to see for the average over the next neighbors and barely visible for the average over the next and second neighbors.On the other hand, parts of the flight where gradients are low show a better agreement with more averaging in the interpolation.This effect gets more prominent if λ c is further increased (not shown).In this sense, the interpolation must be considered as an integral part of the modeling.Hence, for the actual validation, we try to find measures free of any additional interpolation.
Figure 8 demonstrates how the results change if the mixing parameter λ c is varied.The left column shows an example where the mixing in the model is stronger than observed, leading to a mixing curve in the model that is on the concave side of the observed mixing curve.The middle panel shows an example for optimal mixing parameters and the right panel shows an example for too low mixing, where the modeled mixing curve is on the convex side of the observed curve.

Validation measures
The quality of the model runs is judged by calculating a parameter ε for the agreement of the modeled and the observed tracer-tracer relationships for every run.Likewise, a parameter γ for the agreement of the spatial "roughness" of the observed and modeled tracer fields is calculated for every run.
While the parameter ε is a straightforward measure for the mixing intensity of the model, the parameter γ needs some explanation.For low mixing intensities (high λ c ), the model fields show scatter and filaments, where the measurements do not show any.There are two possible reasons for this: first, tracer structures that would be mixed to background values in reality are not mixed in the model due to the low mixing intensity.But additionally, the spurious transport described in Sect.2.3.6 generates small scale structure and leads to the same effect, see also Legras et al. (2005).The two parameters ε and γ give some possibility to distinguish between the two effects: spurious transport generates small scale structure and hence has an effect on γ .But these transport effects do not result in changes in the tracer relationships that are the  basis for ε in first approximation, since the artificial filaments will be mixed with the same intensity as any filament.Any exaggeration of model mixing used to balance the generation of patchyness and to tune γ to values close to observations results in deviations of modeled and observed tracer relationships which show up in ε.
In the following, we will denote the ith methane measurement with m i (counting over all flights without resetting the counter), and the corresponding Halon-1211 measurement with h i .
The parameter ε for the comparison of the mixing intensity is defined by the sum of the distances in tracer-tracer space of the measurement pairs Konopka et al. (2004).That is, A sketch of the method is shown in Fig. 9 (left).Since the absolute magnitude of the sum for methane and the sum for Halon-1211 is very different, relative differences are taken in Konopka et al. (2004) ε where N is the total number of measurements over all flights.Small values mean a good agreement between the mixing intensity of the model and the observations.Corresponding model pairs (m model i ,h model i ) are found by looking for the nearest model neighbor to the measurement location in Konopka et al. (2004) to avoid interpolation to the measurement location.Interpolation would cause additional mixing and may considerably alter the tracer-tracer correlation curve (e.g. by introducing spurious mixing lines, not shown).
There are some possibilities to improve this definition: first, for mixing ratios near zero (as it is the case for Halon-1211 in several parts of the flights), the relative differences go to infinity, and some single measurements can errorneously dominate the overall sum due to measurement error.
Secondly, errors in the wind fields used to drive the model can cause the position of some filament in the model to be shifted from its position in the observations to a slightly different place.This would cause a large distance of the corresponding (m i ,h i ) pairs in tracer-tracer space (where one pair is in the filament and the other is not).However, the mixing intensity could still be right, that is the modeled mixing ratios in the filament would match the observed tracer mixing ratios and the filament would have the right extent and shape.Hence, in a model with perfect mixing but slight systematic errors in the driving wind fields, the mixing curve of the model would still look exactly like the observed mixing curve.Since ε should only measure the mismatch in mixing intensity and not displacements of the filaments, this is not optimal.Konopka et al. (2004) argues that for too low mixing intensities, spurious filaments that would be mixed in the real atmosphere would also cause large distances of corresponding points in tracer-tracer space, which outweigh the effect of displaced filaments in ε.We argue that the key measure of mixing properties of the model is the shape of the tracer-tracer relation and not the point-to-point distances between model and observations in tracer-tracer space and use an approach that results in a clearer separation of the roughness effect (which is not only caused by too low mixing, but also by the random errors in wind fields) and the effect of displaced filaments from the mismatch in mixing intensity in tracer-tracer space.Hence, we try to put the roughness effect only into γ and try to define ε such that it is only affected by the mixing intensity.
A suitable definition minimizes if the shape of the mixing curve of the model gets more and more similar to the shape of the observed mixing curve.We also look for a definition that allows us to take the model values directly to avoid interpolation.This can be accomplished by the absolute differences of the modeled Halon-1211 values to a curve fitted through the observed Halon-1211 as a function of observed methane.A sketch of the method is shown in Fig. 9 (right).The sum of these differences is largely proportional to the area between the mixing curves of the model and the observations.The observed methane is taken from the average of the three instruments.To obtain the absolute difference, we evaluate the fit function at the modeled methane value.That is where k(j ) is the flight number belonging to the modeled value j and ĥobs k is a polynomial fit through the observed tracer-tracer relationship of flight k.Since we use a continuous curve for the observed tracer-tracer relationship, we do not need to use a single corresponding value m model i for every measurement m observed i .Instead, we use all model air parcels that are next neighbors of the measurement locations, from the 2-D triangulation used in the mixing step and with z centered at the measurement location.This method gives us a cloud of air parcels along the flight path.This is indicated here by using j for i as the index and M for N for the number of points.A potential disadvantage of the method is that it can only be applied to data only containing one compact mixing curve (and no additional mixing lines).
In some cases individual methane values in the model are outside of the range of the measured methane values, where the fitted curve is valid.In this case, we do not use these values, so that ε is not based on exactly the same set of points for some model runs.The validity of this approach is verified by extending the mixing curves with a best guess of the likely missing methane values and calculating the parameter again (these tests were not used for the validation).Results are almost identical.The best guess is simply obtained by extending the data of the tracer curves with methane and Halon-1211 data from the third flight for the 5 other flights at the end of January or start of February, since it has the largest range of methane values.The flights at the end of February and in March are extended with data from the eighth flight.
The parameter γ for the roughness of the field in Konopka et al. ( 2004) is defined in two steps: the sum of the absolute differences of successive methane measurements is The same can be done for the model values of methane.The ratio of the measured and modeled is Values near 1 are defined as the optimum.With this definition the magnitude of the roughness will depend on the spatial scale of the structures we are looking at.It will make a difference if the measurement locations are 1 km or 20 km apart.More importantly, it will also make a difference if the model air parcels are 50 km or 300 km apart.In this sense, γ is not a very robust measure of the roughness.Hence, we use a definition of γ that compares the roughness observed and model at the spatial scale of the model resolution and is robust to changes in the spatial distance of the measurements.
Our approach compares the standard deviations of the modeled and observed methane values at the spatial scale of the effective model resolution: we determine the next model neighbors of each measurement i (exactly as described for the other parameter ε) and calculate the standard deviation σ i of the modeled methane values for each set of next neighbors (avoiding any interpolation).Now, all standard deviations are averaged to give a mean standard deviation σ model .In addition, we calculate the mean distance of the points in every set of next neighbors and compute its average d.Now, we look for all pairs (m observed1 i ,m observed2 i ) of measurements which are separated by a distance within [ d− d, d+ d], with d=5 km.It can be shown that the standard deviation of this set of points is given by the following equation: where N is the number of pairs.Now, we define in equivalence to Eq. ( 14).Theoretically, a term for the statistical measurement error in the observed values needs to be considered, since the standard deviation of the observations will not only be caused by real variations of the mixing ratio but also by measurement noise.This would cause an offset σ bias to σ observed and the denominator would be σ 2 observed − σ 2 bias .However, it turns out that this correction is negligible.
It is important to note that our γ gives us only a measure for the deviation of the model variability from what we would expect from the observations at the given model resolution.There is no penalty if small filaments just do not appear in coarser model resolutions.In this sense it gives a relative but not absolute measure for the expected variability.
Note that our definition of ε and γ does not imply a perfect agreement of the observed and modeled time series along the flight path if ε = 0 and γ =1.For example, all filaments could be displaced by a small distance between the modeled and observed time series along the flight path in this case.

Validation results
Figure 10 shows the results for the parameters γ and ε for different model resolutions r 0 and different values of the critical Lyapunov exponent λ c (for a comparison of our definition of ε and γ and the one given in Konopka et al., 2004, see Appendix B).
For a given r 0 , data points show more intense mixing to the left (low values of λ c ) and less intense mixing to the right (high values of λ c ). λ c =10 d −1 is virtually identical with no mixing at all (λ c →∞).
The agreement of the modeled tracer-tracer relationship with the observed one generally improves with increasing λ c .For higher resolutions it does only reach a shallow minimum at roughly λ c =3 d −1 to λ c =4 d −1 (see also Fig. 12), which compares well with the value λ c =3 d −1 given in Konopka et al. (2004).This basically means that the observed mixing is extremely low and very anisotropic: the ratio between r + and r − is larger than 20 for these values of λ c .Switching mixing off in the model is actually not a bad option for optimizing ε.Indeed, the observed tracer-tracer relationships on the first flight on 20 January and the last flight on 12 March are very similar, indicating that mixing in the real stratosphere is extremely low.
Generally, the agreement of the tracer-tracer relationships improves with higher resolution r 0 if λ c is held constant.This is due to the fact that the model gets more diffusive with increasing distance between the mixed air parcels.
Where the triangulation approach of the ATLAS model is applicable (dashed lines), it is generally less diffusive than the CLaMS approach (solid lines), see also Fig. 14.For the same set of mixing parameters, the results of the ATLAS triangulation approach are superior to the CLaMS approach.One likely reason for that is that the number of mixed parcels is lower in the ATLAS approach than in the CLaMS approach for the same λ c (Fig. 11).This is due to the lower average vertical distance between mixed air parcels: the maximum distance of two mixed air parcels is z in the CLaMS approach (if both are at the edges of the layer), but only z/2 in the ATLAS approach (since one of the parcels is by definition at the center of the local layer).A possibility to work against that in the CLaMS approach would be to increase λ c until the percentage of mixed parcels is low enough.However, there would still be qualitative differences between the approaches.E.g., if an air parcel is situated at the boundary of a layer in the CLaMS approach, it is possible that it does not mix with a vertically close parcel just because it is in the adjacent layer.In this respect, the ATLAS approach behaves qualitatively more like in the real atmosphere.In fact, there are some runs with the ATLAS approach for 150 km resolution, which show a lower ε than any run with the CLaMS approach in the same resolution, no matter what λ c is used.
However, note that for coarse resolutions (say 300 km) the mixing process in the ATLAS as well as the CLaMS approach will become unrealistic.The diffusion coefficients of single mixing events will become unrealistically large and that can only be compensated by setting λ c to large values to obtain realistic bulk diffusivities (see Sect. 4).But this will still be a better option for long-term runs in coarse resolutions than using an Eulerian approach.For all resolutions, the roughness ratio γ increases more or less uniformly with decreasing λ c (the model fields are too smooth for low values of λ c and too rough for high values of λ c ), see Fig. 10 (right).The optimal value for λ c based only on γ ≈1 is about 2 d −1 to 3.5 d −1 for all resolutions and both triangulation approaches.
Since the effect of random errors in the wind fields and numerical discretization shows up mainly in the roughness ratio γ and not in the mixing mismatch ε (because the artificial filaments would still be mixed with the right intensity), γ should optimize at lower values of λ c (at higher mixing intensity) than ε.CLaMS triangulation @ min(ε) CLaMS triangulation @ min(γ) ATLAS triangulation @ min(ε) ATLAS triangulation @ min(γ) Fig. 12. Critical Lyapunov exponents where the mixing mismatch ε minimizes (blue) and the roughness ratio γ is closest to one (cyan).Solid lines show the CLaMS triangulation approach, dashed lines show the ATLAS triangulation approach.
(cyan).Note however that the minima are shallow and that some uncertainty is introduced by this.In addition, the AT-LAS triangulation approach is optimal at smaller values of λ c than the CLaMS approach, possibly since it is less diffusive and needs more mixing in the model to obtain optimal results.The values at 300 km should not be taken too seriously, since the minimum of ε is at the highest tested λ c of 10 d −1 here.Often, it will be unavoidable to find a compromise in optimizing γ and ε, i.e. to fight the patchyness by setting the mixing intensity to slightly higher values than suggested by ε (lower values for λ c ).
The results are not very sensitive to the value of the aspect ratio α.Values between 200 and 500 give comparable results to 250.As long as the product λ c t stays constant, the re-sults are also not very sensitive to the mixing time step t.Values of 6 h and 24 h give comparable results to 12 h.These results are also confirmed by the validation of the CLaMS model (Konopka et al., 2004(Konopka et al., , 2005)).
Figure 13 shows that the effect of adding the random vertical coordinate is small in the range of the minimum ε values both for the ATLAS and CLaMS triangulation.However, switching the random vertical coordinate off does decrease the diffusivity for the lowest λ c values in the ATLAS triangulation approach.This is very likely due to the vertical clustering into layers of enhanced parcel density that sets in if the random coordinate is switched off and the mixing intensity is high.This effectively suppresses vertical mixing and lowers the effective diffusion coefficient to values which are normally observed at higher λ c values.That is, while improving the vertical homogeneity of the model, the random coordinate has either negligible effects on the diffusivity (high λ c ) or is just necessary for a realistic behaviour of vertical mixing in the ATLAS approach (low λ c ).

Estimation of diffusion coefficients
It is desirable that the diffusion coefficients obtained with the optimal mixing parameters are in the range of other estimates.However, there are few direct and independent estimates of diffusion coefficients (e.g. from radar measurements, like Fukao et al., 1994, see also Wilson, 2004).Most estimates (including ours) are based on indirect methods as the reconstruction of observed tracer filaments, profiles or tracer-tracer relationships, usually with a diffusion component and sometimes a Langrangian transport component, e.g. by a simple 1-D diffusion model for the mean vertical profile (Massie and Hunten, 1981), numerical models for strain and diffusion of filaments (Balluch and Haynes, 1997;Waugh et al., 1997), or by averaging stochastic back trajectories (Legras et al., 2005).In contrast to e.g.Balluch and Haynes (1997), the information about the diffusion coefficient comes mainly from the tracer-tracer relationship in our method and not from reconstructing fine filaments, which are only resolved in the highest resolutions.A general problem of a comparison are the large differences in the methods.Some are local (radar on the 100 m scale), some give bulk estimates, some are Eulerian, some Lagrangian, some direct and some indirect.
Figure 14 shows the effective mean (bulk) vertical and horizontal diffusion coefficients north of 60 • N at 430 K as a function of λ c and r 0 , calculated as in Sect.2.3.5.Vertical diffusion coefficients obtained from the studies mentioned above are also shown for comparison.The modeled vertical diffusion coefficient increases with lower λ c , which is due to the increasing number of mixing events, and increases with coarser resolution, which is due to increasing distances between mixed parcels.Figure 15 shows the value of the vertical diffusion coefficient at the optimal values of ε and γ in the same manner as in Fig. 12.The vertical diffusion coefficient takes values between 0.07 m 2 /s and 0.36 m 2 /s for the optimal ε and the CLaMS triangulation approach, and values between 0.02 m 2 /s and 0.1 m 2 /s for the optimal ε and the ATLAS triangulation approach.The optimal γ shows vertical diffusion coefficients that are probably too large, if compared with the estimates from the other studies.These values are increased due to errors in the wind fields.In summary, our results suggest a value of about 0.1 m 2 /s for the vertical diffusivity.As noted, the estimated bulk diffusivities will probably be more realistic for high resolutions, since the diffusivities of single mixing events are probably too large for coarser resolutions.
Our estimate agrees well with the value 0.1 m 2 /s given for the lower stratospheric surf zone in Legras et al. (2005).Values from the radar measurements (0.1 m 2 /s to 0.5 m 2 /s for Fukao et al., 1994) tend to be larger than other estimates.Also the Massie and Hunten (1981) value of 0.58 m 2 /s at 22 km seems to be relatively large.As estimated in Sect.2.3.3,values like these would roughly correspond to a minimum width of a filament of 50 km before it disappears.One reason might be that the effects of the large-scale advection by the Brewer-Dobson circulation are not considered in studies like Massie and Hunten (1981).On the other hand, the estimate of Balluch and Haynes (1997) (0.01 m 2 /s to 0.001 m 2 /s) and of Legras et al. (2005) for the inner vortex (0.01 m 2 /s) are on the lower side of our estimate.Differences can at least partly be explained by the likely dependence on season and location (e.g.single filaments in Balluch and Haynes, 1997, or the average north of 60 • N in our approach compared to the separate values for surf zone and vortex in Legras et al., 2005  for the triangulation approach of ATLAS.Shaded areas and solid black lines show vertical diffusion coefficients from other studies for comparison: Massie and Hunten (1981), Fukao et al. (1994), Balluch and Haynes (1997), Legras et al. (2005).Values for λ c =10 d −1 and r 0 =300 km are zero and not shown.Right: same for the effective horizontal diffusion coefficient.Grey lines are the vertical diffusion coefficients from the left panel scaled by α 2 .Solid black line shows estimate from Waugh et al. (1997).CLaMS triangulation @ min(ε) CLaMS triangulation @ min(γ) ATLAS triangulation @ min(ε) ATLAS triangulation @ min(γ) Fig. 15.Value of the vertical diffusion coefficient where the mixing mismatch ε minimizes (blue) and the roughness ratio γ is closest to one (cyan).Solid lines show the CLaMS triangulation approach, dashed lines the ATLAS triangulation approach.Shaded areas and solid black lines show vertical diffusion coefficients from other studies for comparison as in Fig. 14.
The right panel of Fig. 14 shows the horizontal diffusion coefficient, which ideally should be a scaled version of K z , since Haynes and Anglade (1997) claims K h =α 2 K z .The grey lines show the vertical diffusion coefficients from the left panel scaled by α 2 .In general, the agreement is quite good for values smaller than λ c =5 d −1 , although the horizontal diffusion coefficients are somewhat larger than expected from K z (note the logarithmic scale).This could be interpreted as the effective aspect ratio of the model being somewhat larger than 250.The modeled horizontal diffusion coefficient also increases with lower λ c , which is now not only due to the increasing number of mixing events, but also modified by the increasing diffusivity orthogonal to the flow direction and the decreasing diffusivity in flow direction (Sect.2.3.5).Waugh et al. (1997) gives a value of 10 3 m 2 /s for the horizontal diffusion coefficient, which seems to be somewhat lower than our effective diffusion coefficients.The sharp drop in the horizontal diffusion at λ c values of 5 d −1 and 10 d −1 is due to the fact that only very few mixing events occured at these values and that all events were merging events (r − ), while at lower λ c values, also r + events occured, which have a much larger diffusion coefficient (at the same mixing parameters).That hints at the fact that the mixing in the model will get unrealistic at these large values of λ c .

Conclusions
The new global Chemical Transport Model (CTM) ATLAS with Lagrangian transport and mixing was presented.In this paper, we focussed on a basic model description and the presentation and validation of the transport and mixing module of the model.In conjunction with a projected second part concentrating on the chemistry module, the paper is thought to serve as the reference citation for the model.
The ability of the transport and mixing model to reproduce observed tracer data and fine-scale tracer structure has been successfully demonstrated.Excellent agreement to observed structures up to scales of typical polar vortex filaments can be reached with suitable settings of the free parameters of the model, which were extensively tested and validated.
The results of the validation suggest a vertical diffusion coefficient on the order of 0.1 m 2 /s in the high-latitude lower stratosphere, which fits well into other estimates.However, the vertical diffusion coefficient remains a parameter that can only be estimated with large uncertainty, ranging from 0.5 m 2 /s to 0.001 m 2 /s in the literature cited here.
The model makes heavy use of concepts developed for the CLaMS model, and consequently part of the study focussed on comparison with that model.Generally, results seem to be of comparable quality and the changes introduced in AT-LAS were demonstrated to improve the results (new concept for identifying the neighbors for the mixing process) or to have little effect (effect of random vertical coordinate in the CLaMS triangulation, diffusion coefficient climatology versus entropy-based diffusion coefficient).
As always, time and space did not permit to study all aspects of the model and the validation.For example, since the diffusivity of the atmosphere will vary in time and space, it would be interesting to validate the model with data from other geographical regions or seasons, e.g. to see if the model is also able to reproduce mixing lines in cases of strong mixing events (i.e. in the case of fast anomalous mixing, compared to the constant mixing which is slowly lifting the correlation curve on the concave side).It would also be illuminative to compare the validation results of the Lagrangian model to the results an Eulerian transport code would give.
Future plans with the model include the improvement of the representation of the troposphere, and in the long run, the full coupling to a Atmosphere Ocean General Circulation Model.11) (originally from Konopka et al., 2004) with nearest neighbors as corresponding model points, for comparison with ε as defined in Eq. ( 12) shown in Fig. 10 (left).resolution) and the grey lines show the results if the high resolution measurements of Argus and ALIAS are used (about 1 km resolution).While the γ obtained by Eq. ( 16) remains virtually unchanged by resolution changes, the other methods are fairly sensitive to them.

Fig. 3 .
Fig. 3. Effective resolution for different values of the critical Lyapunov exponent λ c (axis) and resolution r 0 (colors).

Fig. 5 .
Fig.5.Vertical distribution of air parcels for the ATLAS triangulation approach (top) or the CLaMS triangulation approach (bottom) and the random vertical coordinate off (left) or on (right).Results of a model run with r 0 =150 km (top) or r 0 =100 km (bottom), λ c =2 d −1 after 4.5 months.The scatter plot shows the distribution of air parcels in the potential temperature-latitude plane (ignoring longitude).Only a subset of all air parcels is shown to make the clustering visible.The blue line shows the vertical density of the points (as the number of points per the local value of the layer depth z).Red lines are the layers of the CLaMS approach.

Fig. 6 .
Figure6shows model results and measurements for the 27 January flight as an example.Results are shown for 50 km resolution and optimal mixing parameters (λ c =3 d −1 , α=250, t=12 h, see Sect.3.5).Only the highest tested resolution of 50 km is able resolve the fine filamentary structures at the edge of the polar vortex.The flight crosses the vortex boundary twice at a mean flight level of 50-60 hPa and contains a dive to 100 hPa around 12:00 UTC.For the comparison, the flight path is transformed to the position of the probed air parcels at 12:00 UTC by calculating short forward or backward trajectories starting at every measurement.Model results at 12:00 UTC are interpolated to the transformed measurement positions by averaging the mixing ratios over the next model neighbors of every measurement location, with a 2-D triangulation performed in the same way as in the mixing step of the model.The model values in the average are weighted according to their distance to the measurement: the weight is 1 for distances smaller than 2r 0 and decreases linearly to zero at 3.5r 0 .Panels a and b show measurements and model results for methane and Halon as a function of time, while Panels c and d show the same as a function of pressure.The agreement of the modeled and observed methane and Halon-1211 values is

Fig. 7 .
Fig. 7. Different methods of the interpolation to the flight path.Measured methane mixing ratios (colored dots) and modeled methane mixing ratios (black line) as a function of flight time (UTC) for the 11 March flight (λ c =3 d −1 , r 0 =50 km): for the nearest model neighbor of the measurement (upper left), the average of the first neighbors (upper right) and the average of the first and second neighbors (lower left).More averaging means a better agreement in areas of low gradients (due to reduced scatter), while less averaging better preserves fine filament structures.Lower right: modeled tracer field for the Northern Hemisphere for all air parcels (dots) between 432-444 K.

Fig. 9 .
Fig. 9. Sketch of the two methods for calculating the parameter ε.Left: in the definition of Konopka et al. (2004), ε is the sum of the distances of corresponding individual points in the Halon-methane tracer-tracer space.Right: in this work ε is the sum of the distances of modeled values of Halon-1211 to a curve fitted through the observed Halon-1211 values as a function of observed methane.

Fig. 10 .Fig. 11 .
Fig.10.Mixing mismatch and tracer roughness.Left: the parameter ε for the difference of the modeled and observed tracer-tracer relationship for different values of the critical Lyapunov exponent λ c (axis) and resolution r 0 (colors).Solid lines show results for the triangulation approach of CLaMS, dashed lines show results for the triangulation approach of ATLAS.Right: same for the parameter γ for the tracer roughness ratio.
Figure 12 shows that this is in fact what is observed in most validation runs.It shows the λ c values where the ε curves minimize (blue) and γ is closest to one

Fig. 13 .
Fig.13.Effect of adding a random number to the vertical coordinate of newly inserted parcels.The solid blue line shows ε for 100 km resolution, the CLaMS triangulation and the random vertical coordinate on, the dashed red line for 150 km, the ATLAS triangulation and the random vertical cordinate on.Dotted lines are the same with the random vertical coordinate off.

Fig. 14 .
Fig.14.Left: effective vertical diffusion coefficient north of 60 • N at a layer around 430 K for different values of the critical Lyapunov exponent λ c (axis) and resolution r 0 (colors).Solid lines show results for the triangulation approach of CLaMS, dashed lines show results for the triangulation approach of ATLAS.Shaded areas and solid black lines show vertical diffusion coefficients from other studies for comparison:Massie and Hunten (1981),Fukao et al. (1994),Balluch and Haynes (1997),Legras et al. (2005).Values for λ c =10 d −1 and r 0 =300 km are zero and not shown.Right: same for the effective horizontal diffusion coefficient.Grey lines are the vertical diffusion coefficients from the left panel scaled by α 2 .Solid black line shows estimate fromWaugh et al. (1997).
Fig. B2.The parameter γ only for r 0 =100 km.The solid lines show the definition for γ from this study (Eq.16), the different dashed lines show the definition of γ inKonopka et al. (2004) (Eq.14)   with different approaches to obtain the model points m model i : nearest neighbor of measurement, average over first neighbors of the measurement, average over first and second neighbors.The blue lines show results if only using the ACATS measurements (about 20 km resolution) and the grey lines if only using the Argus and ALIAS measurements (about 1 km resolution).
tropopause is needed for several purposes in the model.It is set to the 2 PVU potential vorticity surface in the extratropics (dynamical tropopause) and to the 375 K potential temperature surface in the tropics.The model code is platform independent and can run in parallel on a simple personal computer network connected via the internet.Main modules are implemented in Matlab and time critical code in C. Extensive use of free software libraries is made, the most important ones are Qhull for the triangulation, KPP (Kinetic Preprocessor) for the chemistry module and a Matlab implementation of a parallel MPI like interface based on exchange via ssh.
The depth of the boundary layers is again based on z (Sect.2.3.3).All parcels in these boundary layers are deleted in every model step and replaced by newly initialized parcels.The reinitialization of the parcels is necessary to avoid that holes or volumes with parcels outside the model boundary appear due to the large-scale motion.Parcels outside the model domain are deleted.Tracers and chemical species are interpolated from the old parcel positions to the new parcels where possible.Where not possible, tracers and chemical species are set to climatological values.This may e.g.be the case where the downward branch of the Brewer-Dobson circulation produces volumes void of any parcels directly below the upper boundary.A www.geosci-model-dev.net/2/153/2009/Geosci.Model Dev., 2, 153-173, 2009 runs is started with different mixing parameterizations.Each run simulates the period from 1 November 1999 to 15 March 2000.The methane and Halon tracers are initialized at 1 December.Between January and March, the model tracers are compared to in-situ measurement data from several high altitude flights conducted during the campaign that are able to temporally resolve fine tracer structures.