next up previous contents index
Next: 2. User's guide Up: $FILE Previous: Introduction   Contents   Index


1. Scientific documentation

1.1 Lagrangian and Eulerian descriptions of flows

1.1.1 Lagrangian description

The most natural and basic way to describe and model a fluid is to divide it into mesoscopic parcels, i.e., material systems contain a large amount of (microscopic) molecules but with dimensions much smaller than the smallest characteristic scale of the macroscopic flow. So, one can assume that a mesoscopic parcel is a close system governed by the fundamental laws of physics (the fundamental principles of mechanics and thermodynamics).

Such a parcel will be hereafter referred to as a Lagrangian parcel. The Lagrangian description of an evolving fluid (flow) consists in specifying for each Lagrangian parcel, and at any time, its location and any dynamical, thermodynamical, chemical, etc. , quantities that are of interest (hereafter fields with the generic notation $\alpha$).

This requires that each Lagrangian parcel has to be labelled so that its identity is unambiguous. A possible way is to label each parcel with its location, say $\vec{x}_0$, at an arbitrary reference time, $t_0$. The field $\alpha$ is thus written in Lagrangian description as a function of the initial location $\vec{x}_0$ of the considered parcel and time $t$, namely

\begin{displaymath}\alpha_{Lag} (parcel\: present\: at\: \vec{x}_0\: at\: time\: t_0\:,\: t),\end{displaymath}

or more simply


In particular, if the field $\alpha$ is the current location $\vec{x}$ of the Lagrangian parcels, the trajectory of parcel $\vec{x}_0$ is simply given by


1.1.2 Eulerian description

The Lagrangian description do not allow to easily formulate the equations of fluid dynamics, from which numerical models are derived. Moreover most measurement methods in experimental fluid mechanics and meteorology use sensors which are fixed in the study frame.

The Eulerian description consists in specifying the quantity $\alpha$ for the Lagrangian parcel coincidentally present at $t$ at the arbitrary location $\vec{x}$. The field $\alpha$ is thus written in Eulerian description as a function of space coordinates, $\vec{x}$, and time $t$, namely

\begin{displaymath}\alpha_{Eul} (parcel\: present\: at\: \vec{x}\: at\: time\: t\:,\: t),\end{displaymath}

or more simply


It should be stressed that the time evolution of $\alpha_{Eul}(\vec{x},t)$ do not refer to a given Lagrangian parcel but to an ensemble of parcels which have passed at the fixed place $\vec{x}$ (i.e., $\alpha_{Eul}(\vec{x},t+\Delta t)$ gives the quantity $\alpha$ for another parcel than $\alpha_{Eul}(\vec{x},t)$ does).

An interesting particular case of Eulerian field is when $\alpha$ is the initial location $\vec{x_0}$ of the coincidental Lagrangian parcels, namely


In words, this represents the location that the Lagrangian parcel currently present at $\vec{x}$ had at the reference time $t_0$.

1.1.3 Invertibility of the Eulerian initial location field

Since the Lagrangian parcels currently present at different places are all different and so had necessarily different initial locations, there is a one-to-one relationship between the initial and current locations of a the Lagrangian parcel. Hence the Eulerian initial-location field $\vec{x}_{0\:Eul}(\vec{x},t)$ can be inverted into the Lagrangian trajectory field $\vec{x}_{Lag}(\vec{x}_0,t)$.

In other words, the knowledge of $\vec{x}_{0\:Eul}(\vec{x},t)$ at any time since $t_0$ is sufficient to track any Lagrangian parcel along its trajectory, since this only requires at any time $t$ to make a general "identity control" in the whole space in order to find out where is the parcel that has a given value of $\vec{x}_0$.

1.1.4 Equivalence of the descriptions

The above field $\vec{x}_{0\:Eul}(\vec{x},t)$ and its invertibility are the keys to establish the equivalence of the Eulerian and Lagrangian descriptions, in terms of the information they are able to provide on the flow.

In effect, the value of $\alpha$ at any ($\vec{x},t$) is the Lagrangian value of $\alpha$ for the parcel coincidentally present at ($\vec{x},t$). The latter had the initial location $\vec{x}_{0\:Eul}(\vec{x},t)$. Hence


Conversely, the value of $\alpha$ at $t$ for parcel $\vec{x}_0$ can be known by reading at $t$ the Eulerian field $\alpha_{Eul}$ at the current point of the parcel trajectory, i.e., at $\vec{x}_{Lag}(\vec{x}_0,t)$. Hence


The latter expression in particular shows that the knowledge of the Eulerian field $\vec{x}_{0\:Eul}(\vec{x},t)$ at any time since $t_0$, together with an inversion method yielding $\vec{x}_{Lag}(\vec{x}_0,t)$, are sufficient to retrieve the full Lagrangian information from merely the Eulerian fields. The latter statement is the key of the Lagrangian analysis method presented in this document.

1.1.5 Physical variables Lagrangian evolution of physical variables

The evolution of the physical variables (e.g., in atmosphere sciences: water vapour content, potential temperature, potential vorticity, concentrations etc.) attached to Lagrangian parcels can be easily known by comparing the actual (Eulerian) value at given location and time, to the initial value for the parcel under consideration, i.e., the value read at the time origin and at the initial location of the parcel. More precisely, let us consider any physical quantity $\alpha$. The value of $\alpha$ at $t_0=0$ for the parcel located at $\vec{x}$ at time $t$ is given by

\alpha_0(\vec{x},t) = \alpha_{Eul}(\vec{x}_{0\:Eul}(\vec{x},t),t=0),
\end{displaymath} (1.1)

and hence the net change of $\alpha$ experienced by this parcel between $t_0=0$ and $t$ is
\Delta\alpha = \alpha_{Eul}(\vec{x},t)-\alpha_0(\vec{x},t).
\end{displaymath} (1.2) Physically-defined ensembles of Lagrangian parcels

The fields $\alpha_0(\vec{x},t)$ can also be used to select ensembles of Lagrangian parcels defined along a physical criterion at the time origin. Considering the ensemble of Lagrangian parcels defined at $t=0$ by $\{ \vec{x} \: \vert \: \alpha(\vec{x},t=0)=\alpha_{def} \}$, this ensemble later verifies to be $\{ \vec{x} \: \vert \: \alpha_0(\vec{x},t)=\alpha_{def} \}$. For example, this enables one to define true material surfaces that initially coincide with isosurfaces of poorly conserved quantities (e.g., potential temperature or potential vorticity in diabatic atmospheric flows).

1.2 Calculation of the initial location field in a numerical model

1.2.1 Method (continuous formulation)

Most fluid-dynamics numerical model are formulated with differential equations operating on Eulerian fields. It turns out from the previous subsection that the calculation of the Eulerian initial location field $\vec{x}_{0}(\vec{x},t)$ during the run of such a model allows (providing appropriate post-processing) to exhaustively analyse the modelled flow from the Lagrangian point of view.

This calculation can be easily carried out since it is only based on the transport of passive scalar variables - which is a very common task in numerical modelling. The method consists, at an arbitrary reference time $t_0$ of the model run ($t_0$ can be either at the model start or also later during the run), in initializing three fields of passive scalars with the three space coordinates, then driving them along the modelled flow.

More precisely, if Cartesian coordinates are chosen, three scalar fields $x_0$, $y_0$, $z_0$ (hereafter referred to as initial coordinates) are defined this way at $t_0$:

$\displaystyle x_0(\vec{x},t_0)$ $\textstyle =$ $\displaystyle x,$ (1.3)
$\displaystyle y_0(\vec{x},t_0)$ $\textstyle =$ $\displaystyle y,$ (1.4)
$\displaystyle z_0(\vec{x},t_0)$ $\textstyle =$ $\displaystyle z.$ (1.5)

(In the following only Eulerian fields are considered so that the subscript Eul is no longer kept.)

Since these fields are passive scalars driven by the flow, their evolution is governed by the equation

\frac{\partial \vec{x}_0}{\partial t} = - \vec{u}\cdot\overrightarrow{\nabla}\,\vec{x}_0,
\end{displaymath} (1.6)

where $\vec{u}$ is the velocity field and $\vec{x}_{0}(\vec{x},t)$ is the initial location vector-field formed with the components $x_0$, $y_0$, $z_0$.

At the lateral boundaries of the model, the initial coordinates are prescribed depending on the out- or inflow condition. Where outflow occurs, parcels are definitively eliminated from the model. In the inflow case, new parcels enter the model. The initial coordinates of the inflowing parcels are set according to table 1.1. With such values, any parcel which has entered the model during the simulation can be easily identified since one of its horizontal initial coordinates exceeds the model-domain dimension. Moreover, for such a parcel one is able to know where and when it has entered the model. For example, for a model grid-point let $x_0=\xi <
x_{min}$, $y_0=\upsilon$, $z_0=\zeta$. We can conclude that the parcel has entered the model at the point $(x_{min},\upsilon,\zeta)$ and at the time $t=(x_{min}-\xi)/c_{in}$ (where $c_{in}$ is the velocity-like constant defined in table 1.1).

Table 1.1: Values of the 'initial' coordinates for parcels entering the simulation domain at time $t$. $x_{min}$, $x_{max}$, $y_{min}$, $y_{max}$ are the abscissas and ordinates of the lateral boundaries. $c_{in}$ is an arbitrary constant with the dimension of a velocity.
Lateral Boundary Initial Coordinates
$x=x_{min}$ $y_0(x_{min},y,z,t)=y$
$x=x_{max}$ $y_0(x_{max},y,z,t)=y$
$y=y_{min}$ $y_0(x,y_{min},z,t)=y_{min}-c_{in}t$
$y=y_{max}$ $y_0(x,y_{max},z,t)=y_{max}+c_{in}t$

1.2.2 Discretized formulation and implications Discretized formulation

The previous considerations are valid in continuous formulation. However numerical models cannot be formulated this way but rather operate on space- and time-discretized fields (actually 3D matrices). In discretized formulation, the initial coordinates $x_0$, $y_0$, $z_0$, are initialized at $t_0$ with the coordinates of the model grid-points then experience all the modelled transport processes:

The former corresponds simply to the discretized form of the right-hand term of Equation 1.6, while the latter two processes account for the transports by unresolved (subgrid) fluid motions by means of more or less refined parametrizations. In addition, the initial coordinate tracers will also be influenced by any kind of numerical diffusion. Sub-grid mixing

Sub-grid processes as well as numerical diffusion have the general effect to mix the distributions of passive scalars. In the considered special case, where passive scalars are not chemical stuff concentrations (or any other additive quantity) but rather are used to transport the initial coordinates, the mixing has to be interpreted to some extent as a progressive loss of identity for the Lagrangian parcels. For example, the equal mixing of two parcels '1' and '2' originating from altitudes $z_1$ and $z_2$ does not result in a parcel from the altitude $(z_1+z_2)/2$, although this is what is produced by the model.

The mixing of the initial coordinates makes nevertheless sense. While using partial derivative equations (even in discrete form) we indeed implicitly assume that any field, including the initial coordinates tracers, remains continuous. For the example considered, the mixing of parcels 1 and 2 in the same grid cell thus implies that an infinity of parcels with all the intermediate values of the initial altitude between $z_1$ and $z_2$ have also been transported into the considered grid cell. Thus, the information $(z_1+z_2)/2$ for this grid cell is correct, since the mean value of the initial altitude for the parcels present within the grid cell is actually $(z_1+z_2)/2$. However the most precise information for this grid cell should be "parcels from everywhere between $z_1$ and $z_2$", or equivalently "parcels from everywhere around $(z_1+z_2)/2$ within a radius $\vert z_1-z_2\vert$".

This suggest that the information given by the initial coordinates in the context of a discretized formulation should be ideally completed by some kind of accuracy radius indicating the dispersion of the initial locations of all the parcels contained in a grid cell, around the barycenter of these initial locations. This radius would have as initial value the half-size of the grid cell, then increase during the model run, especially in regions of active mixing such as in deep vertical transports, boundary layers, etc. It is also clear that the accuracy radius would be all the greater that the mixed parcels are initially remote. Unfortunately the growth of such an accuracy radius is not so easy to quantify and the authors' reflexion still remain at the stage of this qualitative discussion1.1. Discussion with respect to trajectory methods

It must be however stressed that this progressive loss of identity is not specific of the initial-location method for Lagrangian analyses, but also exists with any trajectory-based method - in a more hidden way. Any inaccuracy in the integration process of a trajectory path make that the considered parcel actually changes every integration step. Hence the parcel considered at the end of the calculated trajectory might be quite remote from the parcel considered at the start, although the calculated trajectory has the appearance of a single, well-defined, parcel path.

Moreover most trajectory methods have the further drawback to ignore the subgrid transports whereas they are taken into account with the initial coordinates method - as well as parametrizations can do.

1.2.3 Post-processing with multiple reinitializations during the model run Run-time reinitializations of the tracers

As discussed above, there may be a great uncertainty on the identity of parcels after they have experienced strong mixing processes. Nevertheless one could be interested in the later evolution of these parcels. For example in the atmosphere, deep convection may transport pollutants from the boundary layer which then chemically evolve in the upper troposphere.

In such cases a possible way to restore a reliable information on parcel identity is to reinitialize the tracer fields $x_0$, $y_0$, $z_0$ along Eq (1.3-1.5) once or more during the model run. This has the added advantage to let the user have the choice in the time-origin of the studied evolution among the various reinitialization times.

Periodically reinitialized tracers will thus give a certain amount of evolution fragments over whiles shorter than the whole run. In fact the Lagrangian information over longer periods than the reinitialization frequency is actually not lost and can be retrieved by simple post-processing (this will be shown in the next paragraph). Therefore reinitializations can be done as often as needed to keep the tracers fields "clean" enough (i.e., not too mixed). Catenation of path fragments

Let $t_0=0<t_1<..<t_n$ be the times when the initial coordinate fields are reinitialized during the model run ($t_0$ corresponding to the model start). Hereafter $\vec{x}_{t_i}(\vec{x},t)$ denotes the 'initial' location at reinitialization time $t_i$ of the parcel currently present at $(\vec{x},t)$ (i.e., the subscript $t_i$ refers to the time origin)1.2.

Let us consider $t_i$ as the time origin for the Lagrangian evolution that is under interest. For $t_{i+1}<t$, the initial coordinate fields have been reinitialised at least once since the time $t_i$. $\vec{x}_{t_i}(\vec{x},t)$ can however be retrieved by joining together fragments of Lagrangian motions. With $t> t_j >
t_i$, the method consists in reading $\vec{x}_{t_j}$ at time $t$ and location $\vec{x}$, then $\vec{x}_{t_{j-1}}$ at time $t_j$ and location $\vec{x}_{t_j}$, and so on, until reading $\vec{x}_{t_i}$ at time $t_{i+1}$ and location $\vec{x}_{t_{i+1}}$. This is sketched in figure 1.1.

Figure 1.1: Sketch for the catenation process or backtrajectories.
\includegraphics[angle=0, width=\textwidth, keepaspectratio=true,clip=true]{figdocu/catenation.eps}

More precisely, let us consider the Lagrangian parcel at time $t$ at the location $\vec{x}$ and denote $\vec{X}(t_j)$ ( $i \leq j \leq i+p$) the successive locations of the parcel at times $t_i<..<t_{i+p}$. If $\vec{X}(t_j)$ is known, the previous location of the parcel $\vec{X}(t_{j-1})$ can be retrieved by reading, at time $t_j$ and location $\vec{X}(t_j)$, the initial location vector transported with the parcel since $t_{j-1}$, that is:

\vec{X}(t_{j-1})=\vec{x}_{t_{j-1}} \left( \vec{X}(t_j),t_j \right).

If $t_j < t \leq t_{j+1}$ with $i < j < n$ (or $t>t_n$ if $j=n$), the location of the parcel $(\vec{x},t)$ at time $t_j$ can thus be retrieved along the following algorithm:

$\displaystyle \vec{X}(t_j)=\vec{x}_{t_j}(\vec{x},t),$     (1.7)
$\displaystyle \vec{X}(t_{j-1})=\vec{x}_{t_{j-1}} \left( \vec{X}(t_j),t_j \right),$     (1.8)
$\displaystyle \ldots$      
$\displaystyle \vec{X}(t_i)=\vec{x}_{t_i} \left( \vec{X}(t_{i+1}),t_{i+1} \right),$     (1.9)

and finally


In particular, one may choose $t_i=t_0=0$ and obtain the initial coordinates field that would have been calculated without reinitialization after the model start. Backward trajectories and along-trajectory evolutions

The algorithm (Eq 1.7-1.9) can be directly used for the computation of individual backward trajectories, since it gives the successive locations of a chosen Lagrangian parcel. Further, the value of any physical quantity $\alpha$ along the Lagrangian parcel trajectory can be immediately obtained from the Eulerian field of $\alpha$ as

\alpha_{traj}(t_j)=\alpha_{Eul} \left( \vec{X}(t_j),t_j \right) .

1.3 Specific implementation in MesoNH

1.3.1 In program MODEL

One-model run

The three passive scalar fields $x_0$, $y_0$, $z_0$, (located at the mass-points of the Arakawa-C grid1.3) are initialized along Eq.1.3-1.5 just before the beginning of the main temporal loop (routine ini_lg.f90 called by ini_modeln.f90). This initialization occurs at least at the simulation START, and also possibly at every RESTART of a new simulation segment (see option "LINIT_LG=.TRUE." in chapter2). For that purpose the horizontal conformal coordinates $\hat{x}$, $\hat{y}$ and the true altitude $z$ are interpolated at the mass-points (see MesoNH's Scientific Documentation - Book 1, §3.1 and 4.2).

During the model run, these fields are treated exactly as any other passive scalar variables by the model dynamics (i.e., they experience all explicit as well as parametrized transports)1.4. Distinction is only made in order to render the turbulent mixing of the Lagrangian scalars desactivable independently from that of the scalar variables (see option LNOMIXLG in chapter2). Other dynamical options are common for all scalar variables.

At the lateral boundaries under inflow condition, $x_0$, $y_0$, $z_0$ are prescribed according to Table 1.1 where: $x$, $y$ are again replaced by the model coordinates $\hat{x}$, $\hat{y}$; $z$ is the true altitude of the mass-point; $t$ is the while since the last (re-)initialization of the Lagrangian tracers. The inflow velocity-like constant $c_{in}$ is set to 10 m/s (parameter XLGSPEED set in modd_lg.f90). The lateral inflow prescription is made in routine setlb_lg.f90.

Run with grid-nesting

In the case of simulation with nested models, the initial coordinates at the boundaries of the most outer model only are prescribed as above.

The initial coordinates at the boundaries of inner models are given by the "father"1.5model for inflow conditions, and evacuated into the father model for outflow conditions. The boundaries of the inner models are thus fully "two-way permeable" for the Lagrangian variables (as well as for any other scalar variable).

1.3.2 In program DIAG

The main purpose of program DIAG concerning the Lagrangian tracers is to perform catenations of path fragments with the algorithm Eq.1.7-1.9 (also Fig.1.1).

The process (compute_r00.f90, called by diag.f90) simply consists in reading data in earlier FM files at locations given by the $\vec{X}(t_k)$'s. In general these do not coincide with grid points. The discretized formulation of the algorithm therefore requires spatial interpolations between the grid points. A simple linear interpolation scheme is used in the present version of the system. The numerical code has been tested by comparing catenated and non-catenated Lagrangian tracers fields covering evolutions over the same while, with an excellent agreement found.

It also may occur during the process that one of the $\vec{X}(t_k)$'s is located outside the model domain. In that case, the process is stopped and we have adopted, for simplicity, to flag the value $\vec{x}_{t_i}(\vec{x},t)$ as "parcel not present in the model domain at $t_i$" without any other information (the point $(\vec{x},t)$ will appear as data empty in plots of the fields $x_{t_i}$, $y_{t_i}$ and $z_{t_i}$, or in any other plots requiring this information). The flag value is -1.E+11 (ZSPVAL in compute_r00.f90).

1.3.3 In program DIAPROG

For plots on iso-$z_0$ surfaces or other iso-surfaces, the same algorithm is used as for plots on isobars, isentropes, iso-PV, etc. . In case of multi-valuation over a given vertical, the point of highest altitude is retained.

The computation of backtrajectories is performed with the algorithm Eq.1.7-1.9 as in DIAG. For a given trajectory, the iteration is stopped as soon as a parcel is found to originate from outside the model domain at the previous step (or of course as soon as the information of the earliest MesoNH file is used).

For the computation of streamlines (case of a stationary flow) the same algorithm is used recursively with the initial coordinate fields of a single MesoNH file (since these fields are assumed to be stationary). The backward computation of the streamline is stopped after 100 iterations or as soon as as the parcel exits the model domain.

For details, please contact Jacqueline Duron ( or Joël Stein

next up previous contents index
Next: 2. User's guide Up: $FILE Previous: Introduction   Contents   Index
Jean-Pierre CHABOUREAU 2005-05-13