Ctraj has been applied to a new method of dynamical tracer reconstruction:
Download a pdf version of this page: description.pdf



Basic algorithm
Tranformed coordinates
Contour advection
Semi-Lagrangian advection

Basic algorithm

The only thing that CTRAJ does is integrate the following pair of ordinary differential equations:

where x=(x, y) is position, t is time and v = (u, v) is a known velocity field. All other codes are built up from this base. The codes take a set of initial points (which we will call the initial "field") and integrate them forward in time. The velocity field is gridded in a Cartesian space and linearly interpoled. The trajectories are integrated with a Runge Kutta fourth order scheme with a constant time step.

Return to top of page

Transformed coordinates

The codes are general enough that any velocity field can be integrated, so long as it is gridded in a Cartesian space. Since the velocities are encapsulated in container classes, adaptation to a non-gridded format or analytical computation of velocities will be quite straightforward. However the main purpose of these codes is to advect velocities over the surface of the Earth (e.g. output from GCMs, assimilation data, etc.) thus they must be transformed to a suitable projection. We use an azimuthal equidistant projection limited to one hemisphere. This has the dual advantage of eliminating the singularity at the pole and having relatively little variation in the area of each grid cell.

The transformations from standard longitude-latitude (constant diameter spherical polar) coodinates are as follows:

where r is defined as follows:

where R is the radius of the Earth and h is the pole upon which the transformation is centred:

The resulting space has the following metric:

while the reverse transformations are given as:

The velocity field transforms as follows:

where u is meridional wind and v is zonal wind. To run the simulation over the entire surface of the Earth we use two fields, one for the Northern hemisphere and one for the Southern. These fit together like two halves of an eggshell. Since equatorial crossings are relatively rare, we only check for them at intervals (typical when the field is output). For this reason, the velocity fields overlap slightly, as seen in the figure. The following transformations convert a coordinate to its equivalent representation in the opposite hemisphere:

Where the subscript h labels the hemisphere.

Return to top of page

Contour advection

Consider a blob of dye stirred by a moving fluid. To first order it can modelled by tracing its outlines and allowing these to move with the fluid. Contour advection is an adaptive method of modelling the contours or isolines of a passive tracer. It works by following the trajectories from points at intervals along the contour. To maintain the integrity of the curve, new points are added or removed as needed. We want the points to lie at roughly constant intervals along the curve even as it is stretched in many directions at once. Rather than use intervals of constant distance, we use fraction of arc:

where is the difference in path between two points and c is the radius of curvature and for a Cartesian space is calculated as follows:

The path along the curve s is first approximated as a set of piece-wise continous line segments. The resulting parametric representation of the contour is interpolated using a cubic spline. A cubic spline returns a set of second-order derivatives with which we can calculate the curvature: New points are added between pairs of points that trace out either too large an arc or too great a distance. Likewise, points are removed if they are too close to their neighbour.

Return to top of page

Semi-Lagrangian advection

A passive tracer simulation aims to solve the following differential equation, termed the advection equation:

where q is the tracer field in Lagrangian space, x0 is the Lagrangian coordinate and S is the source term. While S can accomodate other sources, such as chemical rates, here we concern ourselves with the diffusion term:

where D is the diffusivity tensor. In an Eulerian framework, this becomes:

Since, in the absence of diffusion, tracer values are constant along the trajectories, we can build a tracer simulation using only the trajectory integrator. The idea is to define an Eulerian grid and run back-trajectories from each grid point. Values of the tracer for the current field are interpolated from the previous field.

Because it doesn't matter in what order the trajectories are run (they are independent of the tracer field) and because the interpolation defines a series of coefficients mapping a tracer field from one time step to the next, we write these coefficients into a sparse matrix. The sparse matrix will have sides the size of the number of tracer grid points. A mapping from a two-dimensional tracer field to a one-dimensional vector must be defined.

Expressing this mathematically:

where q(ti) is the tracer field at the ith timestep and Ai is a matrix.

A semi-Lagrangian tracer simulation is not subject to the CFL criterion and is unconditionally stable. The interpolation step provides a small amount of diffusion, thus it can be almost eliminated by using an Eulerian time step (that controls the interpolation between fields, not the trajectory integration) that is as long as the integration time.

Return to top of page


The advection codes spend the majority of their time interpolating the velocity field. We have mitigated this in the case of the tracer simulation by writing coefficients to an array of sparse matrices thus the tracer simulation needs only be run once over a given period of time.


Return to top of page

Next: Manual