Lagrangian advection scheme with shape matrix ( LASM ) v 0 . 2 : interparcel mixing , physics-dynamics coupling and 3-D extension

Abstract. The interparcel mixing algorithm in the Lagrangian advection scheme with shape matrix (LASM) is updated to make the scheme more robust. The linear degeneration criterion is replaced by the maximum deviation of the skeleton points so that the new algorithm is more effective in controlling the shape of parcels, which is vital for long time simulation. LASM is inherently shape-preserving without any complicated filter or limiter, and it is linear. This fact contributes to the ability to preserve the sum of multiple tracers exactly on the parcels in LASM. A newly proposed terminator "toy"-chemistry test is used to test LASM, which shows that LASM can preserve the weighted sum of two reactive species precisely. The physics–dynamics coupling (i.e., tendency evaluation type) is also discussed. A flow generated by a WRF large-eddy simulation is also used to test the 3-D extension of LASM.


Introduction
Lagrangian modeling approaches have long been recognized as a better alternative to the Eulerian or semi-Lagrangian ones for modeling the advection of tracers, due to their ability to reduce the numerical diffusion significantly, though some roadblocks need to be cleared.Therefore, more and more Lagrangian advection (or transport) schemes were proposed during the last two decades from different research communities.Based on the application objects, the research directions of the Lagrangian methods can be divided into four categories: (1) atmospheric dispersion modeling, (2) atmospheric chemical transport modeling, (3) cloud-resolving modeling, and (4) general circulation modeling.Lin et al. (2013) gave a panoramic view of the Lagrangian modeling of the atmosphere.
In the first category, the turbulent diffusion needs to be tackled thoroughly, because the temporal and spatial scales are relatively short and small respectively.These models are called Lagrangian particle dispersion models (LPDMs).The flow is decomposed into resolved mean and unresolved turbulent parts (the molecular diffusion is regularly ignored), where the first part comes from the reanalysis or model output, and the latter part is described by the stochastic process (Thomson, 1987).Two of the representative schemes (or models) are HYSPLIT (Draxler and Hess, 1998) and FLEX-PART (Stohl et al., 2005), which are used widely in the atmospheric emergency response, such as the volcanic ash cloud simulation.It is noteworthy that their turbulent diffusion formulations are different.Although the turbulent diffusions are both formulated in terms of the turbulent velocity components, FLEXPART solves the Langevin equation with the assumption of a Markov process, but HYSPLIT uses an autocorrelation coefficient and a random component.Recently, FLEXPART was extended to WRF in Brioude et al. (2013).Another model is NAME III developed in UK Met Office (Jones et al., 2007).
The chemical transport models (CTMs) consider a large number of tracer species and chemical reactions in the second category.One of the widely used Lagrangian CTMs is STOCHEM (Collins et al., 1997;Stevenson et al., 1997Stevenson et al., , 1998;;Utembe et al., 2010).Chemical Lagrangian Model of the Stratosphere (CLaMS) (McKenna et al., 2002b, a) is a more sophisticated CTM with a featured flowdeformation-based interparcel mixing.Hoppe et al. (2014) implemented CLaMS into ECHAM/MESSy Atmospheric Chemistry model (EMAC) and compared it with the standard flux-form semi-Lagrangian (FFSL) transport scheme in EMAC.In their results, the transport barriers at the edge of the polar vortex are better simulated in CLaMS.Another Lagrangian CTM is CiTTyCAT (Pugh et al., 2012).The boundaries between the above two categories are actually blurred; for example STOCHEM is incorporated in NAME III.
The third category includes the works on the cloud simulation in Lagrangian way.The clouds are divided into warm clouds and ice clouds.Shima et al. (2009) proposed the socalled super-droplet method (SDM) for simulation of warm clouds.Each super-droplet represents a multiple number of real droplets with the same attributes and position.Sölch and Kärcher (2010) introduced a novel Lagrangian cirrus module (LCM), where the ice phase is simulated by a large number of simulation ice particles (SIPs) and mechanisms for an adaptive control of the SIP number are included (Unterstrasser and Sölch, 2014).
The fourth research effort comes from the general circulation modeling (GCM) community.The current research topic is the Lagrangian advection scheme in the dynamical core, which is to replace the Eulerian or semi-Lagrangian schemes by the Lagrangian ones.The turbulent diffusion is not considered, which is left to other physical parameterizations.The pioneer works were done by Reithmeier and Sausen (2002), which proposed Atmospheric Tracer Transport In a LAgrangian model (ATTILA).ATTILA was incorporated in ECHAM4 atmospheric GCM (AGCM) in Stenke et al. (2008), and positive results were gained compared with the standard semi-Lagrangian scheme, because ATTILA is able to maintain steeper and more realistic gradients.In Kaas et al. (2013), a hybrid Eulerian-Lagrangian scheme (HEL) was proposed, which has some similarities to particle-in-cell (PIC) methods.In HEL, one Eulerian scheme runs in parallel with the Lagrangian scheme, and their results are merged each time step.Lagrangian advection scheme with shape matrix (LASM), which was proposed in Dong et al. (2014) and is the research subject of this study, also belongs to this category.LASM conserves the tracer total mass strictly in the sense of parcels, which is similar to ATTILA.Stenke et al. (2009) used ATTILA in a fully coupled chemistry-climate model (CCM), so the fourth category is also blurred with the second one, and may be even with the first one in the future to form a unified scheme.
The major roadblock of the Lagrangian schemes is the aliasing error, which manifests itself as noises in the tracer density distribution.During the advection, the discrete parcels (or called particles, mass packets in other works) are assumed to be isolated from each other.When the flow is nonlinearly deformative, the shape of parcels will be deformed or elongated into filaments, or even be split into parts.Most Lagrangian advection schemes do not explicitly simulate the parcel shape and assume the parcel represents the mean state of a compact volume around it.This assumption of the parcel shape is generally not valid, so large aliasing error will occur when remapping the tracer density onto the mesh of a GCM, if no interparcel mixing is considered.Different researchers designed different interparcel mixing algorithms to reduce or eliminate such errors.For example, AT-TILA redefines the parcel boundaries by simply bringing the mass mixing ratio c of a species in a parcel closer to an average background mixing ratio, or called "interaction by exchange with mean", which is similar in STOCHEM (Stevenson et al., 1998).More advanced algorithms are devised in CLaMS that uses a Lyapunov exponent to measure the deformation in the flow and inserts or merges parcels based on this exponent, and in HEL that incorporates the deformation rate of the flow and mixes tracers in a directionally biased way.Different from them, LASM explicitly describes the parcel shape by a linear deformation matrix, so a parcel can be rotated and stretched, not just translated.By doing this, the flow deformation plays an important role in the remapping process.Meanwhile, an adaptive mixing algorithm incorporating the directionally biased idea in HEL is proposed, but it is not through the mesh.In summary, any successful Lagrangian advection schemes must have an effective interparcel mixing algorithm, due to the limited resolution.
This work also discusses the problems when considering physical or chemical tendencies of tracers.A new terminator "toy"-chemistry test (Lauritzen et al., 2015) for the CTMs is conducted to verify different tendency evaluation types.To ensure the applicability of LASM in 3-D, we make several modifications and use the turbulent flow of a standard WRF large-eddy simulation to drive LASM.

LASM details
The basic formulation of LASM (e.g., linear deformation matrix, skeleton points) was introduced in Dong et al. (2014), and only key ideas are restated below.The interparcel mixing algorithm is redesigned to overcome the potential problems in the old one.The physics-dynamics coupling is also discussed, where the tendency evaluation type is the most important topic in LASM.The necessary changes are made to extend LASM into 3-D.

Key ideas
The continuous fluid is discretized into parcels in finite number N p .The centroids of parcels are {x 0i , i = 1, . .., N p }.The shape of parcel i is approximated by a linear deformation matrix H i , so in 2-D cases the parcels are represented as ellipses that can be translated, rotated and stretched.A spatial coor-dinate x can be transformed into a local coordinate system with x 0i as the origin: where y is the local coordinate for x.To avoid the problems caused by the poles on the sphere, the above calculation is currently transferred onto a local stereographic projection plane with x 0i still as the origin: where x is the projected coordinate for x.In 2-D cases, four skeleton points are associated with each parcel for sensing the deformation of the flow, and H i is determined by those points.The SVD (singular value decomposition) technique is used to decompose H i as where S i is a positive diagonal scaling matrix with the diagonal elements in descending order, and U i and V i are matrices for rotation.The product of the diagonal elements of S i is the determinant of H i .Parcel i carries an amount of mass {m s i , s = 1, . .., N s } and density {ρ s i , s = 1, . .., N s } for each tracer species, where N s is the total number of tracer species.The centroid and volume of parcel i are updated by where V is normally the gridded velocity field from the external models or data.In this study, V is on the latitudelongitude mesh or rectangular Cartesian mesh, and it is interpolated onto the parcels by using bilinear interpolation.On the other hand, the volume V i should be equal to the determinant of H i .Generally, they are different, so the determinant of H i is reset to V i by scaling S i , and H i is reconstructed from Eq. ( 5) each time step.
After the calculation of volume, the densities of all the species on parcel i are gained as where ρ s i and m s i are the density and mass of tracer species s on parcel i.The tracer density is remapped between parcels and model grids with the following formulas: where ρ s I is the density on the grid I , S i is the index set of connected grids for parcel i, S I is the index set of connected parcels for grid I , and ψ I i is the shape function value evaluated on grid I with respect to parcel i (the body coordinate is noted as y I i ).
ψ(y) is defined in Eq. ( 7) in Dong et al. (2014).Since ψ I i is shared by all the tracer species, most of the calculations need to be done once for all the species, so the multi-tracer efficiency is high in LASM, which is the same as in HEL (Kaas et al., 2013).The searching of connected grids for each parcel is based on the cover-tree algorithm (Beygelzimer et al., 2006), and it is the most costly part of LASM.In the future, when parallelizing LASM, the domain is decomposed into small ones (subdomains), and the parcels are dispatched to the corresponding subdomains, so this cost could be reduced since the searching scale for each subdomain is decreased.
At the initial time step, the density is remapped from grids to parcels by using Eq. ( 9), and the total mass on the parcels is rescaled to the one on the grids only for the first step.After rescaling, the parcel volume, the deformation matrix, and the skeleton points are all reset.Then the density on the parcels is remapped to grids each time step by using Eq. ( 10), and the total mass on the grids is rescaled to ensure the consistency between parcels and grids (i.e., a mass fixer is used).It is noteworthy that the total mass is conserved exactly on the parcels during the simulation, which is different from the traditional non-conservative semi-Lagrangian schemes.By contrast, ATTILA in ECHAM4 (Reithmeier and Sausen, 2002) does not consider the mass fixer at all.
The mass fixer will affect some properties on the grids.For example, without the mass fixer, the constant sum of three discontinuous tracers in the deformation test case (Lauritzen and Thuburn, 2012) is exactly preserved on the grids in LASM (not shown), which is hardly achieved by the Eulerian and semi-Lagrangian schemes due to the use of nonlinear shape-preserving filters.When using mass fixer, this property is lost on the grids, but the error is relatively small.It should be noted that these properties will not be affected on the parcels.
By now, the trajectories of parcels and densities of tracers are calculated, and the densities are remapped onto the grids.The next important ingredient of LASM is the interparcel mixing algorithm, which is updated as follows.

Interparcel mixing
Due to the discretization of the continuous fluid, an initial compact parcel cannot keep its integrity under the deformation of the flow.Some part of it needs to be exchanged with other parcels to form some kind of interparcel mixing, which will improve the degree of linear approximation of the parcel shape.
Previously, a parcel with index i is judged to be mixed or not according to a disorder degree D i , which is the ratio of the maximum angle to the mean angle, where the angles are between the major axes of the parcel and its neighbor parcel.When D i is greater than a given threshold, the parcel i is mixed with its neighbors by distributing tracer mass to them, and its major axis is shrunk by a fixed factor of 0.05.The volumes of the neighbor parcels are increased accordingly.D i is a good indicator of the flow condition, but this criterion and mixing rule failed in some cases, and caused extremely large and small parcels to exist in a barotropic test case (see Fig. 4a).In this study, a more effective criterion and mixing rule are designed to better control the shape of parcels.
The parcel shape is approximated by a linear deformation matrix H i as aforementioned.When the flow is linearly or quasi-linearly deformative, or the sizes of parcels are small enough compared with the deformation scale, this approximation is sufficient so that the parcels are free to sense the flow deformation.Meanwhile, the ratio between the major and minor axis length γ i (i.e., the filament degree) is limited to a maximum value γ m to keep the parcel compact, because long parcels will cause extensive computation for searching connected grids.Otherwise, the linear approximation may quickly degenerate, and this causes poor remapping, since the parcel shapes are not well simulated.The most direct symptom of this degeneration is that the positions of skeleton points (say s) deviate from the ones calculated from H i (say s ). Figure 1 demonstrates an example of degeneration, where the current shape of a parcel in red ellipse mismatches the skeleton points in green dots.The deviation is calculated by where function D(., .) is the spherical distance, and L is the half length of the major axis of the parcel.Defined in this way, d is a dimensionless number, independent of the underlying domain (e.g., Earth radius).When the maximum deviation d max of the skeleton points is greater than a given threshold d * , H i is not regarded as a good approximation to the real shape of parcel i.
Another difference from the old mixing algorithm is that the neighbor parcel shapes are not changed, so those parcels will not be disturbed by the mixing.This is similar to Eulerian schemes, where the mesh is fixed.Parcel i is not simply shrunk, and its volume is kept unchanged.The tracer density of the parcels is changed in the following rule, which The red ellipse is the current shape of the parcel, and the green points are the skeleton points of the parcel.The green points not on the red ellipse indicate that the parcel shape is not approximated well by the linear deformation matrix.
is similar to the ones used in STOCHEM (Stevenson et al., 1998) and ATTILA (Reithmeier and Sausen, 2002), but it is not through the mesh (i.e., mix parcels within a grid cells) and respects the local flow properties to a great extent.Note that the tracer species index is omitted for simplicity.Firstly, a weight w j for neighbor parcel j is calculated as before (see Eq. 18 in Dong et al., 2014), and parcel i has a unit weight.The neighbor parcel that is closer to parcel i along the major axis has a larger weight, so the weight is directionally biased as HEL.Then the mean density of parcels weighted by w j is where m and V are the mass and volume respectively, and S i is the parcel index set of parcel i and its neighbors.The density of each neighbor parcel is restored to ρ i as where C mj = w j C m , and C m is a restore coefficient that controls the mixing degree.C m should be related to the time step size and needs to be tuned accordingly.It can also be variant along the vertical spatial direction, and be proportional to the flow deformation, but it is a constant in this study for the sake of simplicity.The density of parcel i is calculated under the constraint of mass conservation as It can be seen that the parcel i is just restored to ρ i with the coefficient C m , and the mass on the parcels is exactly conserved during the mixing.
After the density mixing calculation, the shape of parcel i is changed through reducing γ i (i.e., the filament degree).In the 2-D cases, before mixing, After mixing, where S max and S min are the diagonal elements of S i .The new ones are α is a reshaping coefficient: where α m is the maximum value.This setup will make the shape of parcel i closer to circle in 2-D cases.Then the deformation matrix H i is updated, and the skeleton points are reset accordingly.This reshaping rule is the most artificial part of LASM, which may be revised in the future, but from the results of the following tests, this rule does not have vital defect.
The mixing in the Lagrangian schemes is driven by the flow deformation, whereas the inherent mixing in the Eulerian and semi-Lagrangian schemes is driven by the tracer density gradient, which is similar to the molecular diffusion.Currently, the above interparcel mixing is only for the computational aspect, but when the physical mixing is required, the similar form can be utilized with different parameters (Stevenson et al., 1998).For example, the disorder degree D i can be used as a trigger of the turbulent mixing.In addition, the stochastic turbulent diffusion in LPDM may be incorporated in the future to unify them.

Tendency evaluation
In the real applications, some other processes will calculate the tendencies for different tracer species.There are two types of tendency, which are both implemented in LASM.
The tendency in the first type resides on the Eulerian mesh of an existing model, such as the physics parameterizations (e.g., convection, microphysics) in an AGCM.The gridded tendencies need to be remapped onto the parcels, and the mass tendency is chosen to be remapped to ensure mass conservation: In other Lagrangian schemes, the tendency on a mesh cell is distributed evenly to each parcel contained in that cell (Stevenson et al., 1998;Reithmeier and Sausen, 2002).It is noteworthy that there may be inconsistency between parcels and grids, which will cause negative mass on the parcels, especially when the density is near zero, because the tendency is calculated based on the densities on the grids, and may be excessive for the parcels.When this happens, the parcels will disregard the tendencies in this study.Because the chemical reaction is local, the total mass will be conserved even so.
In the future, extra procedures should be designed to redistribute the problematic tendencies reasonably to ensure mass conservation when the tendency is non-local (e.g., convection).
The tendency in the second type is computed directly on the parcels, which is more natural and no inconsistency will occur.In the CTMs, it is easier to compute the tendencies in this type, because the chemical reactions can be described by the box model (Brunner, 2013).Several Lagrangian CTMs in this type exist, such as STOCHEM, CLaMS (McKenna et al., 2002a), etc., but NAME III calculates the chemistry on a fixed three-dimensional chemistry grid (Jones et al., 2007).Grewe et al. (2014) proposed a theoretical parcel reshaping for the finite mass method (FMM).When applying FMM in the 3-D atmosphere-chemistry models, a reshaping procedure is needed to construct the real density, because the parcel density is not real at the centroid position.In contrast, although LASM is based on FMM, the parcels in LASM carry the real density, so there is no need to reshape.There is also some research (Shima et al., 2009;Andrejczuk et al., 2010;Sölch and Kärcher, 2010;Unterstrasser and Sölch, 2014) on the microphysics, which constructed numerical schemes to simulate cloud and precipitation in the Lagrangian framework.In the future, more parameterization schemes may be designed on the fully Lagrangian framework.
The effects of the two tendency types will be compared in a newly proposed terminator "toy"-chemistry test in Sect.3.3.

Extension to 3-D
LASM can be extended to 3-D by handling the following three factors: (1) the rigid boundary condition in the vertical direction (the horizontal boundary conditions are periodic), (2) the initial shape of the parcels, and (3) the reshaping rule after mixing.
The rigid boundary condition is implemented by ensuring the vertical velocity is zero on the boundary.This is trivial in this study, since the flow is given by the external sources (e.g., WRF-LES simulation).The parcel centroids will never move outside the domain, if the time step size is not too large, but the parcel ellipsoid may penetrate through the boundary.
To avoid potential problems, all the parcels within the boundary grid cells will be reset to a sphere, and mixed accordingly.Similarly, ATTILA relocates air parcels in the boundary layer after one time step, due to the rapid turbulent mixing (Reithmeier and Sausen, 2002;Stenke et al., 2008).The interparcel mixing algorithm described above checks the parcel filament degree and ensures γ < γ m .In the 3-D applications, the horizontal and vertical space scale may differ greatly.For example, the axis span of a sigma terrainfollowing vertical coordinate is 1, but the horizontal axis span is hundreds of kilometers, so the initial parcels are like thin pies, which will confuse the algorithm.Therefore, the simulation domain is transformed to a logical one with equal axis spans in each direction.This transformation involves the flow velocity and the coordinates of grids and parcels.After that, the initial parcels are spheres.
When reshaping the mixed parcels, there is one more element S mid in the scale matrix S in 3-D.S max and S min are changed as in Eqs. ( 18) and ( 19), and S mid is changed linearly as The three elements are scaled to ensure S * max S * mid S * min equals the parcel volume.If γ * = S * max S * min > γ m , the above procedures will be looped until γ * ≤ γ m .
We utilize the 3-D large-eddy simulation (LES) in WRF to construct a turbulent flow to verify this extension.The results are analyzed in Sect.3.4.

Test cases
In this section, four tests are used to verify the updates of LASM.The first two cases are the same as used in Dong et al. (2014).The third one is newly proposed in Lauritzen et al. (2015), which extends the previous inert scalars to reactive species.The fourth one uses a more turbulent 3-D flow generated by WRF-LES.The conventional parameters of LASM are shown in Table 1.

Deformation flow test
The deformation flow test was first proposed in Nair and Lauritzen (2010), and further extended in Lauritzen et al. (2012), which contains four large-scale deformation flow formulas to challenge the advection schemes.In this study, we repeat the case with non-divergent flow to test the basic performance of the new interparcel mixing algorithm.The spatial resolution  Maximum reshaping coefficient 0.5 of the mesh is 1.5 • × 1.5 • , the parcel number is 22 640, and the normalized flow period T is 5.
The effects of changing γ m are shown in Fig. 2 and Table 2.With the decreasing of γ m , the parcels, which are under heavily stretching and shearing, cannot turn to be too long, and their tracer mass will be mixed with the neighbor parcels.Figure 2 shows the changes of mixing ratio distribution of the slotted cylinders at t = T /2 and t = T .When γ m = 2, the filament structures of the slotted cylinders break up at t = T /2.This resembles the aliasing error shown in Fig. 13 in Dong and Wang (2012), and in Fig. 7 in Kaas et al. (2013), but with a different reason.The excessive mixing caused by the small γ m distributes the tracer mass to the wrong places.Thus, γ m should not be too small to let the aliasing error happen.The correlation diagnostics shown in Table 2 are evaluated on the parcels.As expected, the real mixing diagnostics l r is also larger when γ m is small, but the "range-preserving" unmixing diagnostics l u is smaller because of the extensive interparcel mixing.All the "overshooting" unmixing diagnostics l o are zero in LASM, while other high-order Eulerian or semi-Lagrangian schemes need a filter to achieve this requirement.In spite of this, it does not mean γ m should be as large as possible.For example, in the following terminator "toy"-chemistry test, when the tendencies for the tracers are   2014) (the third row), it can be seen that the new algorithm is better to keep the final shape of the slotted cylinders, and l r is 1.97 × 10 −4 in the new algorithm, which is 3.34 × 10 −4 in the old one.
The effects of changing C m are shown in Fig. 3.By increasing C m , the mixing degree is enlarged, especially when C m = 0.1.In the current study, the value 0.001 is chosen, which is the same as STOCHEM and ATTILA.

Barotropic test
The barotropic test provides more realistic flow, though no analytical solution is available.The same finite difference barotropic model is used to drive LASM as in Dong et al. (2014).The spatial resolution is 1.5 • × 1.5 • , and the time step size is 20 s to ensure the stability of the model.The initial conditions and topography are as subcase 2 in Dong et al. (2014), but in this study, the time step size t for LASM is increased to 1800 s to decrease the cost.The mixing coefficient C m is changed to 0.01, due to the enlargement of t.  old and new criterion is shown in Fig. 4. In the old one, many extremely large and small parcels appear, but they are eliminated in the new one.Therefore, the new criterion is better, and should replace the old one.Secondly, two triggers for the interparcel mixing are compared.The first one is minimum mixing, which is the linear degeneration of the parcel shape triggers the mixing.The second one is all mixing, which is the mixing is active each time step for each parcel, but only the shapes of the degenerated parcels are changed.The results at day 4 are in Fig. 5.When C m = 0.001, the differences between minimum mixing and all mixing are imperceptible, due to the large t.After having increased C m to 0.01, the diffusion in all mixing becomes visible.Based on this comparison, LASM manifests its ability of customizable mixing, and more physics-based triggers can be devised with the information available in LASM to fulfill different needs.Another aspect is that C m should be adjusted according to different application scenarios.

Terminator "toy"-chemistry test
This test (Lauritzen et al., 2015) consists of transporting two reacting chlorine-like species (X and X 2 ) in the deformation flow.Although the chemical reactions are nonlinear, the weighted sum X T = X + 2X 2 (i.e., linear correlation) is preserved, even when the reactions are discretized.In addition, one of the reactions (X 2 k 1 → 2X) is photolytic so that the reaction rate k 1 is "terminator-like" and will produce very steep gradients in the species near the terminator.Only the schemes that are semi-linear can preserve the linear tracer correlation (i.e., X T ).Many Eulerian or semi-Lagrangian schemes satisfy the semi-linear constraint when the shape-preserving filters are not applied, but in the real cases, such schemes must be limited by the filters to avoid adverse oscillations, so the linear correlation between the two species is compromised.On the other hand, the Lagrangian schemes have the intrinsic property of being linear, so no complicated (maybe non- linear) filter or limiter is needed to ensure the monotonicity.The linear property is the sufficient condition for the semilinear property (Lauritzen and Thuburn, 2012).For example, LASM can preserve the constant sum of three tracers, which is more challenging than two tracers.
X/4.e-6 2X 2 /4.e-6X T /4.e-6 Figure 8. Cross section of X, X 2 , and X T at day 1 with the tendencies evaluated on the parcels.Results are normalized by 4 × 10 −4 (the initial value of X T ).
Except the scheme itself, the physics-dynamics (or chemistry transport in this test) coupling also affects the results, such as the tendency evaluation type.In LASM, the tendencies of X and X 2 are evaluated on the parcels and grids respectively.The spatial resolution of the latitude-longitude mesh is 1 • × 1 • .The time step size t is 1800 s both for the chemistry and advection processes, because there is no need to use a smaller t in LASM.The results for X, X 2 and X T at day 6 are shown in Fig. 6.Firstly, the constant X T is preserved very well when the tendencies are evaluated on the parcels, whereas there are errors in X T when on the grids.The results of CAM-FV and CAM-SE can be referred to in Fig. 3 in Lauritzen et al. (2015), which indicates that X T is disturbed in those schemes caused by the used shapepreserving filters.Secondly, when evaluating tendencies on the grids, small negative densities will occur near the ter-  minator line.This is caused by the inconsistency between parcels and grids.Currently, LASM ignores the problematic gridded tendencies to avoid the negative densities.In Fig. 8, the cross sections of X and X 2 at 45 • S exhibit that LASM is shape-preserving inherently without any complicated and nonlinear filter, and the sharp gradients are simulated well near the terminator.By contrast, CAM-SE shows Gibbs phenomena when no filter or limiter is used (see Fig. 5 in Lauritzen et al., 2015).In summary, it is better to evaluate the tendencies on the parcels in LASM.For CTM application, the chemical reaction can be easily calculated on the parcels, and the interparcel mixing in LASM will also improve the representation of the parcels relative to the earlier Lagrangian models (Henne et al., 2013).

3-D WRF-LES test
All the previous tests are in 2-D, so there are still concerns about the behavior of LASM in 3-D.This test will use the 3-D flow outputted from WRF-LES, which is extremely turbulent.The parameters for the LES are the standard ones in WRF, with 100 m×100 m×50 m spatial resolution in a Carte-sian domain and t = 1 s.The time duration is 2 h.The advection scheme in WRF belongs to the upwind-based finite difference method.The flux operators are fifth-order operators for the horizontal flux calculations and a third-order operator for the vertical flux calculation (Skamarock et al., 2008).The positive-definite limiter is also used to avoid negative tracer mass (Skamarock, 2006;Skamarock and Weisman, 2009).
Since WRF uses a terrain-following hydrostatic-pressure vertical coordinate η, the span of η is 1, which is too small compared with the horizontal axis span 3900 m, the coordinate system is linearly transformed into a logical one where all the axis spans are 1; so is the velocity field.An idealized tracer with initial cubic distribution is advected by both of the schemes.LASM is driven offline with a larger t = 30 s, and the parcel number is the same as the grid number of WRF.
The results in the 3-D Cartesian domain after 1.5 h for both LASM and WRF are depicted in Fig. 9.It is obvious that LASM has far less numerical diffusion than WRF; see the contour levels of the cloud-like tracer distributions on the upper levels.Meanwhile, there are a lot of structures on the lower levels in LASM, whereas the tracers are diffused heavily in WRF.Nevertheless, the diffusivity of LASM may be lower than what happens in nature, and this should be further verified in the real application where the observation is available.The animation for both of WRF and LASM can be found in the supplement files.
One major concern about LASM is the cluster and rarefaction of parcels.For example, is it possible that a parcel has no or only distant neighbors?It should be noted that the distribution of the parcels is controlled by the flow, which is constrained by the fluid dynamics in turn.In Reithmeier and Sausen (2002), their results showed that the trajectories of the parcels during the 30-year control run remain well distributed (see their Fig.10).A similar figure for LASM is shown in Fig. 10, which shows that the horizontally averaged number of parcels contained in the cells has no significant trend.LASM will also crash when a parcel has no neighbors in this study, but this never happens in the tests.In spite of this, the adaptive refinement of parcels is still valuable for achieving better accuracy and performance, and this will be studied in the future.

Conclusions
The interparcel mixing algorithm in LASM is updated by replacing the criterion of the linear degeneration of the parcel shape and the mixing rule.The new criterion is a more direct indicator of the degeneration symptom, so the parcel shape is better controlled.In the new mixing rule, the densities of the involved parcels are restored to a mean density.
Several tests are utilized to verify the updates, including a newly proposed terminator "toy"-chemistry test.The results reveal that the new algorithm is effective and less arbitrary.It is noteworthy that LASM is a linear scheme, so it can preserve the sum of multiple tracer species and also the weighted sum of two reactive species exactly on the parcels and accurately on the grids.Two triggers for the interparcel mixing are compared, and LASM shows great flexibility so that different triggers can be devised in the future to address different needs.
The two tendency evaluation types are discussed and tested.The first one, which evaluates the tendency on the parcels, is the most natural type.Some processes are trivial to be calculated in this type, such as the chemical reaction and microphysics.The other one, which evaluates the tendency on the grids, can cause negative values on the parcels, which should be carefully handled.
LASM is also extended to 3-D by handling the rigid boundary condition, adjusting the initial parcel shape, and changing the reshaping rule after mixing.The flow outputted by the WRF-LES is used to verify this extension of LASM, and the results of LASM and WRF indicate that LASM has far less numerical diffusion.The parcel distribution also shows that there is no cluster and rarefaction of parcels.

Figure 1 .
Figure1.A case of linear approximation degeneration of a parcel.The red ellipse is the current shape of the parcel, and the green points are the skeleton points of the parcel.The green points not on the red ellipse indicate that the parcel shape is not approximated well by the linear deformation matrix.

Figure 2 .
Figure 2. Deformation test case results with different γ .The left column is the slotted cylinders at t = T /2, and the right column is at t = T .

Figure 3 .
Figure 3. Deformation test case results with different mixing coefficient C m at t = T .
96 × 10 −4 8.78 × 10 −9 0 considered, the larger γ m causes the burr edge along the terminator as shown in Fig.7.Therefore, the moderate value 5 for γ m is used in LASM.By comparing Fig.2(the third row with γ m = 10) with Fig.12inDong et al. (

Figure 4 .
Figure 4. Comparison of the parcel shapes between the old and new interparcel mixing algorithms in a barotropic test case.Two extremely small and large parcels in the old algorithm are spotted by the red ellipses.In the new algorithm, there is no such ill-shaped parcels.

Figure 5 .
Figure 5.Comparison between minimum mixing and all mixing at day 4.All mixing means the interparcel mixing is triggered each time step for each parcel.The mixing coefficient C m has two different values: 0.001 and 0.01.

Figure 6 .Figure 7 .
Figure6.Contour plots of X, X 2 and X T at day 6 for two types of tendency.Solid black line is the location of the terminator line.

Figure 9 .
Figure 9. 3-D volume-rendered advection results driven by the flow outputted by WRF-LES.The left figure is simulated by LASM, and the right one is by the advection scheme in WRF.

Figure 10 .
Figure 10.Evolution of the parcel distributions during the WRF-LES simulation.Shown is the horizontally averaged number of parcels contained in the grid cells.Each line represents a level.

Table 1 .
The conventional parameters of LASM used in the tests.

Table 2 .
The correlation results evaluated on the parcels for different values of γ m in the deformation test case 4.