4.2 The Trajectory Equation




The advection of a particle or puff is computed from the average of the three-dimensional velocity vectors at the initial-position P(t) and the first-guess position P'(t+Δt). The velocity vectors are linearly interpolated in both space and time. The first guess position is

  • P'(t+Δt) = P(t) + V(P,t) Δt
and the final position is
  • P(t+Δt) = P(t) + 0.5 [ V(P,t) + V(P',t+Δt) ] Δt.

The integration time step (Δt) can vary during the simulation. It is computed from the requirement that the advection distance per time-step should be less than the grid spacing. The maximum transport velocity is determined from the maximum transport speed during the previous hour. Time steps can vary from 1 minute to 1 hour and are computed from the relation,

  • Umax(grid-units min-1) Δt (min) < 0.75 (grid-units).

The integration method is very common (e.g. Kreyszig, 1968) and has been used for trajectory analysis (Petterssen, 1940) for quite some time. Higher order integration methods will not yield greater precision because the data observations are linearly interpolated from the grid to the integration point. Trajectories are terminated if they exit the meteorological data grid, but advection continues along the surface if trajectories intersect the ground.

In the horizontal, the integration of the position vector is done in grid units, while in the vertical, a normalized sigma coordinate system is used, where sigma is defined by

  • σ = (Ztop - Zmsl) / (Ztop - Zgl),
where Zmsl is the height in mean sea level to be converted to sigma, gl is the height of ground level, and top defines the scaling height of the model top, the height at which the internal sigma surfaces go flat relative to terrain. The default internal scaling height is set to 25 km. This value is used in conjunction with the default model top setting (10 km for most PBL applications) in the CONTROL file which sets the vertical limit of the internal meteorological grid. If calculations are not required above a certain level, fewer meteorological data are processed thus speeding up the computation. Trajectories will terminate when they reach this level and mass particles will reflect back into the domain. A secondary use of the model top parameter in the CONTROL file is to set the model's internal scaling height to a different value but only if the entry exceeds 25 km.

The animation shows the sequence of model computations of a single trajectory (red) from 750 m AGL with the associated 900 hPa height contours (black). Above the boundary layer, when the winds become geostrophic, the trajectory path is parallel to the height contours.

It should be noted that the variable time step can at times lead to what might appear as inconsistencies in the calculation results. This may be particularly evident when there are strong gradients in the wind field (in space or time) that may not be adequately represented by the resolution of the gridded meteorological data. For instance, when calculating two or more trajectories at the same time compared with their calculation individually, because the maximum wind speeds for each trajectory may differ, the time steps could be different for the individual calculations, resulting in slightly different trajectories. This is not an error but should be treated as a cautionary note that less confidence should be placed on an individual trajectory. Some of the issues relating to trajectory uncertainty will be examined in more detail in subsequent sections.