Lagrangian description of the atmospheric flow from Eulerian tracer advection with relaxation

The Lagrangian representation of fluid flows offers a natural perspective to study many kinds of physical mechanisms. By contrast, the Eulerian representation is more convenient from a diagnostic point of view. This article attempts to combine elements of both worlds by proposing an Eulerian method that allows one to extract Lagrangian information about the atmospheric flow. The method is based on the offline advection of passive tracer fields and includes a relaxation term. The latter device allows one to run the integration in a continuous fashion without the need for reinitialization. As a result one obtains accumulated Lagrangian information, for example, about the recent parcel displacement or the recent parcel‐based diabatic heating, at each point of an Eulerian grid at any time step. The method is implemented with a pseudospectral algorithm suitable for gridded global atmospheric data and compared with the more traditional trajectory method. The method's utility is demonstrated on the basis of a few examples, which relate to cloud formation and the development of temperature anomalies. The examples highlight that the method provides a convenient diagnostic of parcel‐based changes, paving an intuitive way to explore the physical processes involved. Due to its gridpoint‐based nature, the proposed method can be applied to large data sets in a straightforward and computationally efficient manner, suggesting that the method is particularly useful for climatological analyses.

they change with time. However, numerical atmospheric models usually take the Eulerian perspective, which means that they describe the flow in terms of fields on a fixed grid. The Eulerian perspective is very convenient from a diagnostic point of view, since variables are available at known grid points and thus can easily be processed further; moreover, it also conforms to many observing systems. The disadvantage of Eulerian models is that their output usually does not allow direct access to Lagrangian information.
A widely used method for teasing out Lagrangian information from gridpoint-based fields is the calculation of trajectories (Reed, 1955;Danielsen, 1961). Trajectories allow one to study a variety of meteorological phenomena such as heat waves (e.g. Harpaz et al., 2014;Zschenderlein et al., 2019;Spensberger et al., 2020;Zschenderlein et al., 2020), extratropical cyclones and warm conveyor belts (e.g. Wernli, 1997;Eckhardt et al., 2004;Madonna et al., 2014), troposphere-stratosphere exchange (e.g. Wernli and Bourqui, 2002;Birner and Bönisch, 2011;Škerlak et al., 2014), moisture transport (e.g. Stohl et al., 2008;Knippertz and Wernli, 2010;Ryoo et al., 2015), and air pollution (for references see Fleming et al., 2012, their Table  1). A wealth of information can be obtained from trajectories, as they provide time-resolved knowledge about the physical properties along the pathways of individual parcels. For many applications this is useful or even necessary, but for certain applications this may be more than is needed. For instance, it is sometimes sufficient to restrict attention to time-accumulated Lagrangian information, and in such a situation it would be convenient to have a method that is tailored to provide accumulated Lagrangian information only.
In the current article, we propose such a method, which combines elements of both Lagrangian and Eulerian worlds. The method is based on the advection of passive tracers in combination with a relaxation term. A passive tracer can be any kind of substance or physical quantity that is advected by the flow without feeding back to the dynamics. This allows one to track air motion and to visualize the result of atmospheric transport through the use of passive tracers. For instance, chemical species with distinct sources in the atmosphere such as carbon monoxide or ozone have been used as tracers, allowing one to measure (or model) their atmospheric concentration and hence retrieve important information about atmospheric transport pathways. While this method has often been applied in the context of the stratospheric circulation and troposphere-stratosphere exchange (e.g. Brewer, 1949;Dobson, 1956;Fischer et al., 2000;Hoor et al., 2002;James, 2003;Bönisch et al., 2009;Hegglin et al., 2009;Hoor et al., 2010;Ploeger and Birner, 2016;Ploeger et al., 2017), passive tracers have been used much less frequently to study tropospheric dynamics. This is what we aim to focus on in the present article.
Our method builds upon an approach initially proposed by Schär and Wernli (1993) in the framework of a semigeostrophic model and explored further by Gheusi and Stein (2002) with regard to mesoscale numeric modelling. Their approach is based on the advection of three Eulerian passive tracer fields initialized with the coordinates of each grid cell. Likewise, we draw on the ideas inherent in the potential vorticity tracer technique by Davis et al. (1993), which was later applied, for example, by Stoelinga (1996), Gray (2006), or Chagnon et al. (2013). Their technique is based on the advection of potential-vorticity-like tracer fields, each of which is subject to a modification by a certain diabatic process. Similarly to the aforementioned methods, our method uses fields of latitude, pressure, or potential temperature as passive tracers on an Eulerian grid. The novelty of our method is that it incorporates a relaxation term, which turns out to be beneficial for practical reasons.
It is the goal of this article to present the underlying idea of our method and motivate its utility for various types of Lagrangian analysis on the synoptic scale. The article is organized as follows. In Section 2, we outline the basic idea of the method. Section 3 presents an algorithm for the implementation to gridded model data on the synoptic scale. In Section 4, we compare the results of our algorithm with the results from trajectory calculations, before we provide a few short, illustrative examples in Section 5 regarding the method's utility. Finally, a summary is given in Section 6.

Basic idea
We start by presenting the basic idea with the help of an example. Consider the vertical displacement Δ Z of an air parcel, which is defined as where Z(t) denotes the parcel's altitude at time t and t 0 denotes the initial time. The definition implies that the vertical displacement is zero at the initial time, that is, Taking the time derivative of Equation 1, one obtains where represents the vertical motion of the parcel, that is, the vertical wind at the location of the parcel. Note also that we use the notation d∕dt for the time derivative of a variable that depends on time only. If we consider Equation 3 as a differential equation for Δ Z (t), this can be integrated to yield where Equation 2 was used to determine the constant of integration. Note that Equation 5 is not an integral of the vertical wind at a certain location; instead, as denoted by Equation 4, it is an integral of the vertical wind at the position of the parcel, that is, an integral along the trajectory of the parcel. Equation 5 yields the vertical displacement of the air parcel within the time interval of length t − t 0 . The length of this time interval is growing with time t, since t 0 is the (fixed) initial time. However, for some applications it would be desirable to limit this time interval to a finite length like a window that moves with the considered time t. This motives us to add a relaxation term on the right-hand side of Equation 3, yielding the following modified initial-value problem for the evolution of Δ Z (t), namely The relaxation term in itself effectively models an exponential decay of Δ Z (t) with the rate > 0. As we will see shortly, this decay results in a finite memory of the parcel for its vertical displacement. From a mathematical point of view, Equation 6 is an inhomogeneous ordinary differential equation with W(t) representing the inhomogeneous term, that is, the forcing. The equation can be solved using mathematical standard methods such as separation of variables and variation of the constant. The solution of Equation 6 reads which can be readily verified by direct substitution into Equation 6. Obviously, for = 0 this solution recovers Equation 5, which means that the solution takes into account the vertical wind W(t ′ ) throughout the integration with equal weight for any time t ′ . By contrast, for > 0, the contribution of W(t ′ ) at earlier times is exponentially damped away, such that the solution makes substantial use of information from W(t ′ ) only during the recent past t ′ > t − −1 . In other words, a positive value of represents a finite memory of the parcel for accumulating the vertical wind in the expression for the vertical displacement. The equation for Δ Z (t) can be reformulated in an alternative manner, which turns out to be useful in practice (see below). With the help of Equation 4, Equation 6 becomes Defining which satisfies the initial condition Ψ(t 0 ) = Z(t 0 ), and taking the time derivative, we obtain the following initial-value problem: The solution Ψ(t) can be used to retrieve the vertical displacement through Equation 9, that is,

Generalization
Now consider a more general pair of parcel variables (A,Ȧ) satisfyingȦ In the special case of the previous example, these generalized variables can be identified as (A,Ȧ) = (Z, W) owing to Equation 4. A generalized displacement Δ A (t) is then, by analogy with Equation 6, defined through and Equations 10 and 11 turn into from which the desired quantity Δ A (t) can be obtained through Pairs of variables satisfying Equation 12 can be, for example, where X denotes the position of the parcel in the zonal direction, Y is the position of the parcel in the meridional direction, P is the pressure of the parcel, Θ its potential temperature, U is the zonal wind, V is the meridional wind, Ω is the vertical wind in pressure coordinates, and Q = dΘ∕dt is the diabatic heating, defined as the material rate of change of potential temperature, all of these taken at the position of the parcel at time t. The corresponding "generalized displacements" are denoted as Δ X , Δ Y , Δ P , and Δ Θ , and they represent the parcel's displacement in the zonal direction, the parcel's displacement in the meridional direction, the parcel's vertical displacement in terms of pressure, and the accumulated diabatic heating (as defined through potential temperature), respectively. The special case of the (Θ, Q) pair explains why the formulation in Equation 14 may sometimes be more straightforward to apply than the formulation in Equation 13: the former requires only knowledge about the parcel's potential temperature (which usually is readily available), while the latter requires explicit knowledge about the diabatic heating (which is typically much harder to obtain from, or even unavailable in, standard data sets).

Formulation as an Eulerian method
So far we have only considered the evolution of the property of a specific parcel, which corresponds to the Lagrangian perspective. This resulted in a simple inhomogeneous linear ordinary differential equation. We now shift to the Eulerian perspective and consider the parcel displacement as a function of space and time, that is, a (x, t), where x = (x, y, z) denotes the parcel's location in three-dimensional space and t is time (as before). This change in perspective can be achieved through replacing the time derivative along the air parcel, d∕dt, by the total (or material) derivative D∕Dt operating on the respective Eulerian field. The material derivative can be related to partial derivatives in time and space through the following well-known relation: Here, u(x, t) = (u, v, w) denotes the three-dimensional Eulerian wind field, ∕ t is the partial time derivative, and ∇ = ( ∕ x, , y, ∕ z) is the Nabla operator. Effectively, Equation 17 turns Equation 13 into an equation for the evolution of a (x, t), namely Note that we switched consistently from upper-case to lower-case notation in order to clarify that we are now dealing with Eulerian fields, for example, (x, t) instead of Ψ(t) and a (x, t) instead of Δ A (t). In analogy with Equation 14, the problem in Equation 18 is equivalent to solving a problem for (x, t): in combination with Both formulations, that is, Equation 18 and Equations 19 and 20, represent an advection-relaxation problem, which can be solved on an Eulerian grid using standard numerical methods. As pointed out earlier, the formulation in Equations 19 and 20 is potentially easier to use in practice because the inhomogeneous term on the right-hand side is a instead ofȧ. The formulation in Equations 19 and 20 constitutes the basis of what we call our tracer method. This method involves the advection of a passive tracer field (x, t), the values of which are gradually relaxed towards the field a(x, t). The major advantage of the relaxation is that it fades out the past "on the fly" and thus allows one to carry on with the integration uninterrupted over an unspecified and potentially very long period of time; in other words, there is no need for any reinitialization. This avoids any recalculation of overlapping time intervals, making the tracer method potentially very efficient in certain cases. Unfortunately, there is no option to choose the weighting kernel representing the gradual fade-out of the earlier parcel information freely; rather, the inclusion of the relaxation term in the differential equation simply results in an exponential kernel in the solution. However, the decay time associated with this exponential kernel can be controlled by the user via the parameter .
Note that the interpretation of the resulting fields in terms of accumulated Lagrangian information during a specific time interval requires that the integration time is considerably longer than −1 . As a consequence, we consider the time t ≤ t spinup = −1 as a spin-up time. For integration times shorter than t spinup , the effective window size for accumulation from the past increases with t, while for integration times longer than t spinup the effective window size becomes independent of t. The latter property is desirable for a consistent interpretation of the results.

TECHNICAL IMPLEMENTATION
Obviously, there are various options for the numerical implementation of the general method outlined above. The optimum choice is likely to depend on the type of input data, whether they are regional or global, two-dimensional or three-dimensional, and so forth. In the following, we present the algorithm chosen for this work along with the input data we use.

Pseudospectral algorithm
Our algorithm is designed to solve the advection-relaxation problem of Equations 19 and 20 for given gridded data on a three-dimensional, global, atmospheric domain. We apply the pseudospectral method (Orszag, 1972), which uses gridpoint-based as well as spectral field representations. The algorithm's key characteristics are the following (see the pseudocode in Figure S1 for illustration).
• Input fields are given on a global regular latitude-longitude grid with a height coordinate that is terrain-following in the lower troposphere.
• Horizontal spatial derivatives are computed in spectral space using a series expansion in spherical harmonics with triangular truncation. This representation implies a uniform spatial resolution over the entire globe, which is beneficial in terms of numerical stability (Durran, 2010).
• The vertical gradient is approximated by second-order centred differences.
• Product terms of fields are computed in physical space.
• Time integration is performed in spectral space using the explicit third-order Runge-Kutta scheme proposed by Williamson (1980).
• For numerical stability, a hyperdiffusion term of the form D = − ⋅ ∇ 4 (with being the diffusion coefficient) is added to the right-hand side of Equation 19. The term compensates for aliasing errors that are inherent in pseudosepctral methods and can cause nonlinear instabilities (Durran, 2010). In spectral space, the term ∇ 4 can be evaluated in a straightforward manner, as the spherical harmonics are eigenfunctions of the Laplacian operator on the sphere. Generally, we try to keep the amount of hyperdiffusion as low as possible. As a consequence, we choose the magnitude of dependent on the grid spacing, since a finer horizontal grid spacing turns out to require less hyperdiffusion for numerical stability. We determined the minimum necessary amount of hyperdiffusion by trial and error and found = 3 × 10 13 m 4 ⋅ s −1 and = 3 × 10 14 m 4 ⋅ s −1 for horizontal resolutions of 0.5 and 1 • , respectively, to be acceptable for our application.
To be sure, the algorithm could be fine-tuned and improved with respect to accuracy, stability, or efficiency. In particular, our algorithm is limited by the requirement of a relatively small time step to ensure numerical stability. The use of more advanced methods such as implicit or semi-implicit schemes would allow a much larger time step, which presumably leads to enhanced efficiency. However, we did not attempt any further optimization in this regard, since the focus of this article is on proof of concept and illustration of the method.

Application to reanalysis data
We apply our algorithm to data from ERA5 (Hersbach et al., 2020). ERA5 is the latest reanalysis product of the European Centre for Medium-Range Weather Forecasts, providing hourly estimates of a large number of atmospheric, land, and oceanic variables. Unless noted otherwise, we retrieved the data (Hersbach et al., 2017;Hersbach et al., 2018b) with a horizontal resolution of 1 • and a temporal resolution of 3 h. Since required by our algorithm, we interpolated the data linearly to smaller time steps. In the vertical, we use the ERA5 hybrid model levels (Simmons and Burridge, 1981;ECMWF, 2016). In our applications, we restricted ourselves to the lowest 88 model levels, which encompass the entire troposphere as well as the lower stratosphere. To generate the final output, we performed a linear interpolation from model levels to pressure levels.

COMPARISON WITH TRAJECTORIES
The traditional technique to retrieve Lagrangian information from gridpoint-based data involves the computation of trajectories, and there is a large body of literature where this technique has been applied successfully (see references in Section 1). Since our new method also aspires to provide Lagrangian information, it seems desirable to compare these two methods. We do this by first highlighting the conceptual relationship between the two methods, and second comparing the results of our method with results from trajectory calculations for the special case of no relaxation (i.e., = 0). This comparison allows us to test the core of our developed pseudospectral algorithm, namely the advection of tracer fields, and to point out limitations of the methods generally.

Conceptual relationship between the tracer method and the trajectory method
Any trajectory is meant to represent the location X(t) of a parcel in three-dimensional space at time t. As pointed out later, this concept is well-defined as long as the parcel is infinitesimally small. The trajectory is obtained through where U(t) denotes the three-dimensional wind at the parcel's location at time t. In typical applications, U(t) is not known a priori; rather the flow is given in Eulerian representation u(x, t), that is, on a regular spatial grid at discrete points in time. In that case, the integrand U(t ′ ) must be obtained by interpolation in space and time.
The knowledge of a trajectory allows one to retrieve detailed information about the respective parcel. In particular, the zonal, meridional, and vertical displacements are given by respectively, thus providing information about Δ X (t), Δ Y (t), and Δ Z (t) at any time in the history of the parcel. Other Lagrangian displacements can be obtained by integrating the respective variable in time as the parcel travels along its trajectory. Per design, this information is restricted to the points on the trajectory, which is a one-dimensional manifold. To obtain information about an atmospheric volume at a given time, one has to consider a large number of trajectories that terminate in the volume of interest at that time. In addition, since usually the flow field is nonstationary, this computation has to be repeated many times. In theory this is straightforward, and even in practice it has been shown to work well, given modern computing infrastructures and the possibility of parallelisation (e.g. Madonna et al., 2014;Škerlak et al., 2014;Kremer et al., 2020). However, in certain situations we suggest that our method based on Eulerian tracer fields provides a more convenient approach. The main advantage of our tracer method lies in the fact that it provides accumulated Lagrangian information in the form of gridded fields available at any time step. Potentially, this allows one to analyse large data sets such as reanalysis data or output from climate models using standard Eulerian techniques without the need for any data preparation beforehand. The trade-off is that our Eulerian method does not yield time-resolved information along parcel trajectories; rather, it only provides accumulated Lagrangian information about the recent past, where "recent past" is determined by the time-scale −1 . In addition, each Lagrangian property requires its own tracer variable. However, as long as the number of properties is not too large, the increase in computational cost is modest. We thus conclude that our tracer method could facilitate the Lagrangian analysis substantially in those situations where the focus is on accumulated information about a modest number of parcel properties. For specific research questions, our tracer method should, therefore, be a useful tool complementing the trajectory method-even though in principle trajectories allow one to retrieve the same information.

Calculation of pseudotrajectories
As explained above, the tracer method and the trajectory method do not provide exactly the same information.
Thus, an attempt to compare the two methods requires one to either map the results from the tracer method to those of the trajectory method or vice versa. Here, we choose the first way and explain below how we extract information from the tracer method in such a way that it can be compared directly with the trajectory method. Note that this exercise requires quite some effort and is carried out exclusively for the purpose of comparison and validation. The explicit calculation of particle trajectories is and remains the core competence of the trajectory method.
The computation of a backwards trajectory provides time-resolved information about the location of the parcel that eventually reaches a specific grid point at time t. For instance, storing the parcel's location every 3 h during the integration, one obtains the relevant Lagrangian information at times t, t − 3 h, t − 6 h, t − 9 h, and so forth. Figure 1a serves for illustration. We can extract the same information from our tracer method by applying the method in a somewhat unorthodox way. First, as mentioned before, we set = 0, because the computation of traditional trajectories does not involve any relaxation and the associated memory about parcel properties is effectively assumed to be infinite. Using a separate tracer for each of the three geometric dimensions, we are able to describe the net parcel displacement in each spatial direction at any given grid point: x (x, t), y (x, t), z (x, t). However, by design these Eulerian fields only provide aggregated parcel information: they represent the net displacement in each of the variables over the duration of the entire integration (see Equation 5). Of course, this information is less detailed than what can be obtained from the true trajectory with output every 3 h.
Without careful consideration, it may seem as if the issue could be addressed by splitting the Eulerian integration into segments of 3 h duration, since this procedure would yield the displacements during each of the 3-h time intervals, which could then be concatenated. Unfortunately, however, this concatenation would not be straightforward, since the tracer method (only) provides information at grid points, whereas the parcel travels along the true trajectory and is almost never located right at a grid point except (by construction) at the final time.
We circumvent this problem in the following way. Consider a specific grid point for which a backward trajectory was computed (Figure 1a). We now perform several independent integrations of increasing length (3, 6, and 9 h, , and so forth. These integrations are represented as blue arrows in Figure 1b; they provide information about the net displacement of the parcel along the true trajectory during the corresponding time intervals , and so forth. The displacements associated with the 3-h segments along the trajectory can then be obtained from the differences between the net displacements of adjacent intervals (red arrows in Figure 1b). We refer to the procedure described as the calculation of "pseudotrajectories". Note that we calculate backwards pseudotrajectories by integrating forward in time. To the extent that the tracer advection is not fraught by substantial numerical diffusion, the pseudotrajectories should yield very similar results to traditional trajectories.

Comparison of pseudotrajectories with traditionally computed trajectories
We are now going to compare traditional trajectories with pseudotrajectories using data from a specific episode. The chosen episode contains cyclone Vladiana, which crossed the Atlantic in September 2016 (Schäfler et al., 2018) and exhibited a strong warm conveyor belt (Oertel et al., 2019). The latter was associated with large vertical displacements, which makes it a suitable case for our purpose.
First, we calculated an ensemble of traditional trajectories with the well-established Lagrangian analysis tool LAGRANTO (Wernli and Davies, 1997;Sprenger and Wernli, 2015). The ensemble encompasses backwards trajectories initialized on September 25, 2016 at 1200 UTC over central Europe (20 • W -50 • E and 40 • N -90 • N) on multiple vertical levels in the troposphere. We computed the trajectories from three-hourly ERA5 data on a grid with 0.5 • horizontal resolution. We then calculated pseudotrajectories with our tracer method; for that purpose we used a grid of either 0.5 or 1 • horizontal resolution, allowing us to test the sensitivity of our results with respect to grid spacing. The results from the comparison are presented in the following two subsections, in which we focus on relatively short and relatively long integration times, respectively.

Comparison for short integration times
We start the comparison using a situation with a relatively short integration time of two days. As detailed Δ z = -300 m Trajectory which at time t is at a specific grid point: in the next section, on such a short time-scale the concept of an air parcel can generally be considered as sound for synoptic-scale flow; moreover, the error of a trajectory is likely to be small (e.g. Stohl et al., 2001), such that the trajectory method can be regarded as "ground truth".
The right column in Figure 2 depicts two-day LAGRANTO backwards trajectories terminating at a specific location over the Atlantic. The same figure also shows the corresponding pseudotrajectories calculated with the tracer method, using the two different grid spacings. The figure indicates that, for all three variables considered, the This ascent is accompanied consistently by strong cooling (Figure 2d-f), attributable to adiabatic expansion, and diabatic warming (Figure 2g-i), attributable to condensation of water vapour (e.g. Eckhardt et al., 2004;Madonna et al., 2014). Furthermore, both the pseudotrajectories and the LAGRANTO trajectories are very consistent with regard to their horizontal patterns ( Figure 3). For yet another perspective, in Figure 4a-c we show horizontal maps of the net displacement in terms of pressure at 500 hPa after two days of integration. Again, all three panels exhibit similar behaviour, showing consistent displacements in terms of pressure. The overall good agreement between the pseudotrajectories and the LAGRANTO trajectories for short integration times provides confidence that our implementation of the tracer method works adequately from a technical perspective.
Despite the overall good agreement, there are also systematic differences: the results obtained from the pseudotrajectories are generally smoother than the results obtained from the LAGRANTO trajectories, especially when using the coarser grid for the pseudotrajectories. The underlying reason for this difference is the fact that the pseudotrajectories are based on an Eulerian technique, that is, tracer advection on an Eulerian grid, and Eulerian techniques are generally known to be subject to diffusion. Part of the diffusion is numerical diffusion, and another part in our algorithm is due to the hyperdiffusive term that we incorporated for numerical stability. As described earlier, we adjusted the coefficient to the horizontal resolution with the aim of minimizing the amount of hyperdiffusion. We therefore tested the effect of the magnitude of on a fixed grid. Figure 5 reveals that, indeed, the observed ascent within the warm conveyor belt depends on the magnitude of , with stronger ascent for weaker hyperdiffusion. The hyperdiffusion is a drawback of our algorithm, but unavoidable in its current implementation. A less diffusive advection algorithm would be desirable, especially for the coarser grid resolution; however, we will never be able to eliminate effects from numerical diffusion completely, because the tracer method is an Eulerian technique by design.

Comparison for longer integration times including the discussion of limitations
We now proceed with the method intercomparison for longer integration times. For an integration time of 10 days, the displacement fields obtained from LAGRANTO trajectories look spotty and rather detailed (Figure 4f), while the inherent diffusion of our tracer method yields considerably smoother fields (Figure 4d,e). These differences are more substantial than those that we found for short integration times. We do not consider this fact as an indication of failure of either method, but the results from both methods are now likely to be subject to uncertainties and depend on the underlying assumptions. In particular, it is not so clear any longer whether the trajectories can be used as "ground truth", as we did for short integration times.
Assessing and validating Lagrangian information is difficult. In contrast to Eulerian fields produced by a numerical model, Lagrangian diagnostics can hardly be observed or measured directly. Sometimes in situ or satellite observations of chemical tracer concentrations such as ozone or carbon monoxide could be used as a reference; however, it is sometimes difficult to extract the desired Lagrangian information from chemical tracers such as diabatic heating or meridional displacement. Therefore, we refrain from such a comparison here; instead, we aim to sketch the limitations of both methods in order to facilitate a safe interpretation of the results.
Limitations generally arise from computational errors as well as conceptual assumptions, both of which may become critical for long integration times. Computational errors (that is, inaccuracies in the numerical solution) result from errors in the underlying wind field, limited spatial and temporal resolution of the wind field, and errors due to finite-difference approximations. The contribution from the finite-difference approximations is likely to be larger for the tracer method than for the trajectory method, owing to the numerical diffusion inherent in gridpoint-based techniques. Stohl (1998) estimate that the computational error of a trajectory amounts on average to 20% of the distance a parcel has travelled after a few days. In this context, errors resulting from the numerics are generally considered to be much smaller than those resulting from the representation of the wind field (Stohl et al., 2001).
Despite rather small computational errors, the trajectory method sometimes fails to explain observed phenomena without any apparent reason (Stohl et al., 2002). This problem is likely to arise from the implicit assumption that air parcels are infinitesimally small. In reality, any air parcel takes up a finite volume; the latter implies that the value of a variable at the location of the parcel is not properly defined any longer, because it may vary over the volume of the parcel. The issue is highly relevant in flows with deformation, where parcels may be stretched out into thin filaments of exponentially increasing length. This effect is often referred to as "stirring" (e.g. Haynes, 1999;Haynes, 2005) and implies that after some time each parcel samples widely different regions of the atmosphere and cannot be considered as representative of a specific point in space any longer. By that time, the concept of a parcel has lost its proper meaning. The issue is particularly severe in the case of flow bifurcations, as illustrated beautifully by Welander (1955; their figures 2 and 3). On top of that, filamentation resulting from stirring may ultimately lead to "true mixing", which is the homogenization of tracer fields by molecular diffusion (e.g. Haynes, 1999;Haynes, 2005). Thus, the notion of a trajectory becomes-sooner or later-meaningless.
Quantifying trajectory uncertainty owing to the process of filamentation is difficult. Stohl et al. (2002) try to provide an estimate based on trajectory simulations with a particle dispersion model. They find the average uncertainty to be of the order of 20%-40% of the distance a parcel travelled after 10 days and state that "this is more than the errors resulting from wind field analysis, interpolation, and numerical truncation errors." Of course, these values should only be seen as a rough estimate; they may vary strongly with the size of the air parcels as well as the specific flow situation. However, filamentation as a result of stirring seems to be of vital importance when interpreting the significance of trajectories.

The role of the relaxation term
We close the section by highlighting the role of the relaxation term of our tracer method. This term does not usually appear in traditional trajectories methods. As explained in Section 2.1, the relaxation term effectively introduces an exponential weighting kernel, such that the "recent past" is accounted for to a larger extent than the "distant past" in the parcel's memory for its displacement. Translated into the framework of trajectory uncertainty, this means that information from the "recent past", which is generally less prone to error, is accounted for with a larger weight than information from the more "distant past", which is generally more prone to error. This design feature can thus be interpreted as a way of accounting for uncertainty in the parcel's memory for its displacement. Owing to the exponential kernel function, the transition between "recent" and "distant" is gradual, which is desirable, as the growth of errors along trajectories is likely to be gradual rather than abrupt. Of course, the interpretation of the relaxation as a suppression of the accumulating errors has to be taken with some caution, because error growth may be flow-dependent, which is not accounted for with a fixed relaxation constant. A future implementation of our method may account for a flow-dependent value of , for example, as a function of the shear or the deformation of the wind field. However this is beyond the scope of the present article.

ILLUSTRATION OF THE METHOD
In the following, we are going to illustrate the intended use of our tracer method, which is to employ the method with a relaxation constant > 0. We exemplify this on the basis of a few simple examples. The first two examples link displacements in terms of pressure and parcel-based diabatic temperature changes to patterns in cloud fractions (similar to Schär and Wernli, 1993), while the third and forth examples relate meridional displacements to patterns in temperature anomalies and the Lagrangian mean circulation. In all examples we apply a relaxation constant of −1 = 3 days. Finally, we discuss the sensitivity of the results to the choice of this parameter. Our examples are qualitative in nature and meant to be intuitive in order to inspire the reader to further analysis.

Example 1: The vertical displacement and its relation to clouds
A prominent feature of the Earth's atmosphere is the occurrence of clouds. Clouds form in many parts of the atmosphere: for instance cyclones, convective systems, or the vicinity of orography. Figure 6a-c shows an example of the evolution of cloud fractions in the upper troposphere. The question we want to address here is whether we can understand, to some extent, why the patterns in the field of cloud fraction look as they do.
Although the details of cloud formation are complex, involving various processes all the way to the microphysical scale, a major forcing mechanism is simply ascending motion. Ascent is associated with decreasing pressure and adiabatic cooling, which increases relative humidity and eventually leads to saturation. Once formed, clouds often behave like a tracer, that is, they get advected by (j,l) shows the displacement in terms of pressure; blue denotes ascended air masses, whereas red denotes descended air masses. (m-o) depicts the displacement in terms of potential temperature; blue denotes diabatically cooled air masses, whereas red denotes diabatically heated air masses the wind, as long as they do not rain out or dissolve through adiabatic compression in descending motion. Thus, to broadly understand the patterns of cloud fractions in Figure 6a-c, one needs to diagnose which of the air masses have undergone ascent versus descent in the recent past. By contrast, instantaneous vertical wind fields (Figure 6d-f) do not provide this kind of insight, since they only contain information about the instantaneous vertical motion, but not the (recent) past. Similarly, wind fields that have been time-averaged on the Eulerian grid (Figure 6g-i) are inappropriate because they lack parcel-based information. Both deficiencies can be overcome by using a Lagrangian method such as the tracer method.
Applied to the variable pressure, the tracer method yields the accumulated change of pressure along a parcel's trajectory. Usually, this Lagrangian change in terms of pressure can be taken as a proxy for the Lagrangian change in altitude. Correspondingly, the Lagrangian change or "displacement" in terms of pressure indicates which of the air masses have undergone ascent or descent in the recent past (Figure 6j-l). Besides small differences, the displacement in terms of pressure indeed resembles the cloud fractions and their evolution in Figure 6a-c. Cloudy regions mostly coincide with air masses that have undergone major ascent, whereas cloud-free regions mostly coincide with air masses that have undergone only slight ascent or even descent. With our tracer diagnostic we can, thus, demonstrate that ascending motion is indeed a major forcing mechanism for the formation of clouds.

Example 2: The parcel-based diabatic temperature change and its relation to latent heat release
Further insight into Example 1 can be obtained by considering an additional tracer, namely potential temperature. Similarly to the displacement in terms of pressure, the displacement in terms of potential temperature (Figure 6m-o) marks those air masses that have undergone diabatic heating or cooling, respectively, in the recent past. The plots reveal strong signs of diabatically heated air masses mainly in cloudy regions, consistent with the latent heat release associated with the cloud formation process during the recent past. By contrast, diabatically cooled air masses are mostly correlated with cloud-free areas, consistent with the general tendency towards radiative cooling of the free atmosphere.

Example 3: The meridional displacement and its relation to temperature anomalies
In the following, we want to illustrate that temperature anomalies in the extratropics are sometimes strongly associated with horizontal temperature advection. We address this topic from a Lagrangian point of view, linking horizontal temperature advection to meridional parcel displacements.
As an example, Figure 7a-c shows temperature anomalies in the lower troposphere for a six-day period in September 2016, which was associated with quite large warm anomalies over Europe (Zschenderlein et al., 2018). In particular, on September 6 (Figure 7a), a tongue of unusually warm air lay over Europe, extending from France via England to Scandinavia. On September 12 (Figure 7c), almost all of Europe experienced significantly higher temperatures than usual. In between (Figure 7b), somewhat lower temperature conditions prevailed in large parts of Europe.
To find out whether and to what extent the temperature anomalies can be related to meridional parcel displacement, we apply the tracer method to a latitude tracer. The resulting fields then indicate which air masses have originated from the south and which from the north during the recent past. Given that the climatological mean temperature gradient is mostly aligned in the meridional direction, one may generally expect an anomalously strong poleward displacement to be associated with a local warm anomaly and an anomalously strong equatorward displacement with a local cold anomaly.
In order to compare the displacement fields and the temperature anomalies in a meaningful manner, we provide both as standardised anomaly fields computed as deviations from a two-month average (see Section 5.4). By analogy with the temperature anomaly field, the meridional displacement field then shows the extent to which it is unusual compared with the climatology. In other words, the field indicates whether air masses originate more from the north or more from the south than usual. If temperature anomalies are due to meridional advection, the two fields are thus expected to be quite similar.
Figure 7d-f shows that the temperature anomalies from Figure 7a-c can be associated, to a substantial extent, with anomalous meridional displacements over large parts of Europe. More precisely, the fields reveal that air masses originating more from the south than usual mostly coincide with temperatures warmer than usual, whereas air masses originating more from the north than usual mostly coincide with temperatures colder than usual. For example, the elongated warm tongue of air on September 6 is represented quite well by anomalous northward displacement. Likewise, the cold anomaly over the Atlantic is captured well by the meridional displacement. The same applies to the other two examples. This suggests that temperature anomalies in the extratropics can generally be related quite significantly to meridional temperature advection.
The last statement is put on a more quantitative basis by computing the correlation between the corresponding fields. The squared correlation coefficient is given in the grey boxes in each panel; for all three dates it turns out to be well above 0.5, suggesting that at least half of the variance within the temperature anomaly field is described by the meridional displacement.
Of course, there are also differences between the meridional displacement and the temperature field: F I G U R E 7 Utility of the tracer method as a diagnostic for the occurrence of temperature anomalies. The figure shows horizontal cross-sections at a level of 850 hPa for three dates in September 2016. All variables are depicted in terms of normalised anomalies indicating the standard deviation from the two-month average of August and September 2016. R 2 denotes the squared correlation coefficient, which quantifies the quality of the correlation between the respective field and the temperature field. For the computation of R 2 , data points between 40 • N -60 • N and 50 • W -50 • E were considered. (a-c) Temperature from reanalysis. (d-f) Meridional displacement computed with a relaxation constant of −1 = 3 days; blue denotes air masses that have been anomalously displaced southward, whereas red denotes air masses that have been anomalously displaced northward. (g-i) Meridional wind from reanalysis; blue denotes anomalous northerly winds, whereas red denotes anomalous southerly winds. Panel (j-l) Three-day mean meridional wind from reanalysis; colours as in (g-i) for example, over France on September 6, where the meridional displacement anomaly extends further to the east than the temperature anomaly, or over England on September 12, where the displacement anomaly is stronger than the temperature anomaly. These differences are to be expected, since there are two other processes that induce parcel-based temperature changes and can thus contribute generally to the occurrence of temperature anomalies: adiabatic compression in vertical air motion and diabatic processes. Furthermore, differences may arise from the fact that air masses are significantly warmer or colder at their origin than expected from the zonal mean climatology; this is not taken into account in our simple argument.
In addition, our argument does not account for zonal temperature gradients. A more detailed and quantitative analysis along these lines is planned in the near future. Finally, note that the patterns of the temperature field show a closer resemblance to the meridional displacement than either the instantaneous meridional wind (Figure 7g-i) or the time-averaged meridional wind (Figure 7j-l). This illustrates, again, the merit of a Lagrangian method such as the tracer method. The tracer method provides Lagrangian information, namely the accumulation of meridional wind at the location of the parcel, as opposed to purely local Eulerian information,

Example 4: The Lagrangian versus the Eulerian mean meridional circulation
As mentioned before, the gridpoint-based nature of our method allows a straightforward evaluation and post-processing, as, for instance, in the computation of climatologies. To illustrate this particular feature, we computed the two-month time-averaged field of the meridional displacement (Figure 8a,c), based on three-hourly gridpoint-based fields of the meridional displacement of August and September 2016.
The time-averaged meridional displacement field reveals that the outcome of the tracer method is of truly Lagrangian nature, although the method itself is Eulerian from a technical point of view. The meridional displacement (Figure 8a) shows equatorward parcel displacements in the lower troposphere and poleward parcel displacements in the upper troposphere, with the displacement in the winter hemisphere being stronger than that in the summer hemisphere. The time-averaged meridional displacement thus has similar characteristics to the global Lagrangian mean circulation as diagnosed from the transformed Eulerian mean formalism (e.g., compare with Edmon et al., 1980, figure 6;Iwasaki, 1989, figure 4;Juckes, 2001, figure 4) or obtained from averaging the flow on isentropes (e.g., compare with Townsend and Johnson, 1985, figure 9;Juckes, 2001, figure 3;Pauluis et al., 2010, figure 1): It suggests a one-cell structure of the Lagrangian circulation in each hemisphere, with a stronger circulation in the winter hemisphere. The large values of the meridional displacement at the poles result from the fact that the poles constitute an upper and lower boundary for the meridional displacement, so that for high latitudes a displacement from lower latitudes is statistically much more likely; the large values at high latitudes must, therefore, be considered as a geometric effect.
In contrast to the Lagrangian mean, the conventional Eulerian mean of the meridional wind (Figure 8b), exhibits the familiar three-cell structure consisting of the Hadley cell, the Ferrel cell, and the polar cell. The difference between Figure 8a,b highlights the power and beauty of the tracer method: it provides a Lagrangian description of the flow through a simple Eulerian zonal average of a variable that carries Lagrangian information.
Finally, note that the two-month time average is far from being zonally symmetric even in the two-month time average (Figure 8c); rather, it exhibits poleward parcel displacement in northern Europe and equatorward parcel displacement in southern and eastern Europe as well as over the Atlantic, which is probably related to either large-scale stationary waves or frequently recurring weather systems. Again, a qualitative difference in comparison with the mean meridional wind (Figure 8d) is clearly visible.

Sensitivity with respect to the relaxation parameter
A crucial feature of our method is the relaxation term, which determines the time-scale of accumulation to be considered. Based on two examples from above, we now examine the sensitivity of our results with respect to the relaxation parameter .
For this purpose it is instructive to consider, first, the limit of strong relaxation, that is, T ≫ 1, where T denotes the time-scale on which the synoptic scale wind field changes. In this limit, Equation 7 can be approximated as Assuming furthermore (t − t 0 ) ≫ 1, the expression on the right-hand side becomes Translated to the Eulerian perspective, this result means that for strong relaxation the displacement field z (x, t) approaches the vertical wind field w(x, t) divided by . The other limit of interest is the limit of weak relaxation, that is, T ≪ 1. In this limit, Equation 7 turns approximately into Equation 5, meaning that the vertical displacement simply accumulates the vertical wind since the initial time t 0 . For integrations where this time interval is much longer than the synoptic time-scale, that is, t − t 0 ≫ T, we expect that the field z (x, t) does not provide useful information any longer, because the parcels undergo subsequent episodes of upwelling and downwelling in a more or less random fashion, and the displacement field loses any meaningful relation to the instantaneous synoptic-scale flow. Moreover, as we pointed out above, for such long integrations the whole concept of a "parcel" becomes questionable anyway. Therefore, we refrain from taking a closer look at this limit. Figure 9 confirms the above prediction for strong relaxation, based on the example presented in Section 5.2. In the case of strong relaxation ( −1 = 1 h, Figure 9b), the vertical displacement resembles, apart from a more diffusive behaviour, the vertical wind divided by (Figure 9a). Both fields diagnose, for example, distinctive upwelling over Iceland, western Russia, Turkey, and parts of the Atlantic, thereby differing significantly from the displacements obtained with −1 of no more than 12 h (Figure 9c).
For larger values of , the qualitative appearance of the vertical displacements in Figure 9 is rather insensitive to the specific choice of . As long as −1 is smaller than the typical Lagrangian synoptic scale, the tracer essentially samples coherent air streams with steady upward or downward motion. Hence, the vertical displacements do not differ much in their patterns, but mostly in their intensity; the latter is because a weaker relaxation constant accumulates displacement information for a longer time interval. In our example, this is roughly the case for values of −1 in the range between 12 h ( Figure 9c) and three days (Figure 9e). Here, the tracer method mainly detects shorter or longer fractions of the ascent within the cyclone's warm conveyor belt, typically lasting about two days (Wernli and Davies, 1997;Eckhardt et al., 2004). Even for −1 = 10 days (Figure 9f), the vertical displacement field still shows broadly the same structure. However, the overall colour turns towards blue, which means that most of the parcels have experienced upwelling rather than downwelling. Given that the tropopause acts approximately as a rigid lid, this behaviour is plausible from a statistical point of view, since air parcels residing in the upper troposphere must originate, in the long term, from further below.
The foregoing can be made more quantitative as follows. We computed the spatial correlation between displacement fields associated with different values of and show the result in Figure 10a. This plot quantifies the pattern correlation between the vertical displacement field for −1 = 3 days and that for various other values of . This plot generalizes the result from the previous figure by accounting for more than one pressure level and time step. Apparently, for very strong relaxation (small −1 ), the squared correlation coefficients are low, consistent with the earlier result that the vertical displacement fields for −1 = 1 h and −1 = 3 days differ substantially. For values of −1 in the range of 1-10 days, however, the squared correlation coefficients always exceed more than 80%, revealing a fairly strong correlation between the respective displacement fields. This indicates, again, that the displacement fields obtained sustain a similar overall structure for quite a large range of values of . While the overall patterns visible in the displacement fields in Figure 9 are nearly independent of over a broad range of values, this is not true for the typical magnitude of the fields (note the different contour intervals in the different panels). This fact is quantified further in Figure 10b, where we depict the mean absolute vertical displacement for the Northern Hemisphere midlatitudes as a function of −1 . Apparently, for small values of −1 (corresponding to short accumulation intervals) the magnitude of the vertical displacement increases with increasing −1 , but for larger values of −1 (corresponding to long accumulation intervals) this increase becomes more gradual or even saturates. As mentioned above, the general increase in magnitude with −1 is reasonable, since a larger −1 corresponds to a larger accumulation period. Understanding the precise shape of this curve and its dependence on altitude is beyond the scope of this article.
The meridional displacement exhibits a similar behaviour to the vertical displacement in terms of the sensitivity with respect to ( Figure S2). For very strong relaxation, the meridional displacement ( Figure S2b) represents, analogous to Equation 24, the meridional wind divided by ( Figure S2a). For all other values of in that figure, the meridional displacements mostly differ in their intensity but not in their structure. Again, the same conclusion can be drawn from Figure S3 in a more quantitative sense.
The best choice for the relaxation parameter may depend, amongst other things, on the synoptic situation, the region of interest, the variable under consideration, and the specific research question. However, the observed robustness of the patterns obtained with respect to indicates that the exact value is of secondary importance, at least in a qualitative analysis. Furthermore, the observed robustness with respect to suggests that the results should be rather insensitive to the exact choice of the function used as weighting kernel in Equation 7. This finding justifies, a posteriori, our decision to employ an exponential kernel that corresponds to simple Newtonian relaxation in the underlying differential equation.

SUMMARY AND CONCLUSIONS
In this article we proposed an Eulerian method for extracting Lagrangian information about the atmospheric flow. The method is based on the advection of passive tracer fields, which are subject to a judiciously chosen relaxation term. Mathematically, our method requires the solution of an advection-relaxation equation for a tracer field (x, t) on an Eulerian grid. The tracer field (x, t) can be defined through any scalar quantity that can be associated with an air parcel, such as altitude, pressure, latitude, or temperature. As output, the method provides the accumulated Lagrangian change (x, t) of the quantity considered. More specifically, the resulting field (x, t) represents the Lagrangian change of the quantity considered for those air parcels that happen to be located at the respective grid points at the time of interest. The relaxation term guarantees that the accumulation time is essentially restricted to a finite length, which in a typical application is much shorter than the total integration time. The finite accumulation time effectively corresponds to an exponential decay of (x, t) as one moves backward along the parcel's trajectory. This device allows us to run the integration in a continuous fashion without the need for repeated reinitialization and, thus, leads to an efficient and convenient method for specific applications.
We implemented a pseudospectral algorithm to solve the advection-relaxation problem on a three-dimensional, global, atmospheric domain. For the purpose of validation, we switched off the relaxation term and calculated so-called pseudotrajectories, which essentially provide the same information as traditional "real" trajectories. We then compared these pseudotrajectories with real trajectories obtained with LAGRANTO. As long as the length of the integration time is not much longer than the time-scale of the flow phenomenon of interest, our pseudotrajectories show close resemblance to trajectories from LAGRANTO; the only systematic difference is the fact that the pseudotrajectories exhibit a more diffusive behaviour due to the gridpoint-based nature of our algorithm. We interpreted this outcome as showing the correct implementation of our algorithm. Substantially larger differences appeared for time intervals that are considerably longer than the characteristic time-scale of deformation. In that case, the LAGRANTO trajectories provide spotty fields with sharp local gradients, which are familiar in the framework of "chaotic advection" (Aref, 1984;Pierrehumbert and Yang, 1992;Haynes, 1999). By comparison, our tracer algorithm yields considerably smoother fields owing to its inherent diffusion (both explicit and implicit). However, for such long time intervals the parcel approach becomes questionable, because a point cannot be taken as representative of an entire air parcel any longer; both approaches need to be considered as idealizations, and it remains unclear at this point to what extent either of them represents reality. This issue can be avoided by restricting the length of the trajectory to a finite time; in the framework of our tracer method this is effectively achieved by a finite value of . For synoptic-scale applications we suggest to use a value in the range −1 = 1-5 days, since this corresponds roughly to the typical time-scale of midlatitude synoptic-scale flow features.
We then demonstrated the intended use of our method with the help of a few intuitive examples related to synoptic-scale flow, using a finite relaxation constant > 0. The examples illustrated how our method can serve as a diagnostic tool to understand-not predict-the occurrence of certain atmospheric phenomena. For instance, in the first example we showed that observed patterns in cloud fraction correspond fairly well to the patterns in the vertical displacement obtained from the tracer method. By contrast, the correlation with the vertical wind was considerably lower. This result is consistent with the fact that adiabatic cooling by expansion is a major mechanism for cloud formation, and that the amount of cloudiness is related to the accumulated cooling experienced by an air parcel. In another example we pointed out that meridional displacements obtained from the tracer method show some resemblance to the distribution of temperature anomalies. This correspondence indicates that horizontal advection may be an important mechanism for the formation of temperature anomalies. We also briefly considered the patterns of the time-averaged meridional displacement, which indicate a one-cell structure in the meridional plane as opposed to the more familiar three-cell structure of the Eulerian mean meridional circulation, and thus demonstrated the Lagrangian nature of the diagnosed field.
Finally, we showed that the overall patterns of the vertical and meridional displacements in our examples were similar in a fairly broad range of values for , but at the same time the mean amplitude was quite sensitive to . This result suggests that the exact choice of the relaxation constant may be uncritical for a qualitative application of the method, but that care must be exercised in quantative applications. We are currently working to explore the utility of the method in a more quantitative sense. To be sure, a similar caveat applies to any trajectory method, where the results may depend partly on the length of the trajectory.
One major limitation of our method is the fact that it only provides Lagrangian information that is aggregated over a finite time interval. In this sense it provides less detailed information than can be achieved with trajectories. However, for certain applications aggregated information may be enough. These applications might include climatological analyses such as the Lagrangian characterization of heat waves and warm conveyor belts, or studies regarding the waviness of the jet stream. In such applications, our method may facilitate analysis by providing continuous, volume-covering Lagrangian information as direct model output. Due to its Eulerian, grid-based nature, the method is rather easy to handle and can be applied readily to large data sets. Moreover it can easily be implemented as an online tool in traditional Eulerian models, which provides advantages such as the use of highly resolved (in time) wind information (Gheusi and Stein, 2002;Miltenberger et al., 2013). In summary, we think that our tracer method is an interesting alternative to the traditional trajectory method for certain applications and it may stimulate further investigations of atmospheric phenomena from a Lagrangian perspective.

SUPPORTING INFORMATION
Additional supporting information can be found online in the Supporting Information section at the end of this article.