10.4. Lagrangian models

Lagrangian models are based on the idea that pollutant particles in the atmosphere move along trajectories determined by the wind field, the buoyancy and the turbulence effects. Calculating these trajectories means an ODE instead the original dispersion PDE problem that is computationally easy to solve. The final distribution of many particles gives a stochastic estimation of the concentration field. The models either estimate the particle as a single drifting point, and the final distribution of numerous particles is used to estimate concentration fields (trajectory models), or assume a Gaussian dispersion inside each particle and the final concentration field is given as a superposition of these Gaussian distributions (puff models).

The Lagrangian approach has the advantage against Eulerian models that no grid is used for computation, thus spatial discretization errors like numerical diffusion are avoided. The result can be interpolated to any grid which means that the model error is independent of the output resolution. Deposition and radioactivity can be taken into account as a time-dependent decay term within each particle (Stohl et al., 2005), however, the treatment of chemical reactions can only be carried out if concentration field is available on a grid. Therefore reactive materials require an interpolation on a grid at each timestep that introduces spatial truncation errors and numerical diffusion in the model.

Lagrangian models are exceptionally efficient close to the source, where girded computations would require a very fine resolution to handle large gradients. On the other hand, long-range simulations require progressively more particles that causes a much faster increase of computational cost with the domain size than for the Eulerian models. To unite the advantages of both approaches, nested models have been developed that use a Lagrangian model near the source, but interpolate the variables on a grid to perform Eulerian dispersion and chemistry simulations for larger distances.

10.4.1. Calculation of trajectories

The trajectory equation for a single particle is written as an ordinary differential equation



where is the position of the particle, is the grid-scale particle speed including advection, settling and buoyancy, while is the turbulent wind fluctuation vector. Some models present a third term to describe subgrid mesoscale wind fluctuations (Stohl et al., 2005). The turbulent fluctuation term can be estimated as a random walk in a viscous fluid, which is described by the Langevin equation



where TL is a Lagrangian integration time step, is the vertical turbulent velocity fluctuation and dW represents a white noise process with mean zero and variance dt. A similar equation can be written for the horizontal fluctuations. The last component of equation (10.21) describes a random walk, while the first component is a memory term that represents autocorrelation. The TL Lagrangian timescale is a key parameter and is often given explicitly, or computed from velocity fluctuations (Stohl et al., 2005). At large TL values, autocorrelation becomes low, and the turbulent diffusion can be estimated as an uncorrelated random walk, that gives a reasonable estimation only for long distance simulations. However, in some models only this simple random walk process is implemented. Velocity fluctuations are computed from eddy diffusivities, or defined through parameterizations using the Monin-Obukhov-theory (Stohl et al., 2005).

The easiest way to handle anisotropic turbulence is to solve separated Langevin equations for each wind component using different wind fluctuations and Lagrangian timescales (Stohl et al., 2005). This way we neglect cross-correlations between turbulent fluctuations that proved to be a reasonable assumption for meso- and macroscale simulations. The Langevin equation (10.21) assumes Gaussian turbulence, thus it is only valid in case of isotropic turbulence and neutral stratification. Several approaches were presented to modify equation (10.21) in order to take into account mesoscale horizontal wind shear, buoyancy and convective boundary layer turbulence effects (McNider et al., 1988, Luhar and Britter, 1989, Pozorski and Minier, 1998, Stohl and Thomson, 1999).

When solving the trajectory equation numerically, temporal discretization of equations (10.20-10.21) is carried out. Equation (10.21) is usually solved with a first-order scheme (Stohl et al., 2005), however, applications of higher order Runge-Kutta methods are also present. Discretization of equation (10.21) is presented in Pozorski and Minier (1998). We note that for numerical representation of the white noise process in equation (10.21), an efficient random number generator is also needed that can significantly increase the computational cost of the model.

10.4.2. Puff models

Puff models are an intersection between Gaussian and Lagrangian dispersion models. They hold the assumption that the concentration pattern can be well described with a Gaussian distribution; however, the centreline of a plume is not a straight downwind direction, but a Lagrangian trajectory (Figure 10.4). This way they can take into account temporal and spatial wind changes. Puff models split the released mass of pollutants into small units, “puffs”, and then compute the trajectories of all the puffs in a Lagrangian way (through equation 10.19), but keep a Gaussian concentration pattern inside each puff. The final concentration field is given as a superposition of all the puffs’ concentration distributions:



where is the source term, N is the number of puffs, (xk, yk, zk) is the position of the kth puff and is the x-directional deviation of the Gaussian distribution inside the k-th puff. Equation (10.22) is very similar to the equation (10.19) with the difference that at the puff models, a sum of many components is calculated. Reflection from the ground and the inversion layer can be added to the equation in a way presented in equation (10.19), however, it describes only the reflection inside the puff. In a complete 3D model, reflection of the trajectory of the puff also has to be computed. Deposition and decay can be treated within each puff in the same way as in the Gaussian models. The effect of buoyancy can be estimated with the plume rise parameterization, or with both trajectory and in-puff treatment. Taking a time-dependent source strength in Q multiplied with the release frequency of puffs, the temporal variability of the source can be represented in the model.

Comparison of Gaussian plume and puff models.

Figure 10.4: Schematic representation of Gaussian plume and puff models. Puff models still estimate a Gaussian dispersion, but are able to take into account temporal and spatial wind changes.

Puff models calculate the effect of turbulence in two different ways: with a stochastic random-walk approach in the trajectories of the puffs (equation (10.20)); and through the deviation of a normal distribution inside each puff (equation (10.22)). It can be interpreted as a separated treatment of large-scale and sub-puff-scale turbulence. In long-range simulations, as the puffs grow due to diffusion while travelling along their trajectories, initially small puffs will reach a size where the wind shear becomes significant. To deal with this issue, developers introduced a puff-splitting technique that mass-consistently splits a puff into two parts if it reaches either a maximum size or a minimum sub-puff scale wind shear.

There are models like RAPTAD that encounter with both trajectory and puff turbulence. However, most of the softwares simplify the calculation with neglecting one of the turbulent processes. As the simplest approach, CALPUFF and RIMPUFF compute Gaussian turbulence deviations inside the puffs, but neglects large-scale turbulence and only advective trajectories are given. A more sophisticated solution is presented in HYSPLIT, where stochastic random-walk motion of the puff is used only in vertical direction, while horizontally a two-dimensional in-puff Gaussian distribution is assumed. The NAME model can be run both in point-particle (only random walk) and puff mode. We note that if the in-puff Gaussian turbulence is neglected (deviation is zero), the puff can be treated as a single point particle, which leads to the trajectory models.

Puff models were the first tools to simulate long-range atmospheric dispersion processes, as they were developed from existing Gaussian models, but were able to encounter with changing wind and emission data. In Europe, the Danish models, DERMA and RIMPUFF are widely used in risk management for case and sensitivity studies. RIMPUFF is part of the RODOS (Realtime Online Decision Support) system, the European framework for nuclear security. In the US, the CALPUFF was designated (besides AERMOD) as a preferred model of the Environmental Protection Agency. CALPUFF was successfully used worldwide in various environmental protection and public health investigations as well as for regulatory purposes. The RAPTAD model was developed for mesoscale simulations around complex terrain, and was also applicated for nuclear risk management and microscale dispersion problems.

The HYSPLIT model, developed by the NOAA Air Resources Laboratory, is one of the most widely used dispersion softwares. It provides an easy-to-use online interface for single trajectory simulations in order to give a fast estimation of atmospheric dispersion pathways or source regions. It also offers a mixture of a vertical trajectory and a horizontal puff model to determine concentration levels.

10.4.3. Trajectory models

Trajectory models are based on a stochastic simulation of a large number of drifting point-particles that represent a fraction of the released mass of pollutants (Figures 10.5, 10.6, and 10.7). Instead of using the Gaussian assumption, they solve all scales of turbulence explicitly above the TL Lagrangian timescale using the (10.20 and 10.21) trajectory and Langevin equations. They provide a sophisticated treatment of turbulent diffusion without solving partial differential equations, however, to obtain reliable results from a stochastic model, the calculation of a large number of single trajectories is needed that increases the computational cost. When the final particle distribution is available, either the concentration field can be given with a simple box-counting model, or cluster analyses can be performed to determine the most affected areas.

Air mass trajectory in the HYSPIT model.

Figure 10.5: Air mass (or single particle) trajectory in the HYSPIT model.

Dispersion of particles in the HYSPIT model.

Figure 10.6: Dispersion of particles in the HYSPIT model.

Averaged concentration of an air pollutant.

Figure 10.7: Averaged concentration of an air pollutant using particles (from Figure 10.6) in the HYSPIT model.

Modern computers allow the computation of millions of long-range trajectories within a reasonable time that made trajectory models the state-of-the-art tools of regional to global scale atmospheric dispersion simulations. The British NAME and the open-source Austrian FLEXPART are well validated, popular models successfully used worldwide for dispersion simulations. They provided reliable and fast information for authorities and scientists in recent air pollution episodes like the 2001 airborne foot-and-mouth disease epidemic (Mikkelsen et al., 2003), the Eyjafjallajökull eruption in 2010 (Dacre et al., 2011) or the Fukushima nuclear accident in 2011 (Stohl et al., 2011).

Trajectory models are often used in backward mode in order to identify source regions or provide an estimation of the mass of the released material. While it is easy to compute backward advective trajectories using an archived NWP wind field, the irreversibility of turbulent diffusion and deposition has to be carefully treated. In backward simulations, the final distribution obtained from the random walk of many trajectories represents sensitivity, and not a concentration field (Stohl et al., 2005). In other words, while in forward mode random walk describes a physical process (i.e. turbulent diffusion), in backward mode, random walk is used to describe the uncertainty of the trajectory caused by turbulent diffusion and deposition. The output of backward simulations is a source-receptor sensitivity field that provides probabilities of possible source locations. Once the most probable source locations are known, the estimation of the mass of released material has to be carried out (Stohl et al., 2011). It is usually performed with Monte–Carlo methods calculating multiple groups of forward trajectories with different source terms, or Eulerian inverse modelling.

We note that a simplification of trajectory models also exists that compute single trajectories to determine the main direction of the dispersion. These models don’t estimate concentration, but provide a fast and easy-to-use first-guess of the affected areas. Single trajectory calculations can be easily carried out using a numerical weather prediction model’s wind field. A single trajectory gives few information of the turbulence, thus the method provides reliable results only in advection dominated, large distance and free atmosphere simulations. Near-source or boundary layer modelling requires large number of trajectories to estimate the effect of turbulent diffusion.