Evaluation of oceanic and atmospheric trajectory schemes in the TRACMASS trajectory model v 6 . 0

Three different trajectory schemes for oceanic and atmospheric general circulation models are compared in two different experiments. The theories of the trajectory schemes are presented showing the differential equations they solve and why they are mass conserving. One scheme assumes that the velocity fields are stationary for set intervals of time between saved model outputs and solves the trajectory path from a differential equation only as a function of space, i.e. "stepwise stationary". The second scheme is a special case of the "stepwise-stationary" scheme, where velocities are assumed constant between GCM 5 outputs, it uses hence a "fixed GCM time step". The third scheme uses a continuous linear interpolation of the fields in time and solves the trajectory path from a differential equation as a function of both space and time, i.e. "time-dependent". The trajectory schemes are tested "off-line", i.e. using the already integrated and stored velocity fields from a GCM. The first comparison of the schemes uses trajectories calculated using the velocity fields from a high resolution ocean general circulation model in the Agulhas region. The second comparison uses trajectories calculated using the wind fields from an atmospheric reanalysis. The 10 study shows that using the "time-dependent" scheme over the "stepwise-stationary" scheme greatly improves accuracy with only a small increase in computational time. It is also found that with decreasing time steps the "stepwise-stationary" scheme becomes increasingly more accurate but at increased computational cost. The "time-dependent" scheme is therefore preferred over the "stepwise-stationary" scheme. However, when averaging over large ensembles of trajectories the two schemes are comparable, as intrinsic variability dominates over numerical errors. The "fixed GCM time step" scheme is found to be less 15 accurate than the "stepwise-stationary" scheme, even when considering averages over large ensembles.


Introduction
The Lagrangian view of the ocean and atmospheric circulation describes fluid pathways and the connectivity of different regions, which are not readily obtained from a Eulerian perspective.Lagrangian studies often require trajectory calculations using some algorithm that transforms the Eulerian velocity fields, e.g.winds or currents, into trajectories.Although observed velocities can be used, it is much more common to use velocities simulated by a general circulation model (GCM).The purpose of this work is to test the different schemes used in the TRACMASS trajectory model (version 6.0), here named the fixed GCM time step (Blanke and Raynaud, 1997;Döös, 1995), stepwise stationary (Döös et al., 2013) and time-dependent (de Vries and Döös, 2001) schemes.These schemes have previously only been tested using highly idealised velocity fields.Here, we will test the velocity fields simulated by comprehensive GCMs for both the ocean and atmosphere.
The TRACMASS trajectory model (Jönsson et al., 2015) has been continuously updated through the years since it was first introduced by Döös (1995).Version 6.0 represents the latest version, which includes the ability to run TRAC-MASS with the time-dependent scheme by de Vries and Döös (2001) on GCM fields.TRACMASS now also supports Published by Copernicus Publications on behalf of the European Geosciences Union.
many different types of vertical coordinates used in atmosphere and ocean GCMs.The code has also been made more structured and user friendly.
The original feature of TRACMASS and the related Ariane model (Blanke and Raynaud, 1997) is that they solve the trajectory path through each model grid cell with an analytical solution of a differential equation, which depends on the velocities on the faces of the model grid box.This is different from iterative schemes such as the commonly used fourthorder Runge-Kutta (RK4).The TRACMASS schemes have many advantages, e.g.mass conservation within the grid cell in the same way as the GCM itself, as well as fast trajectory computation.Furthermore, as the solution to the differential equation is unique, trajectories can be calculated forward in time and subsequently backward in time to arrive at exactly the original position.This makes it possible to trace the origins of water or air masses as long as stochastic parameterisations (see Döös and Engqvist, 2007) are not activated.
The first trajectory scheme tested here, the fixed GCM time step, is strictly only valid for stationary velocity fields.It can, however, be used with time-varying velocity fields by dividing the time between GCM outputs into intermediate steps and assuming velocities are stationary during the step.The velocities in an intermediate step are found by linear interpolation between two GCM outputs and hence named stepwise stationary.However, using intermediate steps increases the computational cost.The time-dependent scheme does not assume that the fields are stationary and uses instead continuous bilinear interpolation both in space and time.
The fact that the stepwise-stationary scheme uses stepwise-stationary velocities is logical when the scheme is used online, i.e. integrated into a GCM and thus having the same time step as the GCM itself.When the scheme is used offline, i.e. separately from the GCM and after the velocity fields have been stored, the time step is the time between two GCM outputs, which typically is a much longer period than the GCM time step.As the stepwise-stationary scheme assumes that velocities are constant during the time step of the trajectory scheme, processes faster than the GCM output frequency are lost.
An alternative to the stepwise-stationary scheme was introduced by de Vries and Döös (2001), where the trajectory solution was not only solved analytically in space, as was done by Blanke and Raynaud (1997) as well as Döös (1995), but also analytically in time between the GCM outputs.This leads to a more complex differential equation to be solved and integrated as the trajectory progresses through space and time (Döös et al., 2013).The advantage of this time-dependent scheme by de Vries and Döös (2001) is that it does not require any intermediate time steps between the model output times and can instead be integrated analytically between the GCM outputs.This method contrasts the fixed GCM time step scheme by Blanke and Raynaud (1997) and the stepwise stationary by Döös et al. (2013) as well as schemes such as the Euler forward or RK4 methods (Butcher, 2016;Fabbroni, 2009), where the trajectories are integrated forward in time with as short time steps as possible.A comprehensive review of different trajectory codes as well as the fundamental kinematic framework behind these can be found in van Sebille (2016).
In Sect.2, we describe the three different trajectory schemes and how they are integrated in time in both ocean general circulation models (OGCMs) and atmospheric general circulation models (AGCMs).In Sect.3, we test the three trajectory schemes with two different velocity fields, one from an OGCM and one from an AGCM, using various statistics.This study is concluded in Sect. 4 with a summary and discussion of the main results of the trajectory schemes and their tests.

Trajectory scheme theory
The trajectory schemes used in TRACMASS are all mass conserving but make different assumptions regarding the time evolution of the Eulerian velocity and pressure fields.The schemes rely on the assumption that, within a grid cell, the three velocities' components are only linear functions of their corresponding directions, i.e. u = u(x), v = v(y) and w = w(z).An alternative approach is to assume that u = u(x, y, z), v = v(x, y, z) and w = w(x, y, z) inside a grid cell, which might be more realistic in terms of representing unresolved motions.However, no such information is generally provided by GCMs.Furthermore, it would also require that the mass transports through the grid faces are unchanged in order to satisfy the continuity equation of the GCM.
The trajectory schemes integrate the trajectories from the volume or mass transports through the grid-box faces in contrast to many other trajectory schemes that only use the velocity fields.We will first describe how these fluxes are computed and then the three different trajectory schemes.

Mass and volume flux
The TRACMASS trajectory schemes are mass conserving as they, like the GCM, deal with the transport across the grid faces and the transport is only interpolated linearly between the two opposite faces in a grid box.The trajectories will hence never cross a grid boundary.
A GCM mesh is generally spherical or curvilinear.The longitudinal ( x i,j ) and the latitudinal ( y i,j ) grid lengths will hence be functions of their horizontal positions i, j on a curvilinear grid.The vertical coordinate in a GCM has a depth level thickness z n i,j,k , where k is vertical level and n is time step.Note that the vertical resolution can vary not only vertically but also both horizontally and in time, which makes it possible to use any vertical coordinate, e.g.sigma (Marsh and Megann, 2002), z-star, pressure or hybrid coordinates (Kjellsson and Döös, 2012b).The horizontal mass transports through the eastern and northern faces, respec-tively, of the i, j, k grid box at time step n are given by The zonal velocity u n i,j,k and the meridional velocity v n i,j,k are in the above equations on a C grid.It is, however, possible to use the velocities from A-and B-grid models, where the velocities are instead at the corners of the grid cell, leading to This averaging of two horizontal grid points in order to have the perpendicular velocity to the grid box in the middle on the grid face is exactly how a B-grid model discretises the equations when, e.g.solving the continuity equation.
Note that the mass transport can be replaced by the volume transport in models that assume the fluid to be incompressible, which is the case for most OGCMs.In other models (most AGCMs), we may use the hydrostatic approximation to write where g is gravity and p is air pressure.The mass transports through the lateral grid faces in the AGCM expressed by Eqs.
(1) and ( 2) will use Eq. ( 5) to determine z and hence become V n i,j,k = v n i,j,k x i,j p n i,j,k /g.
The vertical mass transport can similarly be computed from the vertical velocity w i,j,k through the upper face of the grid box so that The vertical velocity would in the equation above be taken directly from the stored velocity fields from the GCM.It is, however, in order to guarantee mass conservation, advantageous to instead calculate the vertical transport W n i,j,k from the continuity equation as the TRACMASS trajectory schemes rely on mass or volume continuity.The continuity equation, which expresses conservation of mass, states that Integrating Eq. ( 9) over a finite grid box of volume x y z, we obtain where M i,j,k is the mass of the grid box.The rate of mass change of the grid box ∂M i,j,k /∂t can, on the other hand, be due to (1) compression in an compressible GCM and/or (2) grid-box volume change, which generally in a GCM is due to the time dependence of the vertical resolution so that the thickness of model layers varies in time.
The mass of the grid box is where n is the time level of the stored GCM fields so that time is t = n t G and t G is the time interval between two stored GCM fields.
The vertical mass transport through the top of the grid box is obtained by discretising Eq. ( 10) between two stored time levels: which is computed by integration from the bottom and upwards with the bottom boundary condition W i,j,0 = 0.This is the same way the vertical velocity is computed in the GCM, except that we use the stored horizontal velocities and the grid-size thicknesses to ensure that they satisfy the time dependency correctly.In many OGCMs, the fluid is considered to be incompressible, and thus the density is constant and ρ can be dropped from all equations in order to have volume flux instead of mass flux in the calculations.The vertical volume transport through the top of the grid box becomes If, additionally, the vertical resolution is time independent, the last term can be neglected and thus On the other hand, in many AGCMs, there are both compressibility of the air and a time dependence of the vertical resolution, which is generally expressed in pressure and hence where Eq. ( 5) has been used.Note that in the case of offline calculations, one may instead use centred or forward finite time differences in Eqs. ( 12), ( 13) and (15).naud (1997) and used and developed for ocean mass transport studies by Döös (1995).The velocity inside a grid cell is found by assuming that it is only a function of its direction, i.e. u = u(x), v = v(y), w = w(z).Linear interpolation gives the zonal velocity and similarly for v(y) and w(z).Note that the calculation of the vertical mass transport W n i,j,k through the top face of a grid box, with Eqs. ( 12)-( 15), only involves the velocities on the considered grid box.A 3-D dependency of the velocities (u = u(x, y, z), v = v(x, y, z) and w = w(x, y, z)) would require velocities from other grid boxes, which could potentially break the mass conservation of Eqs. ( 12)-(15).
To calculate the zonal position, x, of a trajectory, we use u = dx/dt and write Eq. ( 16) as the differential equation − u i−1,j,k = 0.The drawback by solving the above differential equation is that x is not constant, and a horizontal grid face is rarely rectangular in a GCM.The solution will hence depend on the position of the trajectory in each grid box.Döös (1995) used therefore a x corresponding to the average latitudinal position of the trajectory in each grid box, which was obtained by computing the trajectories several times in each grid box.Blanke and Raynaud (1997) made this unnecessary by non-dimensionalising the position and used volume fluxes instead of velocities.By substituting x for a nondimensional position r ≡ x/ x i,j and t for a scaled time s ≡ t/( x i,j y i,j z i,j,k ), we get where F = dr/ds is the zonal volume or mass flux, and Its solution describes the zonal displacement within the grid box between the faces and is found using the initial condition r(s 0 ) = r 0 of its zonal position so that The scaled time s 1 becomes where r 1 = r(s 1 ) is given by either r i−1 or r i , when a trajectory enters the western or eastern grid face, respectively.
For a trajectory reaching the grid face r = r i or r = r i−1 , both F (r 1 ) and F (r 0 ) must be of the same sign in order for Eq. ( 19) to have a solution.If F (r 1 ) and F (r 0 ) are of opposite signs, there is a zero zonal transport at a position between r 1 and r 0 , and this position is reached exponentially slow.
The above procedure is repeated for meridional and vertical displacements, where now r = y/ y i,j or r = z/ z i,j,k .This yields non-dimensional position, r 1 , and scaled time, s 1 , for the zonal, meridional and vertical displacements of the trajectory, respectively, inside the grid box under consideration.The smallest transit time s 1 − s 0 and the corresponding r 1 denote through which grid face of the grid box the trajectory will exit and move into the adjacent one.The exact displacements in the other two directions are then computed using the smallest s 1 in the corresponding Eq. ( 18).
Note that Eqs. ( 18)-( 19) are not valid if the transport fields across the grid box are constant, i.e. when (F i−1,j,k = F i,j,k ), since it would imply a division by zero with β = 0 in both equations.The differential equation then simplifies to which has the solution and the scaled time s 1 is If F i−1,j,k = F i,j,k , TRACMASS instead uses Eqs. ( 21) and ( 22).

Stepwise-stationary and fixed GCM time step integrations
The trajectory scheme above is, strictly speaking, only valid for stationary fields.The scheme is, however, possible to use for time-dependent fields by assuming that the velocity and surface-elevation fields are stationary during a limited time interval.The stepwise-stationary method presented here consists of assuming that the fields are stationary during intermediate time steps between two GCM outputs and then updated successively as new fields become available.If this is undertaken online, i.e. in the same time as the GCM is integrated, this time interval will simply be the same as the time step the GCM is integrated by, which is typically between several minutes and a few hours in a global GCM.If instead the trajectories are calculated offline, the time intervals between GCM fields will be at least as often as the fields have been stored by the GCM, at intervals that can be days or even months.
A linear time interpolation of the velocity fields between two GCM velocity fields permits a simple way to have shorter time steps by which the fields are updated in time.The time interval between two GCM velocity fields is t G and the shorter time interval at which the fields are interpolated is t i as illustrated by Fig. 1.The number of intermediate time steps is hence the ratio I S = t G / t i .For any quantity in the GCM output, F , the value at intermediate time step m, located between GCM outputs n − 1 and n, is The coefficients β, δ in Eq. ( 17) are updated when a trajectory moves from one grid box to another.Thus, the time step for the trajectory, i.e. s 1 − s 0 , may be shorter than the intermediate time step, t. t i is hence the maximum possible time step for a given I S but is often shorter if the spatial grid spacing ( x, y, z) is small and t G long.We will therefore test TRACMASS by imposing constant velocities for the entire t G in order to mimic other codes, such as the Ariane code based on Blanke and Raynaud (1997), which do not make any temporal interpolations of the velocity fields.This particular case of the stepwise-stationary scheme with constant velocity fields for the entire period between two GCM outputs will be denoted the fixed GCM time step.These two schemes together with a truly time-dependent scheme, described in next section, will be tested.

Analytical time integration with the time-dependent scheme
The stepwise-stationary integration method presented in the previous section assumes that the velocity and the grid-box thicknesses remain constant throughout the time step, and only spatial variations of velocity are accounted for.Another approach is to interpolate the velocity fields not only in space within the grid box but also in time between the GCM outputs.This approach, introduced in TRACMASS by de Vries and Döös (2001), is more accurate but involves a more advanced differential equation to be solved and integrated along the trajectories.Accounting for both spatial and temporal variations of velocities in the trajectory scheme renders intermediate time steps unnecessary.We will later show that using a large number of intermediate steps, the stepwise-stationary scheme approaches this time-dependent scheme asymptotically.
The time-dependent scheme can be derived in the same way as Eq. ( 17), but instead of a linear interpolation in space, we use a bilinear interpolation in both space and time.As before, we use non-dimensional position r = x/ x and scaled time s ≡ t/( x y z), where the denominator is the volume of the particular grid box.For a zonal volume or mass flux F , a bilinear interpolation in space and time yields s is the scaled time step between two data sets: where t G is the time step between two data sets in true time dimension (seconds).The vertical grid-box spacing is for models with time-dependent grid cell thicknesses replaced with an average between the two time steps: z n + z n−1 /2.Similar expressions for the meridional and vertical directions can be derived.
Connecting the local transport to the time derivative of the position with F = dr/ds, the following differential equation is obtained: where the coefficients are defined by Different analytical solutions exist for the three cases: α > 0, α < 0 and α = 0, which together cover all possible values of α.The acceleration, inside the r − s grid box, is d 2 r/ds 2 = −αr − γ , which is constrained by a linear r-dependent term proportional to α and the constant γ .

√
2α and get where Dawson's integral has been used, as well as the initial condition r(s 0 ) = r 0 .An example of trajectories in this case is illustrated in Fig. 2a, with given values of We see here that α > 0 occurs when the flow changes from divergence in the i direction at time step n − 1 to convergence at time step n.   31) can be re-expressed as where the error function erf dx.An example of trajectories for this case is illustrated in Fig. 2b.We see here that α < 0 occurs when the flow changes from convergence in the i direction at time step n − 1 to divergence at time step n.

The case α = 0
The solution of Eq. ( 26) when α = 0 is This case would normally not occur in a realistic GCM integration since it would correspond to a field constant in time or space, where Note that if the fields are in steady state, Eq. ( 34) is reduced to become identical to the stationary solution of Eq. ( 18).An example of trajectories in this stationary case is illustrated in Fig. 2c.
If instead α = 0 since the fields are constant in space, i.e. the transport across the grid cell is constant (F i = F i−1 ), then we also have β = 0, which leads to a simplification of Eq. ( 26): with the solution An example of trajectories in this case with constant fields in space is illustrated in Fig. 2d.

The transit time
A major difference between the time-dependent and the stepwise-stationary schemes is that in the time-dependent scheme, the transit times s 1 −s 0 cannot in general be obtained explicitly with the time-dependent scheme in contrast to the stepwise-stationary analytical solution of Eq. ( 18).Using the solutions given by Eqs. ( 31)-( 34), the relevant root s 1 of r(s 1 ) − r 1 = 0 (37) has to be computed numerically for each direction.We will now describe how the roots s 1 and the corresponding exiting grid face r 1 can be determined.The displacement of the trajectory inside the grid box under consideration then proceeds as previously discussed for stationary velocity fields.
We now determine the roots s 1 of Eq. ( 37) and the corresponding r 1 needed to calculate trajectories inside a grid box.In what follows, s n−1 s 0 < s n and the relevant roots s 1 are to be in the interval of s 0 < s 1 s n .We also focus on the cases α > 0 and α < 0, since the forthcoming considerations can easily be adapted for the case of α = 0.For numerical purposes, we use As above, s is the scaled time.The coefficient in Eq. ( 38) appearing in Eqs. ( 31) and ( 33) is exactly zero when either the r i−1 or r i grid face represents a solid boundary, so that transport F i or F i−1 is zero for all n, respectively.In these instances, the opposite grid face fixes r 1 , and the root s 1 > s 0 can be computed analytically.If there is no solution, we take s 1 = s n .When all three transit times equal s n , the trajectory will not move into an adjacent grid box but will remain inside the original one.Its new position is subsequently determined, and the next time interval is considered. (a) Examples of how trajectories calculated with the time-dependent scheme evolve as a function of the transport F in the space interval r i−1 < r < r i and in the time interval s n−1 < s < s n , which hence corresponds to an interval between two GCM outputs ( t G ) and of a grid box ( x, y or z).The colour shows the transport values F obtained by the bilinear interpolation between the four corners (F n−1 i−1 , F n−1 i , F n i and F n i−1 ).(a) α > 0 with two corners of transport in the negative direction (F < 0), which correspond to westward, southward or downward directions, and one corner flowing in the opposite direction.(b) α < 0. (c) α = 0 and γ = 0 correspond to the stationary fields, which results in an F field that only changes in the (r) direction.(d) α = 0 and β = 0 correspond to the constant fields in space but which vary in time.Note that the F = 0 line between the red and blue colours corresponds to static flow, which results in "vertical" trajectories in the figures.
The roots of Eq. ( 37) have to be computed numerically if (βγ −αδ)/α = 0.This is also true for locating the extrema of the solutions given by Eqs. ( 31) and ( 33).Alternatively, one can consider the case F (r, s) = 0 using Eq. ( 24) to analyse where possible extrema are located.It follows that in the s−r plane, the extrema lie on a hyperbola of the form r = (as + b)/(c + ds).Obviously, only the parts defined by s n−1 ≤ s ≤ s n and r i−1 ≤ r ≤ r i are relevant.Depending on which parts of the hyperbola, if any, lie in this "box" and satisfy the initial condition r(s 0 ) = r 0 , the trajectory r(s) exhibits no, one or, at most, two extrema.In the latter case, the trajectory will not cross either the grid face at r i−1 or the one at r i (see Fig. 2 for an example).Hence, the trajectories r(s) determining the transit time s 1 −s 0 will have at most one extremum; i.e. there is at most one change of sign in the local transport F .
An efficient way of proceeding is as follows: first, consider the grid face at r i .For a trajectory to reach this grid face, the local transport must be nonnegative, which depends on the signs of the transport F n i−1 and F n i .Four distinct configurations may arise between the model outputs (s n−1 < s < s n ), where the calculation of the trajectory takes place: Figure 3.The four different cases of how trajectories might reach the grid face at r = r i .Note that the trajectories for case 4 can not reach r = r i .The background colours are the same as in Fig. 2 with F > 0 in red and F < 0 in blue.The dashed trajectories outside the grid box denote the necessary computed fictive paths for estimating when s = s 1 and if the trajectories reach r 1 (s 1 ) = r i .grid face at an earlier time, as illustrated by trajectory C in Fig. 3.This is only possible if case 3 applies for the grid face at r i−1 and s # > s 0 , in which case it is determined whether r(s # ) ≤ r i−1 .If this is not the case, there is a solution to r(s 1 )−r 1 = 0 for r 1 = r i and s 0 < s 1 ≤ s n .Subsequently, this root can be calculated numerically using a root-solving algorithm (Press et al., 2007).But if r(s n ) < r i or, if applicable, r(s # ) ≤ r i−1 , we proceed by considering the other grid faces.The arguments for the grid face at r i−1 are similar to those relating to r i .
If case 2 applies and s 0 < s * , we add here to the considerations given in case 1 using s * instead of s n .If there is a root for r 1 = r i , then s 0 < s 1 ≤ s * .This root is illustrated by trajectory D in Fig. 3 with (r 1 , s 1 ) = (r i , s 1D ).
For case 3, we follow the procedure given by case 1.If there is a root for r 1 = r i , then s # < s 1 ≤ s n .This root is illustrated by trajectory E in Fig. 3 with (r 1 , s 1 ) = (r i , s 1E ).
For case 4, no solution of Eq. ( 37) is possible for r 1 = r i , since all trajectories exit through the grid face located at r i−1 , as illustrated by trajectory G in Fig. 3, or will not reach any grid face during the time interval s n−1 < s < s n .We must then instead search for an exit through another of the six grid faces.
All these considerations are applied to each of the three spatial directions in order to determine through which of the six grid faces the trajectory will exit and at which position on the corresponding grid face.
Since the trajectories are unique solutions to Eq. ( 26) and the continuity equation is respected, the TRACMASS trajectories will therefore never hit any solid boundary, such as the coast or the sea floor, unless the sedimentation option is activated, where an extra velocity is imposed, a feature that was introduced in TRACMASS by Corell and Döös (2013).
An example of the evolution of trajectories calculated with the three different schemes within a time-space cell for α > 0 is shown in Fig. 4. The trajectories computed with the stepwise-stationary scheme approach the trajectory computed with the time-dependent scheme with increasing number of intermediate time steps (I S ).The fixed GCM time step trajectory can, however, not follow the time-dependent one since it does not update the velocities between the GCM outputs and consequently deviates immediately as it leaves the initial point (r 0 , s n−1 ).

Tests with different velocity fields
The results obtained from the stepwise-stationary scheme are now compared with those from the time-dependent trajectory schemes using two different sets of velocity fields.The first uses a high-resolution OGCM with z-star coordinates.The second uses a global atmospheric general circulation model Example of how the trajectories differ when computed with the fixed GCM time step method in orange, the stepwisestationary method in blue, purple and green as well as the timedependent method in red.They all start at the same time s n−1 and in the same position r 0 but exit the grid at different locations and times.Note that the stepwise-stationary method needs at least four intermediate time steps (I S ) to be close to time-dependent trajectory.The background colours are the same as in Fig. 2 with F > 0 in red and F < 0 in blue.
with hybrid pressure coordinates.For the stepwise-stationary scheme, five different settings of I S , i.e. the number of intermediate steps, are tested.The fixed GCM time step is also tested for comparison, although it is not a standard feature of TRACMASS.

Ocean trajectories with a high-resolution OGCM
Oceanic velocity fields for this case were obtained from a simulation with version 3.6 of the NEMO ocean model (Madec, 2016) in a global ORCA12 configuration.The horizontal resolution of the ORCA12 grid is approximately 1/12 • , corresponding to x ≈ 6 km at 50 • latitude.Model fields were available as 5-day averages every 5 days.The configuration uses 75 z * vertical levels with partial bottom cells, where z ranges from ∼1 m at the surface to 250 m in the deepest parts of the ocean.The z * coordinate approach permits large-amplitude free-surface variations relative to the vertical resolution (Adcroft and Campin, 2004).In the z * formulation, the variation of the column thickness due to seasurface undulations is not concentrated to the surface level, as in the z-coordinate formulation, but is equally distributed over the full water column.Thus, the vertical levels naturally follow the sea-surface variations, which also implies that they are time dependent, and we therefore have used Eq. ( 12) to calculate the vertical transport in TRACMASS with a time-dependent z n in the equation.The model was forced with 6-hourly atmospheric fields from what is known as the Drakkar forcing set version 4 (DFS4) (Brodeau et al., 2010).Subgrid processes were represented using 125 m 2 s −1 Laplacian isoneutral tracer diffusion and −1.25 × 10 10 m 4 s −1 bi-Laplacian viscosity.TRACMASS has been applied to this specific model integration already by Nilsson et al. (2013), where it was compared with surface drifters in the Agulhas region.This is also the region where we are going to test TRACMASS because of its complex time-dependent dynamics with travelling eddies, known as "Agulhas rings", which "leak" Indian Ocean water into the Atlantic Ocean as part of the conveyor belt.A total of 2193 trajectories were started, evenly spread over four horizontal grid boxes at all depths in the Indian Ocean and followed for 50 days as shown in Figs. 5 and 6.

Atmospheric trajectories with an AGCM
In order to test the trajectory schemes in the atmosphere, we have used the ERA-Interim reanalysis (Dee et al., 2011) from the European Centre for Medium-Range Weather Forecasts (ECMWF) simulated with the IFS (Integrated Forecasting System) model.In this ERA-Interim data set, the vertical coordinate is a terrain-following hybrid coordinate (Simmons and Burridge, 1981), where the pressure at the lower interface of level k is given by p k = A k + B k p s , where p s is the surface pressure and A k and B k are parameters at the level k ∈ [0, 60], with p 60 = p s and p 0 = 0.As in the NEMO ocean model, the grid cell thickness varies in time, and we calculate vertical mass flux from the continuity equation (Eq.15).The ERA-Interim data used here had a horizontal resolution of 1.25 • and is available 6-hourly ( t G ). Trajectories are shown in Fig. 7.They were initiated every 6 h from a grid cell air column over the Eyjafjallajökull volcano eruption during 14-18 March 2010.The trajectories were evenly distributed horizontally and started in exactly same positions for the tests with different time steps using the stepwise-stationary scheme and time-dependent case.

Lagrangian statistics
The average distance between the trajectories obtained with the time-dependent scheme and the five different stepwisestationary cases as well as the fixed GCM time step case are shown in Fig. 8.The distances from the time-dependent trajectories after 50 days for the OGCM case and after 10 days for the AGCM case are presented in Table 1.These average distances have been possible to compute since of all the individual trajectories were started in the exact same positions for the different cases.Results clearly show that the distance between trajectories calculated with the stepwise-stationary scheme and those calculated with the time-dependent scheme decreased as the number of intermediate time steps was increased.The fixed GCM time step case, i.e. when no inter-Figure 5. Agulhas trajectories started evenly distributed in a square of four grid cells and followed for 50 days.Colouring is used to separate the trajectories from each other.
Table 1.The table shows the average distance between the time-dependent integrated trajectories and the stepwise-stationary integrated ones at the end of simulations, which is 50 days for the OGCM and 10 days for the AGCM.I S is the number of intermediate time steps between two GCM outputs.The "maximum time step" stands for the intermediate time step lengths ( t i ), which are used in the different trajectory integrations.The last column is the computational time normalised with regard to the time-dependent case, where theoretical velocity fields are used to compute trajectories, i.e. with no data reading or writing.mediate time steps are used, shows the greatest distance to the time-dependent case.Standard Lagrangian statistics have also been computed for the ocean trajectories (Fig. 9), with the definitions given in the Appendix.The "relative" and "absolute dispersion" as well as the "mean displacement" of the trajectory cluster show how the cluster will disperse and move in time.They reveal a similar pattern, where only the fixed GCM time step case differs from the others.The fixed GCM time step differs already after 3 to 4 days, which should be related to the fact that the GCM velocities are updated every 5 days (= t G ) in this OGCM case.
The Lagrangian velocity autocorrelation, which describes the correlation of the velocity of the trajectories at one time with that of previous times, shows in Fig. 9 how all cases except the fixed GCM time step give nearly the exact same correlation.The Lagrangian timescale, which is computed from the autocorrelation and is a measure of the The timedependent method results are in red as well as those obtained with the stepwise-stationary method with I S = 1, 12, 120, 1200 and 12 000 as well as fixed GCM time steps.Note that these homologous trajectories were selected to illustrate that stepwise-stationary trajectories are closer to time-dependent trajectories when the number of intermediate time steps (I S ) is increased.memory of the trajectories, reflects the same feature with a Lagrangian timescale of approximately 3.9 days for the time step and the time-dependent cases but a slightly shorter timescale of 3.4 days for the fixed GCM time step case.The Lagrangian timescale based on observations with surface drifters is clearly shorter than this both for the global ocean (Döös et al., 2011) and in the Agulhas region (Nilsson et al., 2013).This relatively shorter Lagrangian timescale (hence closer to observations for the fixed GCM time step) is simply due to the abrupt changes in the velocity fields every time these are updated.A realistic shortening of the Lagrangian timescale can only be obtained by incorporating finer scales by increasing the GCM resolution or adding subgrid parameterisations to the trajectories.
The power spectra computed from the Lagrangian velocities show that the fixed GCM time step was more energetic than the other schemes, which all yielded nearly identical results.This is the case for all frequencies.There is also a weak maximum at four cycles per day (6 h), which remains unexplained, although it may be related to the fact that the OGCM uses 6-hourly atmospheric forcing.

Lagrangian stream function and residence time
The mass conservation properties of the used trajectory schemes make it possible to calculate mass transports between different sections in the model domain (Döös, 1995).The approach is that one can associate each trajectory par-ticle with a mass or volume transport.This requires that enough trajectories are computed to fill the model grid in space and time with a sufficient number of trajectories.Lagrangian stream functions can be calculated by summing over trajectories representing a desired path (Blanke et al., 1999;Döös et al., 2008;Kjellsson and Döös, 2012b).The difference between the Lagrangian and the more common Eulerian stream functions is that with the Lagrangian one can isolate a particular path between a starting and an ending section in the ocean or the atmosphere.
The influence of the different trajectory schemes on the inter-ocean exchange of water masses, which takes place in the Agulhas region, has been evaluated by calculating Lagrangian stream functions.Figure 10 shows the Lagrangian barotropic stream function computed from trajectories using the time-dependent scheme and the fixed GCM time step scheme.Lagrangian decomposition has been used to compute two separate stream functions for each scheme, one from trajectories entering the Atlantic and one from those returning back into the Indian Ocean via the Agulhas retroflection region.The time-dependent scheme favours slightly (one additional stream line) the entering into the Atlantic compared to the fixed GCM time step scheme.This is also clearly visible when computing the residence time; i.e. the time trajectories stay within the Agulhas region as shown in the lower righthand panel of Fig. 9.The first particles start to exit the Agulhas region, as defined by the map in Fig. 10, after 50 days.The number of trajectory particles then decays exponentially with an e-folding time of about 210 days.This is rather similar for all trajectory-scheme integrations.There is, however a clear difference in the results where the trajectories exit.The fixed GCM time step scheme results in 38 % flowing into the Atlantic and 59 % into the Indian Ocean after 800 days.All the other trajectory integrations yield very similar results but with 46 % flowing into the Atlantic and 52 % into the Indian Ocean.This suggests that the fixed GCM time step scheme does not capture the same behaviour as the other schemes.
We have repeated the above ocean-trajectory experiment by releasing the particles in other time periods and increasing the ensemble size.The results only changed marginally.

Computational speed
In addition to the higher accuracy of the time-dependent scheme, it was also shown to be computationally faster than the stepwise-stationary scheme with intermediate time steps.In order to quantify this difference, we compared the computational time for the different schemes using analytical velocity fields describing inertia oscillations (Döös et al., 2013), where no data needed to be read nor written since the velocity fields have a known analytical solution and disk storage was switched off.These computational times, shown in the last column of Table 1, have been normalised by dividing with the time obtained with the time-dependent scheme.The stepwise-stationary scheme was only as computationally fast as the time-dependent scheme when no extra intermediate time steps were taken between two readings of the velocity fields (I S = 1) or when using fixed GCM time steps.When the number of intermediate time steps was increased to 12 000, the stepwise-stationary scheme was more than 1000 times slower.The total of 12 000 intermediate steps was also approximately the number of intermediate time steps required in order to obtain as accurate results as those obtained from the time-dependent scheme.

Discussion and Conclusions
The two trajectory schemes available in TRACMASS have here been intercompared by calculating Lagrangian statistics, transports and the distances between the trajectories.This has been done for both oceanic and atmospheric applications.The stepwise-stationary scheme assumed that the velocity fields were stationary for the duration of a userdefined intermediate time step between model output fields.These velocities are, however, updated with a linear interpolation in time when crossing a model grid face.The time-dependent scheme does not assume that the velocity is in steady state during any time interval since it solves the differential equations of the trajectory path not only in space but also in time.This continuous evolution of the time-dependent scheme makes it more accurate than the stepwise-stationary scheme without any significant increase in computational expense.
In addition to these two TRACMASS schemes, we have tested a fixed GCM time step scheme, which is in fact a special case of the stepwise-stationary scheme but with velocity fields always remaining in steady state until a new GCM data set is reloaded in order to mimic the Ariane trajectory model (Blanke and Raynaud, 1997).A consequence of only updating the fields at the GCM output times is that the velocities are assumed to be in steady state for long periods and then changed abruptly with a discontinuity.
The accuracy of the schemes has been evaluated by comparing the distance between particles that have been started from the same positions but with different trajectory schemes and how this distance evolves in time.This distance was shown to depend on the scheme and the number of intermediate time steps for the stepwise-stationary case.The aver- age distance as a function of time between the trajectories obtained from the different schemes is shown in Fig. 8, and their end position distances are shown in Table 1.
The study has shown that the TRACMASS timedependent scheme is likely to be more accurate as well as faster than the scheme with intermediate steps.It remains to be shown how the trajectory schemes used in the present study compare to other trajectory schemes, e.g.Runge-Kutta, which could be used where mass conservation is not important.
The stepwise-stationary scheme needed up to 12 000 intermediate time steps to give as accurate trajectory paths as the time-dependent scheme, which is more than a thousand times as computationally expensive when reading and writ-ing are excluded.The distance between trajectories calculated with the time-dependent scheme and those obtained with the stepwise-stationary scheme decreased as the number of intermediate time steps is increased.The greatest distance was obtained when no temporal variations between GCM outputs at all were considered, i.e. with the fixed GCM time step scheme.We thus conclude that the timedependent scheme is the most accurate of those tested here for two reasons.Firstly, for theoretical reasons since the time-dependent scheme does not assume stationary velocities during any period of time.Secondly, the trajectories computed with the stepwise-stationary scheme converge towards those computed with the time-dependent scheme for an increasing number of intermediate time steps.A future study www.geosci-model-dev.net/10/1733/2017/could be to calculate trajectories first using fields stored at each GCM time step and then using fields stored at longer time intervals.In the first case, trajectories would be very accurate and could represent a "truth", and the second case could be used to evaluate which scheme is the closest to the truth.
The Lagrangian statistics, such as relative and absolute dispersion as well as Lagrangian velocity autocorrelation functions and power spectra, showed almost identical results for the time-dependent and the stepwise-stationary schemes.The fixed GCM time step showed, however, some differences from the other two schemes.For example, the dispersion after 3-4 days was slightly larger for using a fixed GCM time step, which might be explained by an abrupt change every time the GCM velocities are updated compared to the smoother transition of the two other schemes.The results show that the fixed GCM time step method does not capture the same behaviour of trajectories as the other schemes.The Lagrangian statistics are also clearly affected by the model resolution and the time sampling of the GCM fields (Döös et al., 2011;Kjellsson and Döös, 2012a;Kjellsson et al., 2013;Nilsson et al., 2013).Future improvements to the TRACMASS model will involve improvements of the subgrid turbulence parameterisations, which could give more realistic dispersion properties.
The mass conservation of the trajectory schemes in the present study arises from that (1) mass transports across the grid faces are used in the same way as in the GCM itself instead of velocities as in most other trajectory schemes; (2) the mass transport is linearly interpolated within the grid box, where there is otherwise no information of the velocity from the GCM and this enables us to set up a differential equation which has an analytical solution of the trajectory within the grid box.The different trajectory schemes, although mass conserving, will not yield the same results in terms of transports between different sections.The mass transport was tested in the Agulhas experiment, where the fixed GCM time step scheme relatively favoured the Agulhas retroflection with more trajectories returning into the Indian compared to the time-dependent and stepwise-stationary schemes.This difference in mass transport can be explained by the delicate path of the Agulhas leakage, which requires an accurate temporal evolution so that particles can be retained in Agulhas rings.This was better achieved by the time-dependent and stepwise-stationary schemes than by the fixed GCM time step scheme.
The TRACMASS trajectory code with corresponding schemes has been improved and has become more accurate and user friendly over the years.An outcome of the present study is that we strongly recommend the use of the timedependent scheme based on de Vries and Döös (2001) in favour of the stepwise-stationary scheme.We would also like to dissuade the use of the more primitive fixed GCM time step scheme, which is used in other trajectory codes, since the velocity fields remain stationary for longer periods, creating abrupt discontinuities in the velocity fields and yielding inaccurate solutions.We have only tested one OGCM and one AGCM simulation here, but we speculate that at coarser resolution in both space and time, the differences obtained with the two schemes would increase.However, in non-eddying simulations (e.g. 1 • ocean models) this may not be true due to the low variability of the flow.
The TRACMASS strict requirement of mass conservation makes it, however, necessary to have complete velocity fields on the original GCM grid in order to use mass or volume transports in and out of each model grid box.This requirement of mass conservation will always be somewhat more demanding than for other trajectory codes, since it requires a total understanding of the various GCM coordinate systems as well as the incorporation of them in the TRACMASS framework.This state of affairs is in marked contrast to what holds true for various trajectory codes that only require velocity fields with no mass conservation.

Figure 1 .
Figure 1.Schematic illustration of how the transport fields F (t) are updated and interpolated in time between the stored GCM data, which are read in at the time t n and are separated in time by the time interval t G (in red).The fields are then linearly interpolated at the points in blue with intermediate time steps.The number of intermediate time steps between two GCM velocities is in this example I S = t G / t i = 4.

Figure 6 .
Figure6.Example of ocean trajectory paths due to different trajectory schemes and number of intermediate time steps.The timedependent method results are in red as well as those obtained with the stepwise-stationary method with I S = 1, 12, 120, 1200 and 12 000 as well as fixed GCM time steps.Note that these homologous trajectories were selected to illustrate that stepwise-stationary trajectories are closer to time-dependent trajectories when the number of intermediate time steps (I S ) is increased.

Figure 7 .Figure 8 .Figure 8 .
Figure 7. Example of atmospheric trajectory paths starting form the Eyjafjallajökull volcano during its eruption calculated with different trajectory schemes and number of intermediate time steps.The same colour coding of the trajectories as in Fig. 6 is used.Note that the red time-dependent and the blue stepwise schemes with I S = 12 000 trajectories are nearly identical.

Figure 9 .Figure 9 .
Figure 9. Lagrangian statistics of the ocean Agulhas trajectories.The relative dispersion (top left), the absolute dispersion (top right), the mean displacement travelled by the trajectory cluster (middle left), the average Lagrangian velocity autocorrelation of the trajectories (middle right).The average power spectra of the Lagrangian velocities (lower left).The residence time evolution of the trajectory particles in the Agulhas region.Note that all statistics show very similar results, where only those based on the "fixed GCM time step" (orange curves) differ from the rest.

Figure 10 .
Figure 10.The Lagrangian decomposed barotropic stream function based on the particles released as previously but followed until they left the Agulhas region into the Atlantic (a, c) or the Indian Ocean (b, d).The top panels show the time-dependent scheme and the lower panels show the fixed GCM time step scheme.Note that there is more water (one stream line extra) flowing into the Atlantic with the time-dependent scheme than with the fixed GCM time step scheme, which instead favours relatively the flow into the Indian Ocean.Stream line intervals of 8 Sv (10 6 m 3 s −1 ) are shown.