Open Access

Abstract. A new hybrid Eulerian–Lagrangian numerical scheme (HEL) for solving prognostic equations in fluid dynamics is proposed. The basic idea is to use an Eulerian as well as a fully Lagrangian representation of all prognostic variables. The time step in Lagrangian space is obtained as a translation of irregularly spaced Lagrangian parcels along downstream trajectories. Tendencies due to other physical processes than advection are calculated in Eulerian space, interpolated, and added to the Lagrangian parcel values. A directionally biased mixing amongst neighboring Lagrangian parcels is introduced. The rate of mixing is proportional to the local deformation rate of the flow. The time stepping in Eulerian representation is achieved in two steps: first a mass-conserving Eulerian or semi-Lagrangian scheme is used to obtain a provisional forecast. This forecast is then nudged towards target values defined from the irregularly spaced Lagrangian parcel values. The nudging procedure is defined in such a way that mass conservation and shape preservation is ensured in Eulerian space. The HEL scheme has been designed to be accurate, multi-tracer efficient, mass conserving, and shape preserving. In Lagrangian space only physically based mixing takes place; i.e., the problem of artificial numerical mixing is avoided. This property is desirable in atmospheric chemical transport models since spurious numerical mixing can impact chemical concentrations severely. The properties of HEL are here verified in two-dimensional tests. These include deformational passive transport on the sphere, and simulations with a semi-implicit shallow water model including topography.


Introduction
Numerical chemical weather forecast systems and Earth system models include components describing the chemistry, including aerosols, and the interaction of these with cloud and radiation processes (e.g., Grell et al., 2005;Pozzoli et al., 2008).The introduction of many more prognostic variables, sometimes several hundred, representing the concentrations of the individual chemical species, poses some severe challenges regarding computational methodologies.

Conservation of invariant quantities such as mass
Published by Copernicus Publications on behalf of the European Geosciences Union.
7. Consistency between wind and mass fields (minimize the so-called mass-wind inconsistency problem) 8. Compatibility (mixing ratios should be bound by their upstream values) 9. Preservation of constant mixing ratios in deformational flows 10.Preservation of linear correlations between constituents For a brief discussion of these see Machenhauer et al. (2008).
While several of the above-listed desired properties are of particular relevance in atmospheric chemical transport models, most are also highly desirable/necessary when it comes to simulation of geophysical dynamics, not least those components involving water vapor, liquid water droplets and ice crystals in the atmosphere, or e.g.salinity in the ocean.
There is, however, one additional property, not listed above and less discussed in the literature, which is particularly important for chemistry and chemistry-climate applications: 11.Avoidance of spurious numerical mixing/unmixing (Lauritzen and Thuburn, 2012).

Mixing and unmixing in Eulerian-based models
The 11th property above refers to the ability of a scheme to preserve pre-existing functional relations between tracers.
Mixing or unmixing can be divided into three categories (see Fig. 1), real mixing, range-preserving unmixing, and overshooting.
In their native forms, most transport schemes operating on a fixed Eulerian grid (including semi-Lagrangian schemes) will lead to numerical mixing between tracers which, depending on the scheme, can fall into any of the three categories.We will term this implicit mixing.
Obviously, range-preserving unmixing and overshooting are unphysical processes, but even real mixing can be so.At a macro spatial scale corresponding to the grid distance in a fluid dynamical model, mixing in nature can be molecular, or a result of turbulence, i.e., deformations of the flow.In geophysical fluid dynamics (GFD) and for typical grid resolutions the former is several orders of magnitude smaller than the latter, and, thus, molecular mixing/diffusion is neglected in the governing model equations.To the extent diffusion is parameterized explicitly in GFD models it is therefore supposed to represent the mixing associated with unresolved deformations in the flow.In non-deformational flow no mixing should take place and any functional relation between tracers should therefore remain unchanged for inert tracers; i.e., all points in a diagram like that in Fig. 1 will keep their initial positions.Unfortunately, this is not the case for most traditional transport schemes, and therefore they are subject to spurious mixing, even if it belongs to the "real" category.In Eigil Kaas: A hybrid Eulerian Lagrangian numerical scheme.Any new relative concentrations (ξ k ,χ k ), generated by the transport scheme, can be represented as a point.If the point falls within the shaded convex hull, it is classified as real mixing.If within the dashed rectangle but outside the shaded area it is classified as range preserving unmixing, and, finally, if outside the dashed rectangle it is classified as overshooting (adopted from Lauritzen and Thuburn, 2012).
with such schemes, and in general deformational flows, the Lagrangian parcels become irregularly distributed in space and strongly deformed.The irregular distribution is a fundamental problem for the following reasons: 1. Dynamical tendencies due to non-advective processes 170 generally need to be estimated on a regular grid to ensure consistency with and between the governing prognostic equations.
2. In practice, it is problematic to calculate parameterized physical processes on an irregular grid.
175 3.There may be sub-domains with no or very few parcels.
The result of an interpolation from mass centers of Lagrangian parcels to a regular Eulerian grid will typically, after some time, show a completely unrealistic "spotty" distribution because neighboring parcels originate from quite 180 different positions at the initial time of the integration (see e.g.lower right panel of Fig. 11).In nature finite size Lagrangian parcels typically deform into very long, very thin filaments.A Lagrangian transport scheme aiming at resolving the developing shape of individual parcels without any 185 explicit mixing between neighboring parcels would represent the discrete Lagrangian analogy to the pseudo-spectral method.The aliasing of such a Lagrangian scheme is realised when one interpolates from Lagrangian to physical space, i.e., to a fixed Eulerian grid.To avoid such aliasing 190 one can introduce mixing between parcels.One may think of the difference between a mixed and unmixed Lagrangian scheme as an analogy to the difference between spectral and pseudo-spectral schemes.Any new relative concentrations (ξ k , χ k ), generated by the transport scheme, can be represented as a point.If the point falls within the shaded convex hull, it is classified as real mixing; if within the dashed rectangle but outside the shaded area it is classified as rangepreserving unmixing; and, finally, if outside the dashed rectangle it is classified as overshooting (adopted from Lauritzen and Thuburn, 2012).
the case where the tracers are chemically active, this can potentially be a serious problem as spurious chemical reactions are then initiated, or chemical equilibria are displaced.
Note that apart from the initial truncation numerical methods based on orthogonal series expansion functions are the only Eulerian type numerical schemes which do not introduce numerical mixing in the case of non-deformational flow.However, generally filters must be introduced in such schemes to prevent e.g.development of negative values, and this introduces numerical mixing also in regions of (quasi-)linear flow.
In the more realistic case of spatially varying flow Lagrangian parcels will deform into thinner and thinner filaments, which in nature are finally mixed via molecular mixing.An important question is to what extent explicit numerical diffusion/mixing is required as a supplement to that implied by the native version of some numerical scheme in order to mimic the cascade into small scales correctly.For typical grid point/cell-based methods, including semi-Lagrangian schemes, some inherent numerical mixing is almost unavoidable, and this may be sufficient to control the cascade in a statistical sense.In Galerkin methods -e.g. the classical spectral method (Machenhauer, 1979) -the gradual development of non-resolved filaments and structures is controlled by demanding the residual to be orthogonal to the resolved expansion functions (see e.g.Machenhauer, 1979, or Durran, 2010).This gives rise to an implied mixing, which, depending on the chosen expansion functions, is generally non-local in physical space.Explicit horizontal diffusion is required for most Galerkin schemes in order to prevent so-called "spectral blocking" (Machenhauer, 1979), i.e., spurious accumulation of energy on the shortest resolved scales.The situation is different in pseudo-spectral models where considerably stronger scale selective damping is necessary to avoid aliasing.
In conclusion explicit mixing, in terms of filters, diffusion, spectral damping, etc., is needed in both grid point/cell methods and methods based on series expansion in order to ensure shape preservation, and in particular positive definiteness, and to control the cascade into smaller scales in a realistic way.In general the combined implicit and explicit mixing will not represent true physical mixing although it may be real in the sense described in Fig. 1.

The HEL approach
Here we present a numerical methodology, termed the hybrid Eulerian-Lagrangian (HEL) numerical scheme, which has been designed to fulfill as many as possible of the desired properties mentioned in Sect.1.1.The aim is to combine the Eulerian and the Lagrangian approaches in such a way that the main problems related to either of these are eliminated or at least reduced.The ideas behind HEL have been inspired by other Lagrangian approaches, in particular that of ATTILA (Atmospheric Tracer Transport In a Lagrangian Atmospheric model) (Reithmeier and Sausen, 2002;Stenke et al., 2008Stenke et al., , 2009)), which uses a fully Lagrangian scheme for all tracers, including water substances.
Contrary to the problems mentioned in the previous subsection for Eulerian-based schemes, fully Lagrangian schemes are formally exact for the pure advection problem assuming trajectories have been calculated exactly.However, with such schemes, and in general deformational flows, the Lagrangian parcels become irregularly distributed in space and strongly deformed.The irregular distribution is a fundamental problem for the following reasons: 1. Dynamical tendencies due to non-advective processes generally need to be estimated on a regular grid to ensure consistency with and between the governing prognostic equations.
2. In practice, it is problematic to calculate parameterized physical processes on an irregular grid.
3. There may be sub-domains with no or very few parcels.
The result of an interpolation from mass centers of Lagrangian parcels to a regular Eulerian grid will typically, after some time, show a completely unrealistic "spotty" distribution because neighboring parcels originate from quite different positions at the initial time of the integration (see e.g.lower right panel of Fig. 11).In nature, finite size Lagrangian parcels typically deform into very long, very thin filaments.A Lagrangian transport scheme aiming at resolving the developing shape of individual parcels without any explicit mixing between neighboring parcels would represent the discrete Lagrangian analogy to the pseudo-spectral method.The aliasing of such a Lagrangian scheme is realized when one interpolates from Lagrangian to physical space, i.e., to a fixed Eulerian grid.To avoid such aliasing one can introduce mixing between parcels.One may think of the difference between a mixed and unmixed Lagrangian scheme as an analogy to the difference between spectral and pseudospectral schemes.
In the native form of HEL without parcel mixing the density and, optionally, other prognostic variables are known at all times via a fully Lagrangian as well as a traditional Eulerian representation.At each time step a nudging technique is applied where the density information in the downstream translated Lagrangian parcels is used to modify or "repair" an Eulerian-based advection.In this way the non-dispersive, non-diffusive, and shape-preserving advantages of the Lagrangian method can be adopted in an otherwise diffusive and/or strongly dispersive Eulerian-based transport scheme.Physical tendency contributions not related to pure advection are most obviously and accurately calculated in Eulerian space and subsequently interpolated and added to the Lagrangian parcel values.In this way the bulk of the model history is kept in the Lagrangian representation.
With this generic design one is, however, still left with the problem of aliasing in Lagrangian space.To reduce or eliminate aliasing a real mixing, described in detail below, is introduced in Lagrangian space between neighboring parcels.To mimic the physical nature of the cascade into smaller scales, we let the degree of mixing depend on the local flow deformation rate.Not surprisingly, the introduction of such mixing turns out to be instrumental for a proper description of the spectral distribution of prognostic variables; however, this issue is not covered in the present paper.
Although we are not dealing with chemistry here, it is noted that calculations of chemical reactions are naturally performed in the Lagrangian parcels where it is known that only real mixing takes place.

Other relevant Lagrangian type numerical schemes
HEL has some similarities to particle in cell (PIC) methods, which have been used extensively in e.g.plasma physics (Tskhakaya et al., 2007).However, one among several fundamental differences is that in PIC methods the number of Lagrangian parcels far exceeds that of Eulerian grid points/cells, while in HEL these numbers are equal or at least close to the same order of magnitude.The choice of relatively few Lagrangian parcels in HEL was motivated by efficiency considerations since computationally expensive chemical reactions (not dealt with here) involving up to several hundred chemical tracers should be performed in Lagrangian space.
As indicated above the transport and mixing in the HEL scheme is similar to other Lagrangian parcel methods.In Reithmeier and Sausen (2002), and later improved in Stenke  (Frömming et al., 2011).ATTILA does not handle the "dry dynamics" as opposed to HEL; however, besides this the HEL and ATTILA approaches are quite similar when applied for nondry-dynamical prognostic variables.One difference, though, is that ATTILA on average holds more Lagrangian parcels per Eulerian grid cell than the version of HEL presented here, which on average only has one.More importantly, there are some differences in the way horizontal mixing between neighboring parcels is performed in ATTILA and HEL.For both schemes the degree of mixing depends on the horizontal shear deformation rate of the flow; however, in ATTILA this is a simple analytical expression based on Smagorinsky (1963), whereas in HEL the deformation of each parcel is kept as an additional prognostic variable, which is increased each time step in proportion to the shear deformation rate, and attempted to be reduced via realized mixing with neighboring parcels.
The approaches in the CLaMS model (McKenna et al., 2002) are less similar than ATTILA to those applied in HEL.The mixing in CLaMS is based on a dynamically adaptive grid and it becomes active in terms of mass exchange between neighboring parcels when the flow deformation is high.A local, in time and space, Lyapunov exponent is used to determine the degree of mixing that takes place, which in practice takes place via generation of new Lagrangian parcels in strongly deformed flow, or merging of clustered Lagrangian parcels.This is one main difference compared to HEL and ATTILA, where Lagrangian parcels survive throughout the integration.In the ATLAS model (Wohltmann and Rex, 2009) the flow-dependent mixing methodology of CLaMS has been modified with emphasis on better performance in lower-resolution model configurations.Also in FPIC (Kaas et al., 1997) an implied mixing takes place via simple birth and merging of particles.
The so-called "trajectory-tracking scheme" introduced in Dong and Wang (2011) and updated in Dong and Wang (2012) has some similarities to HEL.In two-dimensional problems this scheme treats Lagrangian parcels as polygons with a finite number of edges, and with all Lagrangian parcel polygons spanning exactly the complete integration domain.The Eulerian space representation is obtained via a first-order conservative remapping so that total mass is conserved in the Lagrangian as well as the Eulerian representation.A "curvature-guard" algorithm is applied in order to maintain an accurate polygon representation in deformation flows.However, this algorithm does not lead to any mixing, as in HEL, between neighboring Lagrangian parcels.Therefore, in long-term simulations, one should expect problems equivalent to aliasing.

Overview
The paper is organized as follows: Sect. 2 provides a generic description of the HEL approach, i.e., HEL without any type of mixing between parcels, while Sect. 3 describes how mixing between adjacent parcels is achieved.Section 4 presents passive inert transport tests on the sphere.Various traditional and more recently proposed error measures and evaluation statistics are used to demonstrate the performance of HEL in both solid body rotation flow (Sect.4.1) and in strongly deformational flow (Sects.4.2 and 4.3).Section 5 deals with some initial attempts to implement HEL as the basis for a dynamical core in a geophysical fluid dynamics model.In this case the test bed is a shallow water model in plane geometry.Finally Sects.6 and 7 discuss the results, including some outlooks for future work, and summarize the basic findings.

HEL -passive transport
To introduce the procedure followed in HEL in more detail we first consider the continuity equations for a set of M tracers with densities ρ m : or alternatively in Lagrangian form, where V is the flow velocity vector.For simplicity we have ignored any sources and sinks, and any diffusion in Eqs.(1) or (2).
In geophysical fluid dynamics Eqs.(1) or ( 2) is normally solved via finite volume (FV) methods operating on a fixed Eulerian grid.Two different families of FV methods have been applied: flux based and cell integrated.In flux-based methods the fluxes of mass swept through each face of predefined Eulerian grid cells within a time step are calculated first, and the change in density is then determined from the net inflow of mass into each grid cell.In cell-integrated methods a semi-Lagrangian upstream departure cell is first identified.The estimated total mass in this upstream cell then determines the corresponding arrival (Eulerian) grid cell density.For a more detailed description of the difference between the two families, see Machenhauer et al. (2008), where it is also demonstrated that they are in fact numerically equivalent.For a general review of FV methods, we refer to LeVeque (2002) and Eymard et al. (2000).It should be noted that other non-FV, yet still massconservative, schemes have also been used.For the present work it is of particular relevance to mention the so-called locally mass-conserving semi-Lagrangian method (LMCSL) (Kaas, 2008), which is based on a simple partition of unity principle.
Depending on the chosen order of accuracy any numerical method -Eulerian or semi-Lagrangian type -applied for solving Eqs. ( 1) or (2) on an Eulerian grid will suffer from some degree of dissipation (or possibly anti-dissipation), some numerical dispersion, and generally, for higher-order schemes, the solution will not be shape preserving.
As mentioned above the main motivation behind the present work is to use a fully Lagrangian forecast, run in parallel, to modify the Eulerian grid forecast in such a way that the above-mentioned disadvantages are reduced or eliminated.A purely Lagrangian forecast describes the temporal evolution of the densities of individual Lagrangian fluid parcels as they move around.Formally it is straightforward to integrate Eq. ( 2) for a Lagrangian parcel from time t to some future time t + t: where r(t) is the position vector of the parcel at time t, and D represents the average divergence along the trajectory from r(t) to r(t + t).From Eq. (3) one immediately gets ρ (r(t + t), t + t)) = ρ (r(t), t)) exp tD ; i.e., the effect of divergence over the period from t to t + t is an expansion/contraction factor: multiplied by the original parcel density at time t.
We will now consider the actual numerical discretization of the prognostic equations in Eulerian and Lagrangian space.The Lagrangian parcels are introduced at the initial time, and in the present formulation of HEL they survive throughout the model integration.Also, in the version of HEL presented here, the total number of parcels P and the number of Eulerian grid cells K are equal.At the initial time step, n = 0, Lagrangian parcel densities, L ρ, are initialized by the corresponding values in Eulerian grid cell centroids; i.e., where p, k = 1, . . ., P (= K), and m counts the individual tracers as in Eqs.
(1) and (2).In the following we generally use upper-left superscripts L and E to indicate Lagrangian and Eulerian space representation, respectively.An upperright index denotes the time step.A list of all prognostic variables in HEL to be described below can be found in Appendix A. Assume some numerical scheme has been used to solve Eqs.(1) or (2) in the Eulerian grid cell representation, and let denote the forecast in Eulerian grid cells k = 1, . . ., K at time step n + 1.In general we use the notation (• ) to represent some provisional approximate value.This is also the case here where ρ indicates that the forecast at time step n + 1 is only a provisional Eulerian space forecast to be modified by densities in the Lagrangian representation.
Letting parcels follow downstream trajectories from time step n to n+1 estimated from the actual velocity components in the Eulerian grid, one obtains an approximate Lagrangian solution to the pure advection problem.However, in a general divergent flow the parcel volume density will of course undergo changes.According to Eqs. ( 4) and ( 5) the effect of divergence for parcel p from time step n to n + 1 can simply be modeled as where superscript n + 1/2 indicates that the expansion/contraction factor represents the effect of divergence from time step n to n+1.In practice σ n+1/2 is determined in Eulerian space from the provisional Eulerian space forecast of the "dry air" density, E ρn+1 , i.e., including the effect of divergence, and from a corresponding purely advective forecast, i.e., not including the effect of divergence, which we term E , advρ n+1 : For cell-integrated FV and the LMCSL schemes applied in the present paper the calculation of E σ n+1/2 is straightforward since estimation of E , advρ n+1 is an inherent part of these schemes.
Once new downstream parcel positions L r n+1 have been found, L σ n+1/2 can be obtained via interpolation from E σ n+1/2 .The parcel forecast for parcel p, including divergence, then simply becomes It is important to note that in a dynamical model, in order to prevent numerical instabilities related to fast modes, e.g.gravity waves, E σ n+1/2 p must be based on divergence obtained with a numerically stable scheme.The modification of the Eulerian densities using the parcel densities is done by first interpolating the irregularly spaced parcel densities to the Eulerian grid obtaining certain Eulerian space target values, E T and E T m , m = 1, . . ., M, and then nudging the original Eulerian-based forecast towards these values under the constraints of mass conservation and shape preservation (see details in Sect.2.6).
A generic recipe in an atmospheric multi-tracer application of HEL, not yet considering the mixing, can be summarized as follows at a given time step n (here omitting indices for Eulerian grid cells, k, and Lagrangian parcels, p): 1. Perform a conventional, preferably inherently massconserving, Eulerian or semi-Lagrangian time step of total dry density, E ρn+1 , valid on an Eulerian grid.This is termed the provisional forecast.

Perform a corresponding purely advective time step
in Eulerian space E , advρ n+1 of the dry density, and use this to calculate the divergent multiplication factor, E σ n+1/2 = E ρn+1 / E , advρ n+1 .3. For all tracers, m = 1, . . ., M, perform a provisional Eulerian space forecast, E ρn+1 m .4. Perform a pure downstream displacement of the irregularly spaced Lagrangian parcels; i.e., calculate downstream trajectories and reposition each parcel from L r n to L r n+1 .
5. Interpolate E σ n+1/2 from the Eulerian grid cells to the positions L r n+1 resulting in values L σ n+1/2 .Then calculate the new parcel densities for both the dry air and all the tracers, including the effect of divergence: 6. Interpolate L ρ n+1 from the Lagrangian grid to obtain the target values, E T , in the Eulerian grid.In Eulerian grid cells with no nearby Lagrangian parcels E T n+1 is set equal to E ρn+1 from step 1; see details in Sect.2.5.7.As step 6 but for all tracers m = 1, . . ., M.

Nudge E
ρn+1 towards E T n+1 under constraints of mass conservation and shape preservation for the density; see details in Sect.2.6.The result is the final HEL forecast, E ρ n+1 , for the dry air in Eulerian space.9.As step 8 but for all tracers m = 1, . . ., M. However, now the constraints are mass conservation and shape preservation for tracer mixing ratios.

The underlying numerical scheme
As outlined above, some numerically stable scheme must be chosen in order to obtain the provisional forecast in Eulerian grid space.For applications in HEL it would be reasonable to use a semi-Lagrangian type scheme since trajectory calculations can be partly re-used for estimation of the downstream parcel trajectories; also, it would then be possible to take long steps not subject to classical CFL conditions for advective processes.Therefore relevant schemes include, e.g., flux-based multidimensional schemes such as Lin and Rood (1996) and Leonard et al. (1996), the Departure area Cell-Integrated Semi-Lagrangian (DCISL) scheme (Machenhauer and Olk, 1997), the Conservative Semi-LAgrangian Multi-tracer transport scheme (CSLAM) (Lauritzen et al., 2010), the Semi-Lagrangian Inherently Conserving and Efficient scheme (SLICE) (Zerroukat et al., 2002), or the Locally Mass-Conserving Semi-Lagrangian scheme (LMCSL) (Kaas, 2008).Note, however, that any mass-conserving scheme for solving continuity equations can in principle be used as the underlying scheme for HEL.In fact, relaxing the mass conservation property, any consistent numerical scheme can be used.
In the present paper we have tested the use of first-and third-order versions of the CSLAM and LMCSL schemes to obtain the first guess forecast in Eulerian space.

Estimation of trajectories
The downstream displacement of parcel locations obviously is an essential component in HEL.In simple numerical tests such trajectories can be calculated analytically, or, as in dynamical models, they can be calculated via an iterative procedure, which is equivalent to that used in traditional semi-Lagrangian models for estimating the upstream departure points/cells.In the present applications we have generally used analytical or approximate analytical trajectories for idealized numerical tests, and iteratively estimated trajectories in dynamical model implementations.This is described further in Sects.4.1, 4.2, and 5.1.1.

Update of the parcel volumes
Each Lagrangian parcel represents a certain volume L V of the fluid which, at the initial time step, is simply initialized as the volume represented by the volume of the relevant Eulerian grid cell.
Once L ρ n+1 has been calculated one can update the volume of each parcel.Omitting the upper index "L" the parcel volume V n+1 p at time step n + 1 is determined diagnostically by the constraint that, in the absence of mixing, the total mass of each Lagrangian parcel is conserved, i.e.,  10).The volume of the parcels is a necessary ingredient for mixing between neighboring parcels (Sect.3).

Interpolations from Eulerian to Lagrangian representation
Interpolations between the Eulerian and Lagrangian representations are required as part of the HEL scheme.For example, as explained above, it is necessary to interpolate σ to the Lagrangian representation.Similarly, when HEL is used in a dynamical model, tendencies related to other physical processes must be interpolated.
In the present formulation of HEL all interpolations from the Eulerian grid cells to the Lagrangian parcel locations are fourth-order Lagrange polynomial interpolations.

Target values
Provisional target values, E T n+1 , for the Eulerian space dry air density can be obtained when E ρn+1 and L ρ n+1 have been calculated.The provisional target value in Eulerian grid cell k is composed as a weighted sum of the provisional Eulerian-based forecast, E ρn+1 k , and a parcel-based estimate R k : where R k is defined as where w p,k is a simple bi-linear interpolation weight given to an estimate of the density E ρn+1 p,k in Eulerian grid cell k which is based on the density at the location of parcel number p: g p,k is a second-order numerical approximation to the gradient ∇ E ρ n+1 at the location 0.5( L r n+1 p + r k ), and r k is the position vector of the kth Eulerian grid cell.
The weights w 1 and w 2 in Eq. ( 13) are determined as follows: where w 0 is a small positive number, and H k is a measure of the homogeneity of the distribution of Lagrangian parcels around the Eulerian cell k.Here we have used the following estimate of H k : Eigil Kaas: A hybrid Eulerian Lagrangian numerical scheme.7 in a dynamical model, tendencies related to other physical processes must be interpolated.
In the present formulation of HEL all interpolations from the Eulerian grid cells to the Lagrangian parcel locations are fourth order Lagrange polynomial interpolations.

Target values
Provisional target values, E T n+1 , for the Eulerian space "dry air" density can be obtained when E ρn+1 and L ρ n+1 have been calculated.The provisional target value in Eulerian grid cell k is composed as a weighted sum of the provisional Eulerian based forecast, E ρn+1 k , and a parcel based estimate R k : where R k is defined as where w p,k is a simple bi-linear interpolation weight given to an estimate of the density E ρn+1 p,k in Eulerian grid cell k which is based on the density at the location of parcel number p: g p,k is a second order numerical approximation to the gradient ∇ E ρ n+1 at the location 0.5( L r n+1 p + r k ) and r k is the position vector of the kth Eulerian grid cell.
The weights w 1 and w 2 in Eq. ( 13) are determined as follows: where w 0 is a small positive number, and H k is a measure of the homogeneity of the distribution of Lagrangian parcels around the Eulerian cell k.Here we have used the following estimate of H k where L k is a subset of Eulerian cells including k and its nearest eight surrounding neighbors (for a regular grid as used here).In Eq. ( 16) the maximization of w 1 is introduced to ensure that the target values are always well-defined.The value of w 0 has, somewhat arbitrarily, been set to 10 −5 , i.e., in the target values estimated from Eq. ( 13) there will always be some weight on the Eulerian based forecast.In general this weight is small but in the special case where W k = 0, i.e. there are no neighboring parcels at all, the target value will simply be equal to the provisional Eulerian forecast.
In practice the calculations of R k and W k are not per-480 formed as sums over all parcels as indicated in Eq. ( 14) since this would be very inefficient numerically because it is known that only four of the parcel weights, w p,k , for a given parcel p are different from zero.Instead R k and W k are calculated in a single loop over all parcels where information is 485 distributed (summed) to the neighboring four Eulerian cells, followed by a second loop where the result is divided by the sum of weights for each cell.
Provisional target values for the tracers E T n+1 m , m = 1,...,M are obtained via the same technique as outlined for 490 the "dry air", i.e., E T n+1 above.I.e.all weights are the same.
To obtain an Eulerian space forecast, which is shape preserving and compatible we must identify minimum and maximum permitted mixing ratios, q − and q + , respectively, for 495 each tracer m in each Eulerian grid cell.The mixing ratio for a tracer m is defined as q m = ρ m /ρ, implying that we can always deduce mixing ratios in both the Eulerian and Lagrangian representation when ρ m and ρ are known, and we can always convert a mixing ratio q m back to volume den-500 sity when ρ is known.The upstream Eulerian mixing ratios at time level n are used by selecting the minimum and maximum mixing ratios in the four grid cells k 1 ,...,k 4 surrounding the location of the upstream departure point for the trajectory, which at time level n+1 is located at the cell centroid of 505 cell k, i.e., the departure point in a classical semi-Lagrangian context (see Fig. 2).
Thereby the provisional minimum and maximum mixing ratios for tracer m become: The final limits are obtained by using additional information from mixing ratios in the Q parcels p i ,i = 1,...,Q that have Euclidian distance to the Eulerian cell k at time level n + 1, where L k is a subset of Eulerian cells including k and its nearest eight surrounding neighbors (for a regular grid as used here).In Eq. ( 16) the maximization of w 1 is introduced to ensure that the target values are always well defined.The value of w 0 has, somewhat arbitrarily, been set to 10 −5 ; i.e., in the target values estimated from Eq. ( 13) there will always be some weight on the Eulerian-based forecast.In general this weight is small but in the special case where W k = 0, i.e., where there are no neighboring parcels at all, the target value will simply be equal to the provisional Eulerian forecast.
In practice the calculations of R k and W k are not performed as sums over all parcels as indicated in Eq. ( 14) since this would be very inefficient numerically because it is known that only four of the parcel weights, w p,k , for a given parcel p are different from zero.Instead R k and W k are calculated in a single loop over all parcels where information is distributed (summed) to the neighboring four Eulerian cells, followed by a second loop where the result is divided by the sum of weights for each cell.
Provisional target values for the tracers E T n+1 m , m = 1, . . ., M are obtained via the same technique as outlined for the dry air, i.e., E T n+1 above.I.e., all weights are the same.
To obtain an Eulerian space forecast, which is shape preserving and compatible we must identify minimum and maximum permitted mixing ratios, q − and q + , respectively, for each tracer m in each Eulerian grid cell.The mixing ratio for a tracer m is defined as q m = ρ m /ρ, implying that we can always deduce mixing ratios in both the Eulerian and Lagrangian representation when ρ m and ρ are known, and we can always convert a mixing ratio q m back to volume density when ρ is known.The upstream Eulerian mixing ratios at time level n are used by selecting the minimum and maximum mixing ratios in the four grid cells k 1 , . . ., k 4 surrounding the location of the upstream departure point for the trajectory, which at time level n+1 is located at the cell centroid of cell k, i.e., the departure point in a classical semi-Lagrangian context (see Fig. Thus the provisional minimum and maximum mixing ratios for tracer m become The final limits are obtained by using additional information from mixing ratios in the Q parcels p i , i = 1, . . ., Q that have Euclidian distance to the Eulerian cell k at time level n + 1, which is less than 1.5 grid distances: Once the Eulerian space forecast for dry air density E ρ has been obtained, the minimum and maximum permitted values of volume density for tracer m = 1, . . ., M can easily be obtained:

Nudging of first guess towards target values
After the (shape-preserving) target values have been calculated, the provisional, "first guess", Eulerian forecast can be corrected.The correction, or nudging, is done in two steps.First we calculate the total mass, M > 0: Mass of target field is too large.
In the (extremely unlikely) event that the mass of the target field is exact ( M = 0), the Eulerian field is simply replaced by the target field.In the two remaining possibilities the target field has to be modified to ensure mass conservation.If the mass of the target field is less that the actual mass ( M < 0), we calculate the maximum possible mass of the field, i.e., , where shape preservation is still fulfilled.The target field and the maximum field can then be combined to produce a final corrected Eulerian field, which is both mass conserving and shape preserving.
This is always possible as M + ≥ M. The procedure is the same if the mass of the target field is too large; then the target field is weighed with the minimum mass field, The procedure is then repeated for all tracers m = 1, . . ., M.
The nudging employed in the current version is global, but since all chemistry is calculated in the Lagrangian parcels, it will not introduce errors due to the inevitable numerical mixing in the Eulerian domain.A local nudging method has also been tested, but it leads to somewhat poorer -and in practice less localized -results, as the nudging method's ability to correct the Eulerian values will be lessened by hard locality constraints.The traditional concerns when using global methods, i.e., mass redistribution and unphysical mixing, is fully controlled by the correct values being preserved in the Lagrangian parcels, meaning that HEL is very "local".

Mixing between parcels
As discussed in Sect. 1 densities/mixing ratios or other invariants will generally develop into thin filaments as part of the cascade into smaller and smaller scales in geophysical deformational flows.At some point a model at given resolution cannot represent the spatial scale of the filaments, and explicit horizontal diffusion may therefore be required to prevent spectral blocking.Considered in discrete Lagrangian space, i.e., a model, the analogue to spectral blocking is realized as a gradual development into unrealistically large differences between densities/mixing ratios in neighboring parcels.
An example of this is presented below in Sect. 5. Therefore, due to their non-dissipative nature, explicit mixing must be introduced in Lagrangian models.We introduce a directionally biased mixing as an a posteriori operation applied each time step after the generic HEL forecast, described above, has been obtained.
In the present paper we only consider two-dimensional flow.In this case the degree of mixing between neighboring parcels is based on a modified instantaneous and local two-dimensional rate of deformation: where u n and v n , respectively, are the velocity components in the two directions spanned by coordinates x and y at time step n.It can be seen that the effect of divergence (last term) is subtracted from the traditional expression for deformation rate.This is done because the mixing we want to introduce Tru e as ym pto tic sh ea r dila tat ion ax is a n+1 a n+1

Time step n+1
Normalize at time step n+1 T ru e a s y m p to ti c s h e a r d il a ta ti o n a x is Tru e as ym pto tic sh ea r dila tat ion ax is For each parcel P and P a denotes, respectively, its position and the position of the auxiliary parcel.The associated auxiliary vectors are marked black and gray, respectively, to distinguish between the two parcels.In this example the (true) asymptotic / Lagrangian dilatation axis, indicated with dashed lines, is the same all over the small domain shown, although its direction rotates clockwise from time step n to time step n + 1.Note that it is the relative movements between parcels and auxiliary parcels that are relevant.After several timesteps the auxiliary vectors will align approximately with the true asymptotic dilatation axis.
Asy mp toti c (La gra ngi an) dila tati on axis . This figure describes why we use a factor of 0.25 between the deformation within one time step ∆t • D and the fraction of the parcel area ∆ 2 that should be mixed.The underlying assumption is that the deformation is mainly due to shear, and the figure shows how cells are then deformed within one time step from regular squares into irregular parallelograms with the same area (assuming zero divergence).In order to re-shape a cell into squares the mass in the two shaded triangles must be given off via interpolations to the neighboring cells.The factor comes because the sum of the triangular areas is 0.25∆t where it is again noted that for regular grids the sum over k only includes four grid cells surrounding parcel p.The parcel deformation is finally reduced according to the degree of mixing that has actually taken place: where the factor of "4" is a constant determining how much L δ p is reduced per unit change in relative area.This factor is obtained from the same geometrical considerations mentioned above (see Fig. 4) of the relationship between deformation rate and the fraction of the parcel volume which is actually mixed with neighboring parcels.Note that the factor of 0.25 in Eq. ( 29) ensures that L δ n+1 p cannot be less than zero.
Most computational operations in the above list are common for all tracers and therefore, in multi-tracer applications, the total number of operations is limited.
The HEL-mixing is quite different from the uncontrollable, and in many cases unrealistic, numerical mixing/unmixing (see Fig. 1), which is introduced in most tradi-i.e., Eulerian and semi-Lagrangian type methods.This is be-640 cause in such models the degree of mixing is different from tracer to tracer, because it depends on the spatial roughness of their density fields.The degree of mixing applied in HEL between Lagrangian parcels is the same for all tracers.Thus, in the Lagrangian space representation, i.e., for the parcels, 645 the problem of numerical mixing and unmixing is eliminated.It is noted that the introduction of a flow-dependent mixing based on the degree of deformation is not new.Sadourny and Maynard (1997) introduced a horizontal diffusion which was dependent on the magnitude of the deformation rate of 650 the flow, and later Váňa et al. (2008) used a similar approach to obtain a flow dependent degree of mixing in a semi-Lagrangian model.Also ATTILA (Stenke et al., 2009) and CLaMS (Chemical Lagrangian Model of the Stratosphere) (McKenna et al., 2002) employ mixing depending on hori-655 zontal flow deformation rate.Fig. 3. Downstream advection of two different Lagrangian parcels and their auxiliary parcels in a pure shear flow.For each parcel P and P a denote, respectively, its position and the position of the auxiliary parcel.The associated auxiliary vectors are marked black and gray, respectively, to distinguish between the two parcels.In this example the (true) asymptotic/Lagrangian dilatation axis, indicated with dashed lines, is the same all over the small domain shown, although its direction rotates clockwise from time step n to time step n + 1.Note that it is the relative movements between parcels and auxiliary parcels that are relevant.After several time steps the auxiliary vectors will align approximately with the true asymptotic dilatation axis.
should not lead to excessive damping when the scheme is applied to the full mass field in a dynamical model.Not introducing such a modification would tend to damp dynamically important gravity waves.
The calculation of D is estimated in Eulerian space via centered, i.e., second-order accurate, differences.The expression in Eq. ( 26) is only valid in Cartesian geometry.Therefore, in the applications on the sphere presented below, metric factors have been applied.
In addition to the prognostic variables discussed above a set of three prognostic variables are introduced in Lagrangian space (only) in order to perform the mixing.For each Lagrangian parcel p these include the actual parcel deformation, L δ p , and the two coordinates of a position vector, L r p a , of a passive auxiliary Lagrangian parcel, L p a , that is used to identify the asymptotic dilatation axis for shear, i.e., the direction of "long-term" parcel stretching due to shear considered in a Lagrangian sense; see e.g.Cohen and Schultz (2005).
For each parcel p its deformation is initialized to zero at the first time step: L δ 0 p = 0, and it is updated as part of the mixing procedure described below.The location of the auxiliary parcel p a for main parcel p is initialized one grid distance, , away from p in an arbitrary direction on the integration plane.At each time step the auxiliary parcel is translated downstream using the same trajectory algorithm as for the main parcels.
The mixing operates as follows: 1. Once the modified deformation rate, D n , has been determined in Eulerian space at a given time step, it can be interpolated to the parcel positions enabling calculation of new provisional parcel deformations: 2. Omitting for simplicity the Lagrangian and parcel indices L and p, and the time index n + 1, let ra denote the pure downstream position vector of the auxiliary parcel for a main parcel, which has downstream position vector r.The final downstream position vector r a of the auxiliary parcel is then defined as i.e., the distance between a parcel and its auxiliary parcel is simply normalized to one grid distance.Two examples of the identification of the a vector are shown in Fig. 3.No formal proof is given here that the vectors r a will actually converge toward the Lagrangian shear dilatation axes.Since the additional parcels a are initialized randomly, a certain time will elapse from the model initial state until realistic dilatation axes are identified for each parcel.However this time is proportional to the deformation rate of the flow, and therefore, when the deformation rate is large, the dilatation axes is also quickly approached.Obviously if the flow is linear, i.e., no mixing is needed, then the dilatation axes of the model parcels maintain their initial random orientations.
3. For each parcel, p, let µ p,k represent a fraction of the parcel volume, L V p , which is assigned to a neighboring Eulerian grid cell centroid k. µ p,k is defined as where d p,k is the distance in units of grid distances from the grid cell centroid to the line parallel to L a n+1 p , which passes through parcel p.The value κ determines the degree of directional bias for the mixing.In the present work κ has been set to 10.As an example, if Pa denotes, respectively, its position and the position of the auxiliary parcel.The associated auxiliary vectors are marked black and gray, respectively, to distinguish between the two parcels.In this example the (true) asymptotic / Lagrangian dilatation axis, indicated with dashed lines, is the same all over the small domain shown, although its direction rotates clockwise from time step n to time step n + 1.Note that it is the relative movements between parcels and auxiliary parcels that are relevant.After several timesteps the auxiliary vectors will align approximately with the true asymptotic dilatation axis.where it is again noted that for regular grids the sum over k only includes four grid cells surrounding parcel p.The parcel deformation is finally reduced according to the degree of mixing that has actually taken place: where the factor of "4" is a constant determining how much L δ p is reduced per unit change in relative area.This factor is obtained from the same geometrical considerations mentioned above (see Fig. 4) of the relationship between deformation rate and the fraction of the parcel volume which is actually mixed with neighboring parcels.Note that the factor of 0.25 in Eq. ( 29) ensures that L δ n+1 p cannot be less than zero.
Most computational operations in the above list are common for all tracers and therefore, in multi-tracer applications, the total number of operations is limited.
The HEL-mixing is quite different from the uncontrollable, and in many cases unrealistic, numerical mixing/unmixing (see Fig. 1), which is introduced in most traditional models based on an Eulerian grid/cell representation, i.e., Eulerian and semi-Lagrangian type methods.This is be-640 cause in such models the degree of mixing is different from tracer to tracer, because it depends on the spatial roughness of their density fields.The degree of mixing applied in HEL between Lagrangian parcels is the same for all tracers.Thus, in the Lagrangian space representation, i.e., for the parcels,  that should be mixed.The underlying assumption is that the deformation is mainly due to shear, and the figure shows how cells are then deformed within one time step from regular squares into irregular parallelograms with the same area (assuming zero divergence).In order to re-shape a cell into squares the mass in the two shaded triangles must be given off via interpolations to the neighboring cells.The factor comes because the sum of the triangular areas is 0.25 t • D 2 .d p,k = 0.5, then the exponential factor in Eq. ( 29) becomes about 0.1; i.e., only grid cells close to the ps' dilation axis are assigned an appreciable fraction of volume.For all Eulerian points k with distances to p larger than , µ p,k is set to zero.I.e., in a regular grid µ p,k is only different from zero for a maximum of four individual values of k.It is ensured that the sum of these four weights does not exceed unity.The factor "0.25" is obtained from geometrical considerations of the relationship between deformation rate and the fraction of the parcel volume which should be mixed with neighboring parcels -see Fig. 4. 4. Once µ p,k is calculated, a total mass contribution, µ p,k L ρ m,p V p , for each tracer m is transferred from the Lagrangian parcel p to Eulerian grid cell k.In other words the average density of mass contributions from all parcels "neighboring" k is where it is noted that µ p,k represents elements in a sparse matrix with non-zero contributions from no or only a few parcels.
5. Now the mixing can be realized by transferring the mixed densities in the Eulerian cells back to the Lagrangian parcels.For parcel p the final mixed density L ρ m,p becomes where we have formally used the notation u,L ρ m,p to indicate the unmixed forecasted density in parcel p resulting from the generic HEL recipe.
Note that not only tracers are mixed using Eq.(31).To ensure full consistence between prognostic variables also the dry air is mixed.
It can easily be shown that the total parcel mass for each tracer and for the complete "dry parcel mass of the atmosphere" are conserved when applying Eq. ( 31).
6. Based on the amount of actual mixing, M p , that has taken place for parcel p via the above operations the final parcel deformation is calculated.The actual mixing, not including trivial mixing of parcels with themselves, is where it is again noted that for regular grids the sum over k only includes four grid cells surrounding parcel p.The parcel deformation is finally reduced according to the degree of mixing that has actually taken place: where the factor of "4" is a constant determining how much L δ p is reduced per unit change in relative area.This factor is obtained from the same geometrical considerations mentioned above (see Fig. 4) of the relationship between deformation rate and the fraction of the parcel volume which is actually mixed with neighboring parcels.Note that the factor of 0.25 in Eq. ( 29) ensures that L δ n+1 p cannot be less than zero.
Most computational operations in the above list are common for all tracers, and therefore, in multi-tracer applications, the total number of operations is limited.
The HEL mixing is quite different from the uncontrollable, and in many cases unrealistic, numerical mixing/unmixing (see Fig. 1), which is introduced in most traditional models based on an Eulerian grid/cell representation, i.e., Eulerian and semi-Lagrangian type methods.This is because in such models the degree of mixing is different from tracer to tracer, because it depends on the spatial roughness of their density fields.The degree of mixing applied in HEL between Lagrangian parcels is the same for all tracers.Thus, in the Lagrangian space representation, i.e., for the parcels,

Passive tracer numerical simulations on a cubed sphere
To validate HEL we perform inert passive tracer transport on the sphere driven by both solid body rotation and two types of deformational flow.
For the passive transport tests presented here the underlying Eulerian-based scheme required in HEL is a first-order (i.e., shape preserving by definition) version of CSLAM (Lauritzen et al., 2010).Where necessary this is referred to as CSLAM-1st.The performance of HEL is compared to that of a third-order accurate version of CSLAM in combination with a simple shape-preserving filter (Lauritzen et al., 2010), referred to as CSLAM-M.
In the tests shown here both HEL and CSLAM have been implemented on a so-called cubed sphere grid -see Lauritzen et al. (2010) for details.

Solid body rotation
In standard solid body rotation tests on the sphere a certain spatial distribution returns to its original position after one or more rotations.In this case HEL should perform with high accuracy since the parcels end up in exactly the same grid cell centroids they were initialized in, and no mixing takes place because the deformation rate of this flow is zero.Thus the target values will be almost exactly equal to their initial value, except for the small weight factor w 0 (see Eq. 16).This again implies that the final HEL forecast should be almost exactly equal to the corresponding analytic solution since global nudging is used.
The test example we present here is the solid body advection of a cosine bell with radius R c = R/3, where R is the radius of the Earth.The angle of rotation is π/4 relative to the Earth rotation axis; i.e., the bell passes over the edges of the cubed sphere.One full revolution is completed in 576 time steps of 1800 s each.These settings are identical to those in Putman and Lin (2007).
To validate the results obtained with this simple setup we use the traditional l 2 and l ∞ error norms.Results in terms of l 2 and l ∞ after one revolution are plotted in Fig. 5.As expected the HEL is very accurate due to the low value of the w 0 weight.The l 2 and l ∞ convergence rates for CSLAM-M are 2.82 and 2.04, respectively, while they are 2.24 and 2.16 for HEL.It is noted that the convergence rate for HEL increases significantly if a fixed Courant number is used in such tests (not shown) because this implies that, as the time step decreases, the weight per time unit on the parcels increases relative to that on the underlying Eulerian-based forecast.While the accuracy increases the shorter the time step in HEL, the opposite is generally the case in semi-Lagrangian models such as CSLAM because the number of remappings increases when the time step is reduced.
The temporal growth of error in HEL and in CSLAM-M are quite different: running over several revolutions the error norms continue to grow in CSLAM-M, although slowly, while in HEL the error norms do not grow with time, as one should also expect from the way HEL is designed.

Deformational flow tests
For the deformation flow tests we have used the two types of analytic flow fields, the density shapes, and the validation diagnostics suggested in Lauritzen et al. (2012), and used in a model intercomparison (Lauritzen et al., 2013).For the diagnostics this means that, in addition to the l 2 and l ∞ error norms and related convergence rates used for the simple solid body rotation tests above, we have calculated an additional set of diagnostics, briefly described below.
The two analytical flow fields used have originally been proposed by Nair and Lauritzen (2010), and they include a non-divergent as well as a divergent flow.In both cases the Lagrangian parcels follow relatively complex trajectories, and the flow is composed of a deformational deformation component, which is different in the two cases, and an overlaid translational flow.The translational part is designed to perform exactly one rotation around the sphere (along equator) during the entire simulation.After a half complete period of simulation the deformational flow component goes to zero, and this part of the flow is then reversed so that the final exact solution equals the initial condition.Half way through the simulation, at the time when the deformational flow component goes to zero and starts to reverse, the initial distributions are deformed into thin filaments, particularly for the non-divergent flow.
Three initial, i.e., t = 0, distributions consisting of two isolated Gaussian hills, two slotted cylinders, and two cosine hills are shown in Fig. 6.Details for these distributions are described in Lauritzen et al. (2012).
As an example Fig. 7 shows the result of simulations at t = T /2, i.e., the most deformed time, and at t = T , i.e., the final time, using the slotted cylinder initial condition and the nondivergent deformation flow.The maximum Courant number is 5.5, and the equatorial resolution 1.5 • in these simulations.At this resolution the final distribution at time t = T obtained with HEL with mixing is considerably closer to the analytic solution than CSLAM-M. . 5. Test of convergence for linear advection on cubed sphere using error norms l 2 (left) and l ∞ (right).ncube refers to the numb d cells in each direction on each of the 6 faces of the cubed sphere, i.e., ncube = 24 corresponds to an equatorial angular grid distanc /(4 × 24) = 0.065 radians.
Using the non-divergent flow, basic error norms l 2 and for CSLAM-M and HEL, and two maximum Courant mbers (1.0 and 5.5), are listed in Tables 1, 2, and 3 for ch of the three initial distributions.In general one can nclude that HEL is very accurate at low resolution for more smooth distributions, i.e., the Gaussian and coe hills, whereas CSLAM-M and HEL are comparable at high Courant number and at the two finest resolutions.other feature is that, as compared to CSLAM-M, HEL is s sensitive to the maximum Courant number.This is be-high Courant number for the smooth distributions while difference is small for the rough slotted cylinder distribu where the convergence rates in any case are low.The con gence rates are comparable in the low Courant number ca 780

"Minimal" resolution
Numerical schemes may be constructed to converge fa at least for smooth distributions.However, as pointed by Lauritzen et al. (2012) increases in resolution are o Fig. 5. Test of convergence for linear advection on cubed sphere using error norms l 2 (left) and l ∞ (right).ncube refers to the number of grid cells in each direction on each of the 6 faces of the cubed sphere; i.e., ncube = 24 corresponds to an equatorial angular grid distance of 2π/(4 × 24) = 0.065 radians.
For illustrative purposes Fig. 7 also includes the result of a simulation without any parcel mixing, and where the parameter H k deliberately has been set to unity.In this case, as can be seen, the distribution at time t = T /2 becomes unrealistic since the initial distribution has been cut into small parts (represented by the individual parcels).However, as expected since time is reversed, the final distribution in this aliased model setup is very close to the analytic solution.It is only small errors in the parcel trajectories, and the fact that w 0 is different from zero, that prevents the final field from being equal to the analytical solution.
Using the non-divergent flow, basic error norms l 2 and l ∞ for CSLAM-M and HEL, and two maximum Courant numbers (1.0 and 5.5), are listed in Tables 1, 2, and 3 for each of the three initial distributions.In general one can conclude that HEL is very accurate at low resolution for the smoother distributions, i.e., the Gaussian and cosine hills, whereas CSLAM-M and HEL are comparable at the high Courant number and at the two finest resolutions.Another feature is that, as compared to CSLAM-M, HEL is less sensitive to the maximum Courant number.This is because HEL is influenced much less by the number of semi-Lagrangian remappings needed to finalize each simulation than CSLAM-M.
The convergence rates for each of the three initial distributions and for the two Courant numbers are listed in Table 4.In general CSLAM-M converges faster than HEL at high Courant number for the smooth distributions, while the difference is small for the rough slotted cylinder distribution where the convergence rates in any case are low.The convergence rates are comparable in the low Courant number cases.

"Minimal" resolution
Numerical schemes may be constructed to converge fastat least for smooth distributions.However, as pointed out by Lauritzen et al. (2012) increases in resolution are often computationally expensive, and therefore it is of interest to identify some kind of measure of absolute accuracy for a numerical scheme.A diagnostic designed for this is the "minimal" resolution needed to obtain a certain accuracy for a specific problem.Following the specifications in Lauritzen et al. (2012) we have calculated the minimal resolution as the resolution required to obtain an l 2 error norm for the cosine bell distribution that is less than 0.033.In the specifications (Lauritzen et al., 2012) the result should be obtained for an unlimited scheme, i.e., no shape-preserving limiter should be applied, e.g. on CSLAM, in this case.Minimal resolutions for CSLAM and HEL are listed in Table 5.
The minimal resolution for HEL is coarser than that for CSLAM, particularly for a maximum Courant number of 1.We therefore conclude that the effective resolution for HEL is higher than for CSLAM.
The minimal resolution for HEL is controlled by the strength of the mixing between parcels.If HEL is run without any parcel mixing the minimal resolution goes to infinity in the sense that the numerical solution becomes almost exactly equal to the analytic solution at any resolution, of course depending on the weight w 0 .

Filament diagnostics
The filament preservation diagnostic, l f , describes the transport scheme's ability to preserve thin filaments or gradients in the concentrations.l f is defined as where A(τ, t), the control volume, is a spherical area where the concentration is equal to or larger than a given threshold value τ .Test of convergence for linear advection on cubed sphere using error norms l 2 (left) and l ∞ (right).ncube refers to the number of grid cells in each direction on each of the 6 faces of the cubed sphere, i.e., ncube = 24 corresponds to an equatorial angular grid distance of 2π/(4 × 24) = 0.065 radians.Using the non-divergent flow, basic error norms l 2 and l ∞ for CSLAM-M and HEL, and two maximum Courant numbers (1.0 and 5.5), are listed in Tables 1, 2, and 3 for each of the three initial distributions.In general one can conclude that HEL is very accurate at low resolution for the more smooth distributions, i.e., the Gaussian and cosine hills, whereas CSLAM-M and HEL are comparable at the high Courant number and at the two finest resolutions.Another feature is that, as compared to CSLAM-M, HEL is less sensitive to the maximum Courant number.This is because HEL is influenced much less by the number of semi-Lagrangian re-mappings needed to finalize each simulation than CSLAM-M.
The convergence rates for each of the three initial distributions and for the two Courant numbers are listed in Table 4.In general CSLAM-M converges faster than HEL at high Courant number for the smooth distributions while the difference is small for the rough slotted cylinder distribution where the convergence rates in any case are low.The convergence rates are comparable in the low Courant number cases.

"Minimal" resolution
Numerical schemes may be constructed to converge fastat least for smooth distributions.However, as pointed out by Lauritzen et al. (2012) increases in resolution are often computationally expensive and therefore it is of interest to 785 identify some kind of measure of absolute accuracy for a numerical scheme.A diagnostic designed for this is the "minimal" resolution needed to obtain a certain accuracy for a specific problem.Following the specifications in Lauritzen et al. (2012) we have calculated the "minimal" resolution as the 790 resolution required to obtain an l 2 error norm for the cosine
The l f values for CSLAM-M and HEL are presented in Fig. 8 for simulations with maximum Courant numbers 1 and 5.5.As expected CSLAM-M is more diffusive; i.e., it maintains filaments less well, at maximum Courant number 1 as opposed to 5.5.This is because more remappings are required at the lower Courant number.The HEL values are calculated for both Eulerian and parcel, i.e., Lagrangian, representations.It can be concluded from Fig. 8 that at time t = T /2 of the simulation the Eulerian representation of HEL is generally more diffusive than CSLAM-M for the high maximum Courant number although the highest functional values are maintained to a higher degree than for CSLAM-M.This is because a relatively large weight, w 1 , in Eq. ( 13) is given to the provisional first-order Eulerian forecast, E ρn+1 , due to the highly inhomogeneously spaced parcels at this time.The Eulerian HEL representation does not change significantly with Courant number and is generally better at preserving the maximum values than the CSLAM-M.Therefore the low Courant number CSLAM-M is more diffusive than the corresponding Eulerian HEL representation.The Lagrangian l f values are generally much closer to 100 than the corresponding CSLAM-M and Eulerian HEL representations, and there is almost no dependency on the Courant number.The last observation is completely as expected since the parcel mixing takes place at the same time interval in the two cases.It is noted (not shown) that the l f values are quite insensitive to the mixing frequency between Lagrangian parcels.In a more general application, using nonanalytic trajectories, there could in theory be some spatial overlapping between parcels.However, such overlaps cannot be quantified due to deformation of the Lagrangian parcels into filaments, and HEL does not include information about their exact shape.The total area, however, is conserved.
It is important to note that without parcel mixing the Lagrangian l f values would all be exactly 100.As argued above, introduction of mixing is fundamental in all Lagrangian models in order to avoid long-term unphysical accumulation of energy at the smallest resolved scales.The same applies to such Eulerian-based models (in principle, including semi-Lagrangian models) where the inherent numerical mixing is "too weak" to properly represent non-linear scale interactions and prevent spectral blocking for a given model resolution.This is typically the case in e.g.spectral or
pseudo-spectral models.So, although we know that explicit mixing must be introduced in some undiffusive models for purely physical reasons, we do not know exactly how much mixing is required.I.e., the optimal l f values are, unfortunately, also unknown.The fundamental idea behind the parcel mixing applied here has been to base it on simple geometric considerations and thereby obtain a simple first principle guess on the required amount of actual physical mixing.

Pre-existing functional relations and mixing
To evaluate the mixing properties discussed in Sect.1.2 the statistics proposed in Lauritzen and Thuburn (2012) and Lauritzen et al. ( 2012) have been calculated for initial conditions consisting of two tracers: cosine bells -corresponding to χ in Fig. 1 -and corresponding non-linearly related bells -corresponding to ξ .The flow is the same non-divergent deformation flow as above.
These mixing statistics include real mixing, l r , rangepreserving unmixing, l u , and overshooting, l o .The more precise definitions of l r , l u , and l o are provided in Appendix B.
The error norm, l o , should always be zero, indicating that the scheme in question is shape preserving.However, the second norm, l u , which ideally should be zero as well, will generally not be zero, unless the scheme is semi-linear and monotone (Thuburn and McIntyre, 1997).This was one of the motivational factors for the development of HEL.The first norm, l r , should be a non-zero value, since real mixing is always present; it should, however, as described in Sect.3, not be artificial numerical mixing but physically based mixing.
The mixing diagnostics for CSLAM-M and HEL are listed in Table 6.Mixing diagnostics are important indicators for the influence of transport schemes on chemical reaction rates and equilibria.Since chemistry calculations (not dealt with here) are performed in Lagrangian space only the parcel values of l r , l u , and l o are listed for HEL.It can be seen from Table 6 that HEL behaves according to its construction: only real mixing and no range-preserving unmixing or overshooting takes place.In CSLAM-M a weak range-preserving unmixing can be seen, and the general level of real mixing is larger than for HEL, particularly at low resolution and low Courant number. .8. Diagnostics for filament preservation, l f , for the two resolutions 1.5 • (left) and 0.75 • (right).The curves show the shape preserv LAM-M (crosses), HEL in Eulerian grid (squares), and the HEL in Lagrangian space (circles).Each scheme is shown with a maxim urant number of 1.0 (dotted line) and 5.5 (solid line), respectively.
Deformational passive advection with divergence aditional l 2 and l ∞ error norms for the strongly divergent w are presented in Table 7 for the cosine hill initial condin.In this case the two maximum Courant numbers tested 0.6 and 3.2, respectively, and it can be seen that HEL is w considerably more accurate than CSLAM-M, particuly for the small Courant number tested.As for the nonergent case HEL is relatively insensitive to the Courant mber.
Implementation and tests in a shallow water dynamical model this section we demonstrate a dynamical application of L namely as transport scheme also for the dynamical core an existing geophysical type model, namely a simple shalwater model in Cartesian geometry (Kaas, 2008).Dete the simple geometry applied in this model it does inde the fundamental processes and thereby potential probs in the shallow water framework.
Model setup e governing differential equations are is the geopotential thickness of the flow, and h s the stationa surface geopotential height ("topography").h m is a quant we can think of as the contribution to geopotential heig from the m'th tracer.The mixing ratio for the m'th tra can be evaluated from i.e., formally the volume density for this tracer is ρ m ρh m /h, where ρ is the density of the "dry" fluid in the sh low water model.It has been assumed that h m h sin 925 otherwise, the expression in Eq. ( 39), which formally re resents specific density would not approximate mixing tio.The term F h represents a weak globally mass conse ing Newtonian relaxation towards the initial "zonal" avera profile of h.F h mimics the effect of diabatic processes.

930
nally the C, D, and S terms represents possible chemist diffusion/mixing, and possible source and sink terms (i emissions, depositions, sedimentations), respectively.Although u, v, and h do not depend on the h m values the tracers the model is formally set up as "online couple 935 (see, e.g., Grell and Baklanov, 2011), i.e., all equations solved each time step.
The integration domain covers a domain defined by x [x min ,x max ] and y ∈ [y min ,y max ], with x min = 0 km, x max 20000 km, y min = −10000 km, and y max = 10000 km.T boundary conditions are periodic in both directions and w enforced symmetry around the line y = 0 ("Equator") the variables u, h, h s , F h , and h m ,m = 1,...,M , and an symmetry around the same line for v. Also, the Coriolis p rameter y

Deformational passive advection with divergence
Traditional l 2 and l ∞ error norms for the strongly divergent flow are presented in Table 7 for the cosine hill initial condition.In this case the two maximum Courant numbers tested are 0.6 and 3.2, respectively, and it can be seen that HEL is now considerably more accurate than CSLAM-M, particularly for the small Courant number tested.As for the nondivergent case HEL is relatively insensitive to the Courant number.

Implementation and tests in a shallow water dynamical model
In this section we demonstrate a dynamical application of HEL namely as a transport scheme also for the dynamical core of an existing geophysical type model, namely a simple shallow water model in Cartesian geometry (Kaas, 2008).
Despite the simple geometry applied in this model it does include the fundamental processes and thereby potential problems in the shallow water framework.

Model setup
The governing differential equations are where u, v are the flow speed components in the x-y plane, f is the Coriolis parameter, g the gravitational acceleration, h is the geopotential thickness of the flow, and h s the stationary surface geopotential height ("topography").h m is a quantity we can think of as the contribution to geopotential height from the m'th tracer.The mixing ratio for the m'th tracer can be evaluated from i.e., formally the volume density for this tracer is ρ m = ρh m / h, where ρ is the density of the dry fluid in the shallow water model.It has been assumed that h m h since, otherwise, the expression in Eq. ( 39), which formally represents specific density, would not approximate mixing ratio.The term F h represents a weak globally mass-conserving Newtonian relaxation towards the initial "zonal" average profile of h.F h mimics the effect of diabatic processes.Finally the C, D, and S terms represent possible chemistry, diffusion/mixing, and possible source and sink terms (i.e., emissions, depositions, sedimentations), respectively.Although u, v, and h do not depend on the h m values of the tracers, the model is formally set up as "online coupled" (see, e.g., Grell and Baklanov, 2011); i.e., all equations are solved each time step.

Geosci
The integration domain covers a domain defined by x ∈ [x min , x max ] and y ∈ [y min , y max ], with x min = 0 km, x max = 20 000 km, y min = −10 000 km, and y max = 10 000 km.The boundary conditions are periodic in both directions and with enforced symmetry around the line y = 0 ("Equator") for the variables u, h, h s , F h , and h m , m = 1, . . ., M, and antisymmetry around the same line for v. Also, the Coriolis parameter f = 2 sin π y y max − y min is anti-symmetric around y = 0 ( is the angular velocity of the Earth rotation).
As for the inert tracer applications tested above the strategy followed for applying HEL in the shallow water model is to use parcel densities/geopotential thicknesses to modify an existing solution in Eulerian space.As the underlying solution we have used a locally mass-conserving, semi-Lagrangian transport scheme LMCSL (i.e., not the CSLAM as above) with a semi-implicit treatment of the gravityinertial wave terms, in combination with an Arakawa C-grid staggering; see Kaas (2008) for details.
In the present implementation we have only applied the HEL technique on the mass fields (h and the h m 's) of the model.The wind field forecast is based on the same thirdorder semi-Lagrangian scheme as in Kaas (2008).
As for the case of passive/inert advection the divergent expansion/contraction factors, σ n+1/2 , that are needed to include the effects of divergence in Lagrangian space are first calculated in Eulerian space: where E h is the complete provisional Eulerian space forecast including semi-implicit adjustments, and E , advh is the corresponding purely advective semi-Lagrangian forecast, i.e., also an Eulerian space forecast.The Eulerian space values E σ n+1/2 are now interpolated to each parcel location (see Sect. 2.4) at time level n + 1 and subsequently multiplied on the parcel values of L h and L h m , m = 1, . . ., M. Mixing depending on the flow deformation rate is then performed in Lagrangian space as described in Sect. 3 The final Eulerian space forecasts of E h and E h m , m = 1, . . ., M are obtained via exactly the same nudging procedure as described above in Sect.2.6 For the inert and passive advection tests with prescribed analytical velocities in Sect. 4 the underlying Eulerian space forecasts were all based on a first-order numerical scheme providing good numerical efficiency.Here in our dynamical model implementation we have found that it is necessary to keep third-order accurate remappings in the semi-Lagrangian scheme in order to obtain sufficiently accurate estimates of pressure gradient terms and to ensure sufficiently accurate coupling between the momentum equations (Eqs.35 and 36), and the continuity equation (Eq.37).
Although we have used third-order remappings it could of course still be possible to run with first-order remappings for all the tracers, m = 1, . . ., M,i.e.,when solving Eq. (38).In this way one can retain the same high multi-tracer efficiency as in the transport tests in Sect. 4. One may argue, though, that this violates the mass-wind consistency property, i.e., the 7th of the desired properties listed in Sect. 1.It is noted, however, that since the bulk of the model memory in HEL is kept in Lagrangian space, this is now only a theoretical problem: if an inert and passive tracer m is initialized as h m = h, the Lagrangian space density of this variable will continue to be exactly equal to that of h; i.e., at any time step n we have L h n m = L h n , independent of the order of accuracy of the underlying Eulerian-space-based forecast of E h m .

Estimation of trajectories
For the passive advection tests in Sect. 4 all trajectories were calculated from analytically defined velocity fields.In the shallow water model, as well as in any update to a threedimensional model, it is necessary to estimate trajectories.

Discussion
A number of issues related to the introduction of the new HEL scheme deserve some discussion provided in the following subsections.

Number of parcels
In the present applications of HEL the total number of parcels is equal to the number of Eulerian grid cells.In principle one could, however, easily increase the number of parcels, al-ber of parcels in the shallow water model is four times that 1130 of the number of Eulerian grid cells has been carried out, and the results were very similar to those in the lower panels of Fig. 10 and the middle panels of Fig. 11, although the tracer mixing ratio fields in these additional simulations, as expected, are slightly smoother.1135

Mixing
As noted previously, the mixing between Lagrangian parcels introduced here is based on simple geometrical principles obtained with a semi-Lagrangian model, as here, upstream trajectories are needed; i.e., one has to identify all departure points at time n t for trajectories which at time (n + 1) t end up in each Eulerian grid cell centroid.In HEL it is furthermore required to calculate downstream arrival points at time (n + 1) t for trajectories beginning at the irregularly spaced locations for each Lagrangian parcel at time n t.
For the upstream trajectories we use the approach described in Kaas (2008).This is a conventional iterative procedure using two iterations.However, here we use third-order accurate bi-cubic Lagrange interpolations of the upstream velocities at time level n as opposed to the first-order interpolations in Kaas (2008).As explained in Kaas (2008), in order to ensure satisfactory behavior of the LMCSL scheme, it necessary to include the effect of accelerations in the estimate of the trajectories.For the estimation of downstream parcel trajectories an equivalent procedure has been followed (see Appendix C for details).

Shallow water model results
The initial state of the shallow water model is chosen rather arbitrarily as wavy structure in the geopotential surface height field, h, shown in Fig. 9, with the velocity field (not shown) simply initialized to be in geostrophic balance with this mass field.Figure 9 also shows the bottom topography, h s , consisting of a "sharp" isolated "mountain", and an initial field of mixing ratio for a single tracer field, which here simply is a step function in a background of zero mixing ratio.
The results obtained with HEL are compared to those obtained with the third-order semi-implicit LMCSL scheme (Kaas, 2008) without introduction of any shape-preserving filters.It is noted (not shown) that one could just as well have verified HEL against a traditional semi-implicit semi-Lagrangian (SISL) time stepping scheme, since SI-LMCSL and SISL turn out to produce almost identical forecasts.All simulations presented below have been performed with a horizontal resolution of 128 × 128 points cells −1 .The time step was one hour, which gives a maximum Courant number slightly below one.The reader is informed that LMCSL and HEL (both with and without parcel mixing) can be run stably with maximum Courant numbers up to about 3, and thereafter numerical mountain wave resonance (Rivest et al., 1994;Lindberg and Alexeev, 2000) becomes visibly destabilizing.No off-centering or other techniques to control mountain wave resonance were introduced.
Figure 10 shows the surface height field and the mixing ratio field after 48 h of simulation for the initial conditions plotted in Fig. 9.It can be seen that the geopotential height fields for each of the two forecasts are almost indistinguishable, although the HEL result is based on densities at the locations of the irregularly spaced Lagrangian parcels.The mixing ratio fields are also similar and it can be seen that HEL via its design is shape preserving.
Corresponding height and mixing ratios after 20 days of simulations are shown in Fig. 11.It can be seen that HEL and LMCSL continue to produce very similar results for the height field despite the underlying non-linear model equations (non-linear chaotic error growth does not become visible until around day 30-40).The tracer mixing ratios are also similar although the LMCSL field is smoother than that obtained with HEL.In the lower panel of Fig. 11 we have also -for illustrative purposes -plotted the height and mixing ratio field obtained in an additional HEL simulation where the parcel mixing was switched off.The mixing ratio for the passive inert tracer is now highly unrealistic, demonstrating the importance of the mixing we have introduced.In this unmixed version the initial step function type mixing ratio has been "cut into bits and pieces" determined by the location of the individual Lagrangian parcels.It is noted (not shown) that the density of parcels is quite homogenous in the sense that there are no larger regions without any parcels.Surprisingly, the height field continues to be smooth and very similar to that of LMCSL and of the parcel mixed version of HEL, and the unmixed HEL continues to be numerically stable.We suspect that dynamical adjustments between the velocity and mass fields in the model tend to prevent development of local small-scale "lows" and "highs", which one could envisage due to errors in parcel trajectory calculations.As an example, the flow around a parcel with anomalously low density, i.e., height, as compared to its nearest neighbors, tends to be unbalanced in a way leading to development of local convergence, which eventually increases the parcel density via the σ factor (see Eq. 40).
Admittedly, the shallow model we have used here is based on a simplified cartesian geometry.It is planned to implement HEL in a shallow model in spherical geometry, and to apply it for standard validation tests of shallow water model dynamics.However, since we have here obtained almost indistinguishable results in the HEL version and the underlying purely Eulerian-based model (left panels of Figs. 10 and 11), it is likely that the result obtained in spherical geometry will also be very similar to the underlying Eulerian-based scheme -whatever that might be.
By introducing a passive tracer with initial density equal to that of h, i.e., the mixing ratio is equal to one, we have tested to what extent the HEL is "wind-mass" consistent, and it has been found that this is indeed the case since the density field (not shown) continued to be almost completely identical to that of h even in long simulations.It is noted that for this test it was necessary to switch off the Newtonian cooling F h since otherwise one field would have been forced while the other was not.

Pre-existing functional relations in the shallow water model
It was shown in Sect.4.2.3 that the mixing diagnostics l r , l u , and l 0 obtained with HEL in strongly deformational flow were quite acceptable.We have performed a corresponding calculation for transport in the shallow water flow using exactly the same initial functional shapes as in Sect.4.2.3, and in Fig. 12 these are shown as time series for the nonshape-preserving LMCSL scheme and for HEL.As expected LMCSL produce non-zero values of both range-preserving unmixing, l u , and of overshooting, l 0 .For the Lagrangian representation in HEL l u and l 0 are both zero via construction, while the same is the case for l 0 in the Eulerian HEL representation.The level of l u in the latter is very small as compared to that in LMCSL.  is about half of that.The fact that the overall levels of real mixing in HEL and LMCSL have the same order of magnitude indicates that the ways we have obtained the mixing and Eulerian space smoothing in HEL are not completely unrealistic.

Discussion
A number of issues related to the introduction of the new HEL scheme deserve some discussion provided in the following subsections.

Number of parcels
In the present applications of HEL the total number of parcels is equal to the number of Eulerian grid cells.In principle one could, however, easily increase the number of parcels, although this has a computational cost, particularly due to increased costs of mixing.A test (not shown) where the number of parcels in the shallow water model is four times that of the number of Eulerian grid cells has been carried out, and the results were very similar to those in the lower panels of Fig. 10 20 Eig mixing in other ways, which in practice can lead to a stronger or weaker mixing between parcels.
One potentially controversial issue is the degree of directional bias of the mixing.As described in Sect. 3 our mixing is biased so that it is dominated by mixing with neighbors 1145 that are aligned along the asymptotic dilatation axis.This

Mixing
As noted previously, the mixing between Lagrangian parcels introduced here is based on simple geometrical principles minimizing the need for empirically based tuning.However, it is of course possible to formulate such physically based mixing in other ways, which in practice can lead to a stronger or weaker mixing between parcels.
One potentially controversial issue is the degree of directional bias of the mixing.As described in Sect. 3 our mixing is biased so that it is dominated by mixing with neighbors that are aligned along the asymptotic dilatation axis.This approach is based on the geometrical principle illustrated in Fig. 4: remapping of the parcels into regular squared shapes filling the integration domain only requires mixing along the asymptotic dilation axis.This remapping is needed to obtain an unaliased representation of the Lagrangian parcel densities -a problem that is quite different from that of molecular mixing, which is generally isotropic in nature.We have tested the effect of performing the mixing with a fake directional bias, which is not along the asymptotic dilation axis but instead along an axis perpendicular to it.The result (not shown) is an excessive damping and considerably larger error norms for all the inert passive transport tests reported above.Our actual choice of κ = 10 was based on a compromise: a much smaller value would be too isotropic and too damping, and a much larger value would result in too little realized mixing; i.e., the parcel deformations L δ would grow to unrealistic values in strongly deformational flows.
The present paper does not investigate the influence of parcel mixing on the distribution of energy on different wave numbers.This is the subject of ongoing research.

Local versus global nudging towards the target values
The fourth desirable property listed in Sect.1.1 states that a transport scheme should be transportive and local.While this is fully achieved in the parcel (Lagrangian) representation of HEL, the locality property is formally not fulfilled in the Eulerian representation since the nudging we perform on the Eulerian-based forecast towards the target values is performed globally (see Sect. 2.6).We have formulated and tested a local nudging which only involves mass reorganizations between neighboring cells.This locally massconserving version of HEL performs satisfactorily with results (not shown) that are almost comparable or somewhat degraded as compared to the standard HEL version.However, the local nudging is quite expensive from a numerical point of view, and therefore this version has not been investigated further.In practice since the bulk memory in HEL is in Lagrangian space, the standard version of HEL is highly local, as can also be seen directly from all the plots presented above.In fact, since the limits within which local Eulerian values can change due to the nudging are defined from Eigil Kaas: A hybrid Eulerian Lagrangian numerical scheme.form on the Eulerian based forecast towards the target values is performed globally (see Sect. 2.6).We have formu-1175 lated and tested a local nudging which only involves mass reorganizations between neighboring cells.This locally mass conserving version of HEL performs satisfactorily with results (not shown) that are almost comparable or somewhat degraded as compared to the standard HEL version.How-1180 ever, the local nudging is quite expensive from a numerical point of view, and therefore this version has not been investigated further.In practice since the bulk memory in HEL is in Lagrangian space the standard version of HEL is highly local as can also be seen directly from all the plots presented 1185 above.In fact, since the limits within which local Eulerian values can change due to the nudging are defined from locally defined values, the global nudging can be considered a localized process.

Computational efficiency 1190
The purpose of the present paper has been to describe HEL, and to demonstrate its accuracy.A careful investigation of the computational cost of HEL would require multiple tests on a massively parallel computer system.Here we have only performed single processor CPU tests using an Intel 1195 Core2 Duo, E6550 @ 2.9 GHz processor, and the Intel Fortran 13.0.0compiler with flags: -ipo, -O3, -no-prec-div,static, and -xHost.The tests reveal that HEL is considerably faster than CSLAM-M for the passive tracer test presented in Sect.4, particularly when many tracers ara considered.As an 1200 example Fig. 13 shows the CPU timing required to perform the non-divergent deformation test in Sect.4.2 with an equatorial spatial resolution of 0.75 • , and a maximum Courant number of 5.5.The number of passive inert tracers was varied from 2 to 20 and it is seen that the multi-tracer efficiency 1205 of of HEL is better than that of CSLAM-M.
A main reason for the faster performance of HEL is that only first order accurate re-mappings are needed in the underlying Eulerian forecast for passive transport, while in CSLAM-M third order re-mapping have been used.locally defined values, the global nudging can be considered a localized process.

Computational efficiency
The purpose of the present paper has been to describe HEL, and to demonstrate its accuracy.A careful investigation of the computational cost of HEL would require multiple tests on a massively parallel computer system.Here we have only performed single processor CPU tests using an Intel Core2 Duo, E6550 @ 2.9 GHz processor, and the Intel Fortran 13.0.0compiler with flags -ipo, -O3, -no-prec-div,static, and -xHost.The tests reveal that HEL is considerably faster than CSLAM-M for the passive tracer test presented in Sect.4, particularly when many tracers ara considered.As an example, Fig. 13 shows the CPU timing required to perform the non-divergent deformation test in Sect.4.2 with an equatorial spatial resolution of 0.75 • , and a maximum Courant number of 5.5.The number of passive inert tracers was varied from 2 to 20, and it is seen that the multi-tracer efficiency of HEL is better than that of CSLAM-M.
A main reason for the faster performance of HEL is that only first-order accurate remappings are needed in the underlying Eulerian forecast for passive transport, while in CSLAM-M third-order remapping have been used.
Although we have not performed parallel efficiency tests, the code has been prepared somewhat for parallelization.The most important issue relates to the way individual parcels are transferred from the memory of one CPU/node to another.As usual in geophysical fluid dynamics each CPU is reserved for a certain number of horizontal Eulerian grid cells including a halo zone.Corresponding to this, the actual physical location of each Lagrangian parcel determines in which part of the memory it is stored.Since the number of parcels in a given domain can vary significantly due to the divergence of the flow, the actual memory allocation for Lagrangian parcel information required for each CPU must be somewhat higher than that corresponding to the average parcel density.
Although our first tests suggest that HEL is computationally efficient, particularly for multiple tracers, there is an important memory penalty.Considering, e. passive transport using a traditional Eulerian-based scheme with K Eulerian grid cells/points, and with M different tracers, the total number of prognostic variables is K × M.
For HEL, however, the corresponding number is typically K × (M + M + 6), where "6" reflects the additional variables of parcel deformation ( L δ), parcel volume ( L V ), components of the parcel position vector ( L x and L y) and the position components of the auxiliary parcel ( L x a and L y a ).

Application in a three-dimensional model
Arguably the most important issue is how the HEL scheme is extended for use in three-dimensional applications.At least two fundamentally different options may be investigated: -using a Lagrangian or quasi-Lagrangian vertical coordinate as in Sørensen et al. (2013) each Lagrangian parcel could stay within the same layer and have an instantaneous vertical extension equal to that of the layer.The vertical extension of Lagrangian parcels may in fact be used as an additional prognostic variable and used to determine the layer thickness just as it was demonstrated with the shallow water model in Sect. 5 (a one layer model).
-letting parcels float freely between the Eulerian vertical levels as in e.g.ATTILA (Stenke et al., 2008(Stenke et al., , 2009)).For a fully dynamic model implementation instantaneous parcel densities and temperatures may be used to modify or nudge the parcel vertical positions so that they are more consistent with the vertical structure of the local atmosphere as represented in the Eulerian grid.It is speculated that such nudging may be used as tool to stabilize the model with regard to fast gravity and sound waves.
In a three-dimensional application the mixing between parcels must be reconsidered, and presumably it is necessary to separate horizontal and vertical mixing since they represent quite different processes in stratified fluids.

Prognostic variables in a dynamic model
In the dynamic tests in Sect. 5 the HEL scheme was only used for transport of the mass field.In a future application in a baroclinic model in spherical geometry one may obviously use HEL as a transport scheme for other quantities and invariants such as momentum, total energy, and potential vorticity.

Using HEL in an atmospheric chemical transport model
HEL should be suited quite well as a numerical fundament for an atmospheric chemical transport model, and has in fact been designed with this goal in mind.Since the mixing we have introduced between neighboring parcels is purely real, it is logical to perform chemical calculations in the parcel space.In an accompanying paper (Hansen et al., 2012) it is described and demonstrated how HEL performs when it is used as the underlying numerical scheme in a simple socalled online atmospheric chemistry transport model.

Conclusions
The original motivation for developing HEL was to set up a numerical method for use in fluid dynamics that fulfills all the 11 desirable properties listed in Sect.1.1.The passive transport tests described in Sect.4, the testing in the dynamical model (Sect.5), and the efficiency investigations in Sect.6.4 clearly show that HEL fulfils all 11 properties, and that it is very multi-tracer efficient.At high Courant numbers, the convergence rates of HEL are lower than those of CSLAM-M, with which it has been compared.However, the absolute level of accuracy in HEL is very high.
A fundamental component of HEL is a directionally biased mixing between neighboring cells, which is proportional to the deformation rate of the flow.This mixing has been formulated in such a way that the 11th property (avoidance of spurious numerical mixing/unmixing) continues to be fulfilled in the Lagrangian representation of HEL.Thus HEL should be ideal as the underlying scheme for chemical transport in the atmosphere particularly if chemistry calculations are performed in Lagrangian space.

List of prognostic variables
The set of prognostic variables for solving the twodimensional transport/continuity problem (passive transport) in HEL includes

Fig. 1 .
Fig. 1.Illustration of numerical mixing categories.The thick curve is the pre-existing functional relation between tracer ξ and tracer χ.Any new relative concentrations (ξ k ,χ k ), generated by the transport scheme, can be represented as a point.If the point falls within the shaded convex hull, it is classified as real mixing.If within the dashed rectangle but outside the shaded area it is classified as range preserving unmixing, and, finally, if outside the dashed rectangle it is classified as overshooting (adopted fromLauritzen and Thuburn, 2012).

Fig. 1 .
Fig.1.Illustration of numerical mixing categories.The thick curve is the pre-existing functional relation between tracer ξ and tracer χ.Any new relative concentrations (ξ k , χ k ), generated by the transport scheme, can be represented as a point.If the point falls within the shaded convex hull, it is classified as real mixing; if within the dashed rectangle but outside the shaded area it is classified as rangepreserving unmixing; and, finally, if outside the dashed rectangle it is classified as overshooting (adopted fromLauritzen and Thuburn, 2012).
and the total mass of the target values, M T = K k=1 E T k • E V k , as well as the discrepancy, M = M T − M, between the two.This will lead to three different possibilities: M < 0: Mass of target field is too small.M = 0: Mass of target field is correct.

Fig. 3 .
Fig.3.Downstream advection of two different Lagrangian parcels and their auxiliary parcels in a pure shear flow.For each parcel P and P a denotes, respectively, its position and the position of the auxiliary parcel.The associated auxiliary vectors are marked black and gray, respectively, to distinguish between the two parcels.In this example the (true) asymptotic / Lagrangian dilatation axis, indicated with dashed lines, is the same all over the small domain shown, although its direction rotates clockwise from time step n to time step n + 1.Note that it is the relative movements between parcels and auxiliary parcels that are relevant.After several timesteps the auxiliary vectors will align approximately with the true asymptotic dilatation axis.

Fig. 4 .
Fig.4.This figure describes why we use a factor of 0.25 between the deformation within one time step ∆t • D and the fraction of the parcel area ∆ 2 that should be mixed.The underlying assumption is that the deformation is mainly due to shear, and the figure shows how cells are then deformed within one time step from regular squares into irregular parallelograms with the same area (assuming zero divergence).In order to re-shape a cell into squares the mass in the two shaded triangles must be given off via interpolations to the neighboring cells.The factor comes because the sum of the triangular areas is 0.25∆t • D∆ 2 .
645the problem of numerical mixing and unmixing is eliminated.It is noted that the introduction of a flow-dependent mixing based on the degree of deformation is not new.Sadourny and Maynard (1997) introduced a horizontal diffusion which was dependent on the magnitude of the deformation rate of 650 the flow, and laterVáňa et al. (2008) used a similar approach to obtain a flow dependent degree of mixing in a semi-Lagrangian model.Also ATTILA(Stenke et al., 2009) and CLaMS (Chemical Lagrangian Model of the Stratosphere)(McKenna et al., 2002) employ mixing depending on hori-655 zontal flow deformation rate.

Fig. 4 .
Fig.4.This figure describes why we use a factor of 0.25 between the deformation within one time step t • D and the fraction of the parcel area 2 that should be mixed.The underlying assumption is that the deformation is mainly due to shear, and the figure shows how cells are then deformed within one time step from regular squares into irregular parallelograms with the same area (assuming zero divergence).In order to re-shape a cell into squares the mass in the two shaded triangles must be given off via interpolations to the neighboring cells.The factor comes because the sum of the triangular areas is 0.25 t • D 2 .

Geosci
Fig.5.Test of convergence for linear advection on cubed sphere using error norms l 2 (left) and l ∞ (right).ncube refers to the number of grid cells in each direction on each of the 6 faces of the cubed sphere, i.e., ncube = 24 corresponds to an equatorial angular grid distance of 2π/(4 × 24) = 0.065 radians.
igil Kaas: A hybrid Eulerian Lagrangian numerical scheme.13 ig.7. Simulated distributions at times t = T /2 (left panels) and t = T (right panels), based on the slotted cylinder initial distributio om top to bottom the plots show results obtained with CSLAM-M, HEL and HEL without any parcel mixing, respectively, all run wi maximum Courant number of 5.5 and an equatorial resolution of 1.5 • .ble 1. Statistics for the Gaussian hill problem.The columns show maximum Courant number (C), Equatorial resolution in degrees (∆λ d the error norms l 2 and l ∞ , respectively, for both CSLAM-M and HEL.

Fig. 7 .
Fig.7.Simulated distributions at times t = T /2 (left panels) and t = T (right panels), based on the slotted cylinder initial distribution.From top to bottom the plots show results obtained with CSLAM-M, HEL and HEL without any parcel mixing, respectively, all run with a maximum Courant number of 5.5 and an equatorial resolution of 1.5 • .

Fig. 8 .
Fig.8.Diagnostics for filament preservation, l f , for the two resolutions 1.5 • (left) and 0.75 • (right).The curves show the shape-preserving CSLAM-M (crosses), HEL in Eulerian grid (squares), and the HEL in Lagrangian space (circles).Each scheme is shown with a maximum Courant number of 1.0 (dotted line) and 5.5 (solid line), respectively.

Fig. 9 .
Fig. 9. Initial height field h (upper left) and surface topography (upper right) for the shallow water model.The lower panel shows an example of an initial field of a single inert tracer.Only the "Northern Hemisphere" part of the fields are plotted.

Fig. 10 .
Fig.10.Results after 48 h of simulation based on the initial conditions and bottom topography in Fig.9.Panels to the left show the surface height field, while those to the right show the mixing ratio.The upper panels show the result obtained with the LMCSL scheme, while HEL is shown below.

Fig. 9 .
Fig. 9. Initial height field h (upper left) and surface topography (upper right) for the shallow water model.The lower panel shows an example of an initial field of a single inert tracer.Only the "Northern Hemisphere" part of the fields are plotted.

Fig. 10 .
Fig. 10.Results after 48 h of simulation based on the initial conditions and bottom topography in Fig. 9. Panels to the left show the surface height field, while those to the right show the mixing ratio.The upper panels show the result obtained with the LMCSL scheme, while HEL is shown below.

Fig. 11 .
Fig. 11.As Fig. 10 but after 480 h of simulation.The lower panels show results in a HEL configuration without any parcel mixing.

Fig. 12 .
Fig.12.Time series of mixing diagnostics l r , l u , and l 0 obtained in 48 h simulations with the shallow water model using the LMCSL and HEL schemes.For HEL both the Eulerian and the Lagrangian representations are plotted.It is noted that in the Eulerian representation of HEL l 0 = 0 (not shown), while in the Lagrangian representation both l u = 0 and l 0 = 0 (not shown).
g., two-dimensional www.geosci-model-dev.net/6 Density of dry air in Eulerian space E ρ m , m = 1, . . ., M Density of tracers in Eulerian space L ρ Density of dry air in Lagrangian parcelsL ρ m , m = 1, . . ., M Density of tracers in Lagrangian parcels variable is represented by K Eulerian cells or P Lagrangian parcels.Note that for passive transport it has here been assumed that velocity components are known (nonprognostic) variables.
Stenke et al. (2009)ed in a general circulation model to simulate transport of water vapor and cloud water.This was extended inStenke et al. (2009)to an interactively coupled chemistry-climate model version of ATTILA; i.e., water vapor, cloud water, and mixing ratios of all chemical tracers are known for each Lagrangian parcel.ATTILA is able to maintain steep gradients, is mass conserving, numerically non-diffusive, and has been used e.g. for studying the climate radiative forcing related to aircraft contrails www.geosci-model-dev.net/6/2023/2013/Geosci.Model Dev., 6, 2023-2047, 2013 et al. (2008),

Table 1 .
Statistics for the Gaussian hill problem.The columns show maximum Courant number (C), equatorial resolution in degrees ( λ), and the error norms l 2 and l ∞ , respectively, for both CSLAM-M and HEL.

Table 2 .
As Table1but for the slotted cylinder problem.

Table 3 .
As Table1but for the cosine hill problem.

Table 7 .
As Table 1 but for the cosine hill problem in strongly divergent flow.Note that the maximum Courant numbers C in this case are different from those of Table 1.