The Lagrangian Analysis Tool (LAGRANTO) is adopted and applied to ECMWF’s latest ocean reanalysis. The primary motivation behind this study is to introduce and document LAGRANTO Ocean (LAGRANTO.ocean) and explore its capabilities in combination with an eddy-permitting ocean reanalysis. The tool allows for flexibly defining starting points, within circles, cylinders, or any user-defined region or volume. LAGRANTO.ocean also offers a sophisticated way to refine a set of computed trajectories according to a wide range of mathematical operations that can be combined into a single refinement criterion. Tools for calculating—for example, along-trajectory cross sections or trajectory densities—are further provided. After introducing the tool, three case studies are presented, which were chosen to reflect a selection of phenomena on different spatial and temporal scales. The case studies also serve as hands-on examples. For the first case study, at the mesoscale, ocean trajectories are computed during the formation of a Gulf Stream cold-core ring to study vertical motion in the developing eddy. In the second example, source waters are traced to the East Greenland Spill Jet. This example highlights the usefulness of a Lagrangian method for identifying sources or sinks of buoyancy. The third example, on annual time scales, focuses on the temporal evolution of extreme potential temperature anomalies in the South Pacific and the memory of the involved water parcels. Near-surface trajectories reveal that it takes approximately 5 months after the highest temperature anomaly before the involved water parcels cool to their climatological mean values at their new positions. LAGRANTO.ocean will be made publicly available.
The Lagrangian analysis of fluid motion is the natural framework for studying the temporal evolution of fluid parcels in time and space and the associated evolution of pertinent variables—for example, salinity or potential density—along the flow. The study of quasi-Lagrangian (i.e., constraint to the horizontal) ocean floaters or drifters has a long and rich history—see Davis (1991) and Rossby (2007) for comprehensive reviews—that dates back at least to the float experiments conducted by Swallow (1955). The Argo profiling floats are among the most prominent examples of quasi-Lagrangian observations of the state of the ocean (Roemmich et al. 2009; Riser et al. 2016). One drawback of drifting floats is that they are advected with the horizontal flow and hence the float trajectories do not follow the vertical motion of parcels along sloped isopycnals or vigorous vertical motion during convection. For example, Riser and Rossby (1983) reported that the horizontal motion of quasi-Lagrangian floats—in their case, the SOFAR float—can substantially differ from that of a three-dimensional trajectory. Three-dimensional trajectories can provide additional insights into the nature and the dynamics of the flow, particularly if combined with observation-constraint data, such as gridded reanalysis data.
Fully three-dimensional Lagrangian trajectories are difficult to achieve in the real ocean, which is one reason why Lagrangian analyses are commonly applied to ocean models, and several Lagrangian tools exist in the ocean modeling community and are worth acknowledging. Among such tools are ARIANE (http://stockage.univ-brest.fr/~grima/Ariane/; Blanke and Raynaud 1997; Blanke et al. 1999, 2001) and TRACMASS (tracmass.org; Döös 1995; de Vries and Döös 2001; Döös et al. 2011, 2013). Below, we present a brief background survey, which is meant to provide an impression of the diversity of previous studies based on computed parcel trajectories. This overview is not intended to provide a complete overview of the existing literature.
One of the first studies based on water parcel trajectory modeling using 4D gridded data is the work by Döös (1995), who studied the flow and mass exchange in the Southern Ocean. TRACMASS emerged from this study. Later, Blanke et al. (2001) extended this work and investigated the interocean exchange of mass on a global scale. Water mass transformation has been studied from a Lagrangian perspective, for example, in the southeast Atlantic by Rimaud et al. (2012), in the Denmark Strait by Koszalka et al. (2013), and in the Solomon Sea by Melet et al. (2013). Lagrangian modeling allowed Drijfhout et al. (2003) to reveal the influence of ocean eddies on the structure of the thermohaline circulation. The dynamics underlying the Deacon cell, the wind-driven meridional overturning component of the Antarctic Circumpolar Current, was studied based on water parcel trajectories by Döös et al. (2008). On a local scale, Jönsson et al. (2011) used water parcel trajectories to identify the water mass sources in the Gulf of Finland. More recently, Thomas et al. (2015) applied a Lagrangian approach to decompose the Atlantic meridional overturning circulation (AMOC) according to the contributions from different subduction zones. Nutrient advection in the Arctic Ocean (Popova et al. 2013), biological production (Lachkar and Gruber 2013), jellyfish transport (Berline et al. 2013), and Agulhas leakage (Rühs et al. 2013) are further research questions that have been addressed using water parcel trajectories. Note that Lagrangian modeling is of importance in the planning of search and rescue missions at sea—see Breivik et al. (2013) for a review—and the modeling of oil spills (Jones et al. 2016).
Because of the chaotic nature of the ocean circulation, trajectories are highly sensitive to their starting position in time and space, which makes a single trajectory unique and statistical analyses of large trajectory ensembles a necessity (LaCasce 2008) to obtain robust results. Nevertheless, during a limited period, the flow can organize not only with coherent flow features in space but also with respect to temporal changes in physical properties along the flow. During these periods trajectory ensembles can be identified and analyzed, and the results are more robust compared to analyses based on single trajectories. The strategy that we adopt in this study is widely used in the atmospheric science community, and it involves the computation of a large ensemble of trajectories from which coherent flow features are identified with posterior selection criteria, for example, the rate of change of potential density within a period.
The tool used to compute ocean trajectories in this study is the Lagrangian Analysis Tool Ocean (LAGRANTO.ocean, available online at www.lagranto.ethz.ch). LAGRANTO.ocean is based on the core code of LAGRANTO (Wernli and Davies 1997a; Sprenger and Wernli 2015), which thus far has primarily been used to study atmospheric phenomena across a wide range of spatial and temporal scales and in different datasets, for example, ERA-40 (Uppala et al. 2005), ERA-Interim (Dee et al. 2011), and the Twentieth Century Reanalysis (20CR; Compo et al. 2011). Although it is primarily an offline tool—that is, it makes use of model output data—an online version of LAGRANTO that computes trajectories simultaneously while the simulation is performed has recently been implemented in the regional numerical weather prediction model COSMO (Miltenberger et al. 2013). For studies with a focus on very short time scales—that is, convection in the ocean or atmosphere—online trajectory modeling appears to be a promising future direction. The adaptations made to LAGRANTO required to handle ocean data, the bathymetry, and further numerical details of LAGRANTO.ocean are briefly described in section 2 of this study. A tutorial and more detailed documentation are provided online (at www.lagranto.ethz.ch), and specific numerical details and the organization of the code are documented by Sprenger and Wernli (2015).
The example case studies aim to highlight the potential benefits and complementary insights that can be obtained using (offline) water parcel trajectories. Specifically, we aim to shed light on the following questions, which are chosen to reflect a wide range of spatial and temporal scales to help us document LAGRANTO.ocean:
Case I: Lagrangian flow through a Gulf Stream cold-core ring. What does the three-dimensional mesoscale flow through a forming Gulf Stream cold-core ring look like, and are parcel trajectories based on Ocean Reanalysis System 5 (ORAS5) depicting the mesoscale flow? (section 3a).
Case II: Sources of dense water on the East Greenland shelf. What is the temporal evolution of overflow water that forms on the East Greenland shelf, and does it have the potential to contribute to the East Greenland Spill Jet? (section 3b).
Case III: Persistence of seasonal sea surface temperature anomalies. How long is the memory for a water parcel that was anomalously warmed in the South Pacific? (section 3c).
We argue, similar to Foukal and Lozier (2016), that water parcel trajectories are particularly well suited for illuminating the underlying mechanisms of, for example, sources and sinks of salinity or temperature anomalies and their propagation through the oceans, and that a truly Lagrangian approach is a useful complementary approach to lead–lag correlations as used in earlier studies (e.g., Sutton and Allen 1997; Eldevik et al. 2009).
This study has a descriptive nature and is designed to document the tool and several applications, serving as future reference for more detailed studies. The remainder of this paper is organized as follows. First, the underlying dataset is introduced in section 2. Next, the structure of LAGRANTO.ocean, numerical details, and the boundary conditions are documented (section 3). Three case studies with three different temporal and spatial scales are discussed in section 4: the first example in section 4a, the second example in section 4b, and the third example in section 4c. The case studies are also showcased to illustrate the handling of LAGRANTO.ocean. A summary and outlook are provided in section 5.
The database for this study is the latest eddy-permitting ocean reanalysis product ORAS5 of the European Centre for Medium-Range Weather Forecasts (ECMWF). This reanalysis has been developed based on Ocean Reanalysis Pilot 5 (ORAP5) (Zuo et al. 2017; Tietsche et al. 2017), using the same ocean and sea ice model, method, and resolution. The NEMO ocean model (Madec 2008) is used to compute the eddy-permitting reanalysis. The grid spacing is approximately 0.25° outside the tropics. In the tropics, the grid is magnified to approximately 3 times this resolution to better resolve the narrow equatorial current systems. The output is available as 5-day averages. The standard model has been further modified to account for three wave effects (Breivik et al. 2015). These effects include a momentum flux estimated from the dissipation term (accounting for the intensity of breaking waves) of the ECMWF Wave Model (EC-WAM) component of the atmospheric reanalysis (ECMWF 2013). The surface boundary condition of the turbulent kinetic energy equation is likewise modified to account for the energy flux from breaking waves (Craig and Banner 1994). Finally, the Coriolis–Stokes forcing term is introduced in the momentum equation (Breivik et al. 2014, 2016).
The assimilation system uses more up-to-date observational datasets than ORAP5. ORAS5 uses the recently released quality-controlled EN4 in situ dataset (Good et al. 2013) with better vertical resolution and extended coverage than the previous version EN3 used in ORAP5. The altimeter sea level data have also been updated to use the latest version [Data Unification and Altimeter Combination System 2014 (DUACS2014)] from Archiving, Validation, and Interpretation of Satellite Oceanographic Data (AVISO). The underlying sea surface temperature (SST) analysis before 2008 has also been changed, and in ORAS5 it is from the HadISST2 dataset—the same used in the currently produced atmospheric reanalysis ERA5. The sea ice concentration product is the same as in ORAP5.
Concerning the quality of the reanalysis data, a weak spurious trend has been found in the wind field due to changes to the assimilation data. Furthermore, there is a trend in the wave field related to the introduction of assimilated wave height from altimeters after 1992 (Aarnes et al. 2015). However, these effects are considered to be too weak to substantially affect the ocean reanalysis. For the purpose of this study, the performance of the reanalysis is of secondary importance because our method can be applied to any type of gridded 4D data for which 3D currents are available. Because of its comparably high spatial and temporal resolutions, ORAS5 is particularly suited for the purposes of this study.
3. The structure of LAGRANTO.ocean
The LAGRANTO.ocean code is a modular Fortran code that includes input/output (I/O) routines for handling netCDF and ASCII data. The atmospheric components are made publicly available via an ETH Zürich–hosted web server (www.lagranto.ethz.ch), where a detailed user manual and case study examples are provided. The structure of the code is documented in Sprenger and Wernli (2015). We plan to distribute LAGRANTO.ocean also via this web page. One of the strengths of LAGRANTO is its flexibility in defining starting positions or starting volumes (e.g., circles, boxes, lines, or user-supplied polygons on model levels or any type of isosurface), computing back and forward trajectories and tracing 2D and 3D fields along the flow not only at the current parcel position but also in its vicinity (e.g., potential density can be followed along the flow but also, for example, 100 km to the east of a trajectory). Furthermore, LAGRANTO.ocean allows for computing vertical cross sections along the trajectories or outputting the results at every time step as a parcel density (i.e., number of trajectories per area) as either ASCII or netCDF data. Such functionalities to refine the computed trajectories make LAGRANTO.ocean a powerful tool for studying ocean dynamics from a Lagrangian perspective. Future developments and improvements made to the atmospheric components of the code will be incorporated into LAGRANTO.ocean.
a. Numerical details
Water parcel trajectory calculations are performed by building on the core code of the Lagrangian analysis model LAGRANTO 2.0 (Wernli and Davies 1997a; Sprenger and Wernli 2015), which has been widely used in its original version by the atmospheric community for the past 20 yr and has been applied to ECMWF’s ERA-40, ERA-Interim, 20CR, and to regional numerical weather prediction model output (e.g., WRF and COSMO) (e.g., Wernli and Davies 1997a,b; Sprenger and Wernli 2003; Knippertz and Wernli 2010; Martius et al. 2013; Schemm et al. 2013; Miltenberger et al. 2013; Schemm and Wernli 2014; Schäfler et al. 2014; Schäfler and Harnisch 2015; Schemm et al. 2016). The method solves the trajectory equation
where u denotes the three-dimensional current velocity at the parcel position x = (λ, ϕ, z). The parcel is assumed to be a fluid particle of negligible size (Provenzale 1999; Pasquero et al. 2007). The trajectory location is obtained through an iterative procedure. This procedure starts with a forward Euler time step using the velocity at the parcel position and is followed by three iteration steps using an adjusted mean velocity in Eq. (1). The adjusted mean velocity is obtained from averaging between the velocity at the current parcel position and the interpolated next position [for details, see, e.g., Figs. 2, 3 in Sprenger and Wernli (2015)]. For higher accuracy, the user can manually increase the number of iterations. A fourth-order Runge–Kutta scheme to solve Eq. (1) is also available. Bilinear interpolation is used in the horizontal, and linear interpolation is used in the vertical. For the linear interpolation in time, a 6-h time step is used, which is 1/20 of the input data’s time resolution and is adequate for complying with the Courant–Friedrichs–Lewy criterion for numerical stability (Seibert 1993). The 5-day-averaged time resolution itself is only marginally adequate for rapidly overturning deep convection (if the underlying model is capable of resolving convection), but it is sufficient to characterize the large-scale flow. Trajectories based on time-mean flow are sufficiently accurate to conserve mass (Nehrkorn et al. 2010). Future developments will include an additional eddy diffusivity term in the trajectory equation.
b. Code structure, documentation, and handling
LAGRANTO is controlled via the command line, and it consists of four core components and several post- and preprocessing tools. Additionally, we provide routines for plotting the output for NCAR Command Language (NCL), MATLAB, and Python. Each of the four core components serves a single purpose and is controlled individually from the command line. Every component comes with its own manual page (man page), which provides additional examples and general guidance. The following overview parallels parts of the LAGRANTO core documentation (Sprenger and Wernli 2015) but with an additional focus on the ocean-specific aspects. The four core components and their purposes are as follows:
startf: defines the trajectory start positions.
caltra: computes trajectories.
select: refines a set of trajectories.
trace: traces additional fields along a set of trajectories.
If the full LAGRANTO software suite is installed, then a call to compute water parcel trajectories based on ORAS5 data is caltra.ocean to help distinguish between input data from different sources. For example, based on COSMO output data, it is caltra.cosmo, and for ERA-Interim, it is caltra.ecmwf (similar for select, trace, and startf). Next, the individual components are introduced in greater detail.
The general command structure is
startf.ocean date outputfile 'horizontal specification @ vertical specification @ unit’
Date (format: YYYYMMDD_HH) and “outputfile” is the name of the output file where the starting points are stored as ASCII formatted text. To define the horizontal and vertical specifications of the starting region, LAGRANTO.ocean offers several geometric shapes, for example, circle, box, point, and line or user-supplied polygons. In the vertical, the user can specify a single level or a number of levels between two layers.
For example, if we wish to release trajectories at 0000 UTC 12 January 2010 from every grid point inside a box in the central North Atlantic encompassing 45°–50°N, 25°–30°W in the horizontal and in the vertical between a depth of 100 and 500 m in steps of 100 m, the corresponding one-line command is
startf.ocean 20100112_00 start.txt . .. . 'box.grid(25,30,45,50)@profile(100,500,5)@m',
where box.grid takes minima and maxima latitudes and longitudes as input parameters, and profile takes the top and bottom of the layer and the number of equidistant steps. Rather than selecting every grid point inside the box as a starting point, we may wish to choose equidistant steps of, for example, 100 km. In this case, we replace 'box.grid (minlon,maxlon,minlat,maxlat)' with 'box.eqd (minlon,maxlon,minlat,maxlat,ds)', where ds denotes the step size (km). Rather than box we may choose from circle, line, and point, or specify an input file with user-defined arbitrary coordinates. The latter is useful if trajectories are to start from a curved polyline. Furthermore, it is possible to specify vertical start points between layers of potential temperature or potential density, or on a single isosurface or model level. To start the computation from a single potential temperature level, we simply replace 'profile' with 'level(18)@K'. The man page of startf.ocean, which is called by the command lagrantohelp startf, provides a more complete list and additional examples of how to define starting regions.
The general command structure is
caltra.ocean startdate enddate inputfile outputfile options.
If the end date is before the start date, then backward trajectories are computed. The input file is the output from startf.ocean, and output is the file name of the text file where the trajectories are stored. Furthermore, options can be used to specify, for example, the output frequency in minutes. For example, if we wish to output the interpolated trajectory positions every 24 h we use “-o 1440”. Moreover, there are different options to handle the boundary conditions or the time step used in the trajectory computation. The options are documented on the man page, which can be viewed using lagrantohelp caltra. For example, if we wish to compute trajectories backward for 2 months from the previously defined box in the central North Atlantic and output their positions every 24 h, then the corresponding command reads as follows:
caltra.ocean 20100112_00 20091112_00 startpoints.txt traj.txt -o 1440.
The general command structure is
trace.ocean inputfile outputfile.
The trace component of LAGRANTO.ocean reads in the output from caltra.ocean and traces additional fields along the computed trajectories. By default, caltra traces only time, latitude, longitudes, and depth along the parcel position. The user can specify additional fields by creating a simple text file called tracevars in the working directory. This file contains the name of the field to be traced and the filename where the additional field is stored. For example, if we wish to trace potential density and the surface heat flux, the content of tracevars is
where the first column contains the variable names to be traced, the second column is a scaling parameter that can be used to transform the units of the variable, and the third column indicates whether the variable exists (0) or if LAGRANTO needs to compute it (1). The atmospheric component of LAGRANTO is able to compute several variables—for example, potential vorticity—even if the variable is not contained in the model output. In this case, the third column is 1. For LAGRANTO.ocean, this option is currently not fully implemented. The last column contains the prefix of the file in which the variable is stored. In this example, trace.ocean assumes that the variable sigmat is stored in PYYYYDDMM_HH and shflx is stored in SYYYYDDMM_HH. If a 2D variable is traced, then the trajectory’s depth information is ignored.
The general command structure is
select.ocean inputfile outputfile 'operator:variable:threshold:time'.
The select.ocean component of LAGRANTO.ocean offers the possibility to further refine the previously computed set of trajectories. For example, to select all trajectories that have potential density values greater than 27.8 kg m−3 after 72 h, the corresponding command reads as follows:
select.ocean traj.txt traj.filter.txt 'GT:sigmat:27.8:72'
Currently, the following set of operators can be applied: GT (greater than), LT (lower than), EQ (equal), IN (within a certain range), and OUT (outside a certain range), as well as TRUE and FALSE. The latter two operators are useful to check for a specific value. Additionally, it is possible to apply additional operators to a variable before the selection criterion is checked. In the abovementioned example, we may choose to check and select only those trajectories that, for example, reach a maximum in potential density greater than 27.8 kg m−3 during a certain period in hours after start using 'GT:sigmat(MAX):27.8:72,144'. Currently available as operators applicable to the variable of interest are ABS (absolute value), VAR (variance over a period), MIN (minimum over a period), MAX (maximum over a period), SUM (sum over a period), CHANGE (absolute value of the difference between two time steps), DIFF (difference between two time steps), and MEAN (mean over a certain period). Several filter criteria can be combined using the logical operators and (and) and | (or), allowing for sophisticated selection criteria. Furthermore, LAGRANTO.ocean can check whether a trajectory crosses a certain region during a specific period or travels a minimum length. Again, the man page provides additional examples and a complete list of available parameters.
c. Boundary conditions
The interaction with the bathymetry poses a new challenge and is the major new component added to LAGRANTO. In LAGRANTO.ocean, we use no-normal and no-slip boundary conditions. No-normal flow conditions are implemented by setting the component normal to the bathymetry to zero at grid points adjacent to the bathymetry [u(i, j) in Fig. 1]. No-slip boundary conditions are implemented by setting the first grid point outside the valid domain to opposite of the velocity at the adjacent grid point inside the valid domain. In Fig. 1, for a north–south boundary, this corresponds to replacing the meridional velocity υ(i, j −1) = υ(i, j). Consequently, a linear interpolation in the zonal direction between the meridional flow υ(i, j −1) and υ(i, j) yields zero current speed at the boundary. Together with the no-normal conditions this yields zero flow at all boundaries. The procedure can be applied before the original C grid is destaggered (see below). After destaggering the grid we may replace u(i, j −1) = u(i, j) to obtain zero velocity at the boundary for linear interpolation in the zonal direction.
The combination of no-normal and no-slip boundary conditions consequently yields zero flow at all boundaries. The boundary conditions apply only to the first grid point inside the valid domain. If one trajectory still interacts with the bathymetry, then it is removed from the output and analyzed only until the time it interacts with the boundaries. An alternative treatment, used by Koszalka et al. (2013), for example, would be to start a new trajectory at the intersection point using only the velocity component parallel to the bathymetry. No-normal boundary conditions might also be extended into the valid domain, where they steadily weaken with increasing distance to the bathymetry (i.e., damping of the normal flow). At this stage, however, we decided to stop the trajectory computation in the case of their interception by the boundaries, as it is not clear how reasonable such a damping boundary condition is for the normal flow. However, we plan to implement various options to address this issue from which the user can choose. For a detailed discussion on the role of the boundary conditions in ocean modeling in general, see, for example, Moore (2004).
d. Preprocessing ORAS5 data
Theoretically, the trajectory computation in LAGRANTO is independent of the input data’s coordinate system (e.g., ORCA grid) or the staggering (e.g., Arakawa C or B). However, LAGRANTO.ocean currently requires the input data to be on a regular grid, that is, the latitude and longitude arrays are one-dimensional and monotonically increasing. In this study, the data are interpolated and rotated from the native tripolar grid to a regular geographical grid with a horizontal resolution of 0.25° everywhere. At the same time, we also destagger the data. The vertical spacing of the model levels is kept constant. All fields are available as 5-day averages. Döös et al. (2011) noted that water parcel trajectories based on 5-day averages lead to Lagrangian time scales that are nearly 2 times greater compared to higher-resolution-model output. However, it remains unclear how much of this difference is due to the ocean model they used in question and how much is due to an underestimation of the Lagrangian time scale by ocean drifters, which are not perfect Lagrangian devices (for a discussion on the observed differences between observed and computed Lagrangian data, see Garraffo et al. 2001; McClean et al. 2002; van Sebille et al. 2011; Nilsson et al. 2013).
4. Case-study applications
The following case studies are chosen to reflect a selection of phenomena on different spatial and temporal scales while retaining a degree of simplicity in order to facilitate a hands-on experience.
a. Case I: Lagrangian flow through a Gulf Stream cold-core ring
The study of vertical motion in ocean eddies and along ocean fronts is vital (e.g., Nardelli 2013; Pidcock et al. 2013; Pallàs-Sanz et al. 2010; Mahadevan et al. 2008; Nagai et al. 2006) and several simplified linearized conceptual models for vertical transport in ocean eddies exist (McGillicuddy et al. 1998, 2007; Siegel et al. 2011). However, studies that take the nonlinear submesoscale dynamics into account have highlighted the complex nature of vertical motion inside ocean eddies (Martin and Richards 2001; Mahadevan et al. 2008; Nardelli 2013) that is not entirely captured by the simplified conceptual models. For example, Nardelli (2013) identified, based on a combination of satellite observations and semigeostrophic diagnostics, vertical oscillations along parcel trajectories flowing through a closed vortex as the dominant source of vertical motion in ocean eddies. Hence, it can be insightful to study the mesoscale dynamics of ocean eddies from a Lagrangian perspective. Suppose that we are interested in the flow through a Gulf Stream cold-core eddy and how the water parcels encircle the eddy and change their depth. For such a scenario LAGRANTO can help us to study vertical motion in ocean eddies in great detail.
In this example, from January to March 2010, we show that the 5-day-averaged mean flow allows us to depict aspects of the mesoscale circulation in an ocean eddy. The backward trajectories start on 2 March 2010 at the center of the cold-core eddy at 38°N, 62°W from all grid points in the horizontal inside a circle with a radius of 50 km and on all grid points between the surface and a depth of 2000 m. The corresponding LAGRANTO.ocean command to define the starting cylinder is
startf.ocean 20100302_00 startpoints.txt 'circle.grid(-62,38,50)@grid(0,2000)@m'.
Afterward, we use caltra.ocean to compute the trajectories backward from the starting positions defined in startpoints.txt for a period of 60 days:
caltra.ocean 20100302_00 20100101_00 startpoints.txt traj.txt.
From the ensemble of trajectories (~500 trajectories) now stored in traj.txt, we select those that change their depth between the first and the last time steps by at least 100 m. The corresponding LAGRANTO.ocean command is:
select.ocean traj.txt traj-refined.txt 'GT:depth(CHANGE):100:FIRST,LAST'.
Here, CHANGE identifies a change in the depth in any direction, that is, it is not clear whether the trajectories ascend or descend. If we want to limit the backward trajectories to those that descend, we could use DIFF instead. Rather than using the entire time span (FIRST,LAST), it is possible to specify specific time steps (h). Let us now have a closer look at the selected trajectories.
The refined bundle of trajectories is located on 6 January 2010 at a depth of ~1000 m (Fig. 2a) and moves southward while descending along the westward flank of the cold tongue along steep isopycnals (Fig. 2b). On the eastward flank of the elongated tongue, the trajectories ascend (Fig. 2b) before wrapping cyclonically around the center of the developing cold core. Because the sloped isopycnals form a “peninsula” with a step front toward the warmer side of the Gulf Stream, parcels that move cyclonically around the tongue necessarily slant upward when moving northwestward on its eastern flank. Similarly, on the westward flank southward-moving parcels slant along the downward-tilted isopycnals (Fig. 2c). When the cold-core ring forms a closed system, the trajectories circle cyclonically around its center while ascending, consistent with the near-surface divergence associated with cold-core rings (Richardson 1983). An animated version of the trajectories is provided in the online supplement (JTECH-D-16-0198.s1). A follow-up research question could study similar trajectories in a fully developed and closed eddy.
Another branch of eddy-related research centers on the two-way interaction between ocean eddies and the atmosphere (e.g., Williams 1988; Small et al. 2008; Cerovĕcki and Marshall 2008; Chelton 2013; Frenger et al. 2013; Gaube et al. 2015), if we are interested in changes in buoyancy and the role of air–sea fluxes. The heat loss at the surface above an ocean eddy is an important contribution to the heat budget and can cancel the effect of heat convergence in the mixed layer (Cerovĕcki and Marshall 2008). To better understand buoyancy changes in ocean eddies, it can be insightful to supplement a budget analysis with a trajectory computation. After tracing potential density (named sigmat on the preprocessed ORAS5 data) and other variables along the previously computed large ensembles of trajectories, we select those that exhibit a strong change between the first and the last time steps, that is,
select.ocean traj.txt sigma-refined.txt 'GT:sigmat(DIFF):0.4:FIRST,LAST'.
Using this command, we select from the large ensemble of backward trajectories those that experience an increase in potential density of more than 0.4 kg m−3. The near-surface trajectories identified in this way wrap cyclonically around the developing cold-core ring (Figs. 3a,b). The motion is quasi horizontal, which is consistent with a weaker slope of the isopycnals. The trajectories encircle the cold core (Fig. 3c) one additional time during the same period as when the previously selected deeper trajectories located at ~1000-m depth continued directly south (Fig. 2). Between January and March 2010, the water parcels experience an increase in potential density from 26.2 to 26.6 kg m−3 (Fig. 4a). The first step increase between 6 and 16 January is accompanied by a strong reduction in potential temperature from 17.7° to 16.5°C (Fig. 4b) with almost no change in salinity (Figs. 4c,d). The net (downward) heat flux, which is traced at the horizontal position of the trajectories, indicates a pronounced heat loss of the ocean during this period (Fig. 4e), which ultimately is the diabatic effect that modifies the potential temperature along the trajectories. A second pronounced increase in potential density between 31 January 2010 and 20 February 2010 is driven predominantly by the increasing salinity (Fig. 4c) from 36.10 to 36.32 psu. This is reflected in the way that the temperature–salinity (T–S) trajectory moves to higher salinity and potential density values while changes in potential temperature are weak (Fig. 4d). The net (upward) water flux during this period, computed from evaporation minus solid and liquid forms of precipitation, indicates evaporation from the ocean surface (Fig. 4f), which ultimately drives the increase in salinity along the trajectories.
This example highlights the strength of the trajectory computation because it allows us to pinpoint the different mechanisms that are responsible for the changes in potential density along the flow through the cold-core ring. Furthermore, parcel trajectories based on 5-day-averaged current are able to represent the mesoscale flow through the developing Gulf Stream eddy. Potential future research may address the interplay between the Gulf Stream eddy and the atmosphere from a Lagrangian perspective.
b. Case II: Sources of dense water on the East Greenland shelf
The dense waters that form in the Nordic seas contribute to the lower limb of the global thermohaline circulation, affecting the northward ocean heat transport and the high-latitude climate (Bacon 1998; Rahmstorf 2002). However, the linkage from the dense water formation to the Greenland–Scotland ridge overflows and the deep western boundary current transports is not one to one. The latter two show variability on short time scales but no seasonal or longer-term variability in observations (Jochumsen et al. 2012, 2015; von Appen et al. 2014; Harden et al. 2016), whereas the heat loss fueling the deep-water formation varies on seasonal and interannual time scales. Furthermore, the East Greenland shelf, west and southwest of the Denmark Strait, is an additional source region of waters with potential density anomalies larger than 27.8 kg m−3, the common definition for the Denmark Strait Overflow Water. Rather than flowing through the Denmark Strait, these dense waters contribute to the Denmark Strait Overflow plume through the East Greenland Spill Jet (Pickart et al. 2005; Koszalka et al. 2013; von Appen et al. 2014). In this example, we highlight the formation of dense waters on the Greenland shelf with an ensemble of 180-day backward trajectories released in the vicinity of the Denmark Strait.
To this end, we compute 180-day backward trajectories released at every grid point in a cyclinder centered on the Denmark Strait between a depth of 200 and 700 m on 25 February 2010 (Fig. 5). The corresponding LAGRANTO.ocean command is
startf.ocean 20100225_00 startpoints.txt 'circle.grid(-29,66.25,500)@grid(700,2000)@m'.
After computing the trajectories using caltra.ocean and tracing additional variables along the flow using trace.ocean, we select from the large ensemble of trajectories those that on 25 February 2010 have potential density greater than 27.8 kg m−3. We recognize that this threshold for the overflow water may be different in the model compared to what is found in the real ocean. However, for simplicity we retain this widely used threshold of 27.8 kg m−3. The corresponding one-line command reads:
select.ocean trajectories.txt overflow.txt 'GT:sigmat:27.8:FIRST'.
The analysis reveals three distinct regions at the East Greenland shelf with waters dense enough to contribute to the Denmark Strait Overflow plume (Fig. 5). The first region is located to the north of the Denmark Strait (Fig. 5 at A) and has its source region further upstream along the East Greenland Current. To the south of the Denmark Strait, there are two more regions (Fig. 5 at B and C). The northern region is located west of the Denmark Strait sill along the East Greenland shelf (Fig. 5 at B). Here, the trajectories indicate predominately recirculating flow. Excluding with a regional identification criterion the recirculating trajectories at B (Fig. 5), it becomes clear that this area receives waters from north of the Denmark Strait and from the Irminger Sea (Figs. 6c,d). Note that we analyze trajectories started at only one time step. Accordingly, the recirculating trajectories may have arrived at B from north of the Denmark Strait or from the Irminger Sea prior to the analyzed 180-day period. The southern region is located south of the Ammassalik Fjord around 65°N and 38°W (Fig. 5; C), with source waters derived primarily from the Irminger Sea (Fig. 6b at C). In T–S space, the water mass in region B lies on the line connecting its two source water masses (Fig. 6a; B1 and B2), which indicates that the final water mass properties are a result of mixing of the two source water masses. In region C, the relatively homogeneous source water mass properties are modified by winter cooling, which increases the mean density to above 27.8 kg m−3 (Figs. 6a,b; C). Region B is similar to the East Greenland Spill Jet source region recognized in previous studies (Koszalka et al. 2013; von Appen et al. 2014), whereas region C has not previously been reported as a source region for dense waters. A further assessment of the importance of the water mass transformations in the southern region is beyond the scope of this study, but we note that the water mass properties in this region are likely to have a seasonal signal because they depend on surface cooling and might therefore be present only at the shelf during winter and spring.
In this example, we have illustrated some of the aspects of the dense water circulation that can be easily studied using LAGRANTO.ocean and the related processes thereof. However, the usefulness of the trajectory calculations of course depends on the appropriateness of the underlying dataset. For example, in our case the 5-day averaging leads to much slower-than-observed propagation time scales (Döös et al. 2011). Furthermore, for more robust results, the trajectory calculations should be performed over longer times rather than the one-time release (here, 25 February 2010) shown in this example.
c. Case III: Persistence of seasonal sea surface temperature anomalies
Next, we move to even longer time and spatial scales and suppose we are interested in studying air–sea interactions and the process involved in the formation of SST anomalies. Studies have shown that strong SST anomalies may persist over seasons or reemerge and moderate the atmospheric circulation (Alexander and Deser 1995; Kushnir et al. 2002; de Coëtlogon and Frankignoul 2003; Deser et al. 2003; Ciasto and Thompson 2009). During the 2009/10 Southern Hemispheric summer season (DJF), the South Pacific experienced a record SST anomaly (five times the local standard deviation) that was formed by a persistent anticyclone, which itself is believed to be enhanced by the strong central-Pacific El Niño event that occurred during this period (Lee et al. 2010). In this section, we use this event to address (at least partially) with the help of LAGRANTO.ocean the following open question: How long can we measure surface temperature anomalies, not locally but rather along the history of the affected near-surface water parcels?
To this end, forward and backward trajectories are calculated for 1 yr starting on 27 December 2009 from all grid points inside a 100-km radius centered at 45°S, 125°W at a depth of 0.5 m, which is where the SST anomaly during this period is strongest. Further, we refine the trajectories to those that are located throughout the integration period (one year forward and one year backward) between the surface and a depth of 1 m to make sure they stay in contact with the atmosphere throughout this period. The Lagrangian analysis allows for monitoring the evolution of temperature without the need to account for horizontal advection.
Along the flow we trace potential temperature, salinity, potential density, and the surface heat flux. Additionally we trace along the flow the corresponding monthly climatology computed from 30-yr ORAS5 data. At every time step the climatological values are interpolated to the parcel’s position similar as for the full fields.
The selected near-surface water parcels are located to the south and southeast 1 yr before arriving at 45°S, 125°W (Fig. 7) from where they move northward from January 2009 to April 2009 and then eastward until December 2009. Along-track potential temperature and heat flux evolve in agreement with the seasonal cycle (Figs. 8a,b,c,d). Both are strongly correlated at a one-month lag to a one-half-a-month lag with the heat flux leading the temperature evolution. The heat flux, which is one forcing term of the along-flow temperature tendency (Figs. 8c,d; green contour), and the temperature tendency are highly correlated. During early March 2009, the heat flux and the temperature tendency turn negative and the ocean starts to lose heat. Shortly after, we observe the maximum in the along-parcel potential temperature during the first half of the year 2009 (+11 °C) before the temperature reacts to the surface heat loss and along-flow potential temperature starts to decrease. By early September the heat flux and temperature tendency turn positive again (Fig. 8c), and we observe the minimum in along-parcel potential temperature (+8 °C) during the year 2009 (Fig. 8a). The Lagrangian analysis would allow us to compute the integrated heat gain along the water parcels and hence allows for quantifying the heat flux’s contribution to the observed warming. The result (not performed in this study) could then be used to compare the heat flux contribution to other diabatic processes that are able to modify potential temperature along the flow, for example, mixing.
Next we compare the along-flow evolution of potential temperature to the monthly climatology. The interpretation of along-flow anomalies, however, is more complex than the interpretation of the full fields because the anomalous flow speed or direction emerges as an along-flow temperature anomaly. Throughout the year 2009 the water parcels are close to the climatological value until November 2009 (Fig. 8a). From November onward the ocean parcels warm strongly above the local climatology (Fig. 8a). The trajectory computation allows for quantifying the along-flow heat gain plus the anomalous heat gain. The along-flow integrated heat flux anomaly between September 2009 and the end of 2009 indicates an additional heat gain of 0.7 GJ m−2 compared to normal conditions. Because the evolution of the heat flux leads the temperature evolution, and the temperature tendency and the heat flux are strongly correlated, it appears plausible that the corresponding along-flow anomaly in potential temperature is driven by this additional heat flux. The corresponding T–S diagram (Fig. 8c) shows that the decrease in potential density from the end of October 2009 to the end of the year 2009 (~1 kg m−3) is solely due to the potential temperature increase at constant salinity.
After being warmed in late 2009 to potential temperatures of approximately 15°C, the ocean parcels advect in a relatively confined stream northeastward during the following year, 2010 (Fig. 7). During this period they are warmer than the climatology for 5 more months until around May 2010. (Fig. 8b). The heat flux fluctuates during this period around the monthly climatology (Fig. 8d). In particular during the period between January and May 2010, when the ocean parcels adjust slowly their temperature to the climatology, there is no clear additional heat gain or heat loss compared to normal conditions. Consequently, the heat flux is unlikely the main driver underlying the temperature adjustment toward the climatology and rather drives only the seasonal cycle.
After the parcel temperature adjusted back to normal values during May 2010, the parcel’s potential temperature remains close to the climatological values until November 2010 (Fig. 8b). Afterward, the trajectories rapidly warm to potential temperature anomalies of approximately 3°C above the local climatology between November and December 2010 (~19°C). This is in agreement with an along-flow integrated heat flux anomaly during the period from October to December 2010 of approximately 0.8 GJ m−2 (Fig. 8d). Because the heat flux again leads the development of the temperature, and because of the high correlation between the heat flux and the along-flow temperature tendency (Figs. 8c,d; green contour), it is may also be the driver underlying the observed temperature anomaly. During both anomalous warming events in late 2009 (Fig. 8a) and late 2010 (Fig. 8b), an anomalous heat gain of approximately 0.7–0.8 GJ m−2 is associated with an anomalous temperature increase of approximately 2°–3°C above the local climatology though the water parcels are located at different locations.
Noteworthy are also two warming events during June 2009 and June 2010, both of which are reflected in an increase in along-flow temperature tendencies (Figs. 8c,d; green contour). Absolute temperature changes are, however, weak but are observable (Figs. 8a,b; note the two small kinks in the purple temperature curve). Overall, we should not expect a similar response of the temperature during summer compared to winter when the ocean is cooler and the mixed layer is deeper.
In summary, we argue that ocean trajectories can provide, in general, useful and complementary insights into the persistence of potential temperature anomalies and other corresponding anomalies (in the above-mentioned example, surface anomalies) and can be used to supplement, for example, a heat budget analysis that quantifies anomalous advection in the ocean and in the atmosphere. In this example, it takes about 5 months after the maximum in the along-flow temperature anomaly for the parcels to adjust their temperature back to the climatology at their new position with only a little contribution from the anomalous heat flux. Further studies may link this Lagrangian view with the mixed-layer response time to surface heat exchange or shed more light on the interplay with the atmosphere.
5. Concluding remarks and outlook
This study has introduced and highlighted some examples of computing Lagrangian water parcel trajectories in a high-resolution ocean reanalysis dataset and documented the structure and handling of LAGRANTO.ocean. Three, albeit brief, examples have been presented as showcases to illustrate that LAGRANTO.ocean allows, in a straightforward way, the following:
Identifying coherent flow features associated with the mesoscale ocean circulation. The presented case study accurately depicts vertically alternating cyclonic motion through a mesoscale Gulf Stream cold-core ring.
Following water mass transformation, and the underlying processes (salinity vs temperature change), along the trajectories. Here we presented a case study where we identified water masses that potentially feed the East Greenland Spill Jet.
Tracing water parcels that form an extreme seasonal sea surface temperature anomaly and quantify the ocean parcel’s amount of received heat by integrating heat flux (anomalies) along the flow.
Based on the three case studies, we argue that parcel trajectories provide insights that are complementary to Eulerian-based approaches. For example, in case study II, ocean trajectories allow us to better pinpoint the origin of dense water, its sources, and sinks. From case study III, we could quantify the time that it takes for water parcels in this specific example to “recover” from anomalous surface temperatures. Accordingly, we argue that water parcel trajectories based on a high-resolution ocean reanalysis can provide helpful insights into the nature of persistent temperature or salinity anomalies. Finally, the 5-day-average current fields archived in ORAS5 are—as noted earlier—only marginally adequate for Lagrangian studies of deep convection. Moving to daily averages would greatly enhance the usefulness of future ocean reanalyses.
The current bottleneck in LAGRANTO’s performance is the I/O. To prepare the code for the next generation of climate model (CMIP6) and reanalysis data (ERA5), we plan to parallelize the netCDF interface of the code. The actual trajectory computation is much less demanding compared to the read and write of netCDF data. Furthermore, as discussed earlier, we plan to provide the user with the possibility of choosing between different boundary conditions.
Wernli and Davies (1997a) introduced LAGRANTO approximately 20 years ago. The code has been intensively revised by Sprenger and Wernli (2015), and with this publication, it is now moving into the ocean (www.lagranto.ethz.ch). Sebastian Schemm acknowledges funding from the Swiss National Science Foundation (Project P300P2_167745). Øyvind Breivik acknowledges funding from the Research Council of Norway through the projects RETROSPECT (Grant 244262) and CIRFA (Grant 237906). The ORAS5 reanalysis data used in this study are available through the ECMWF ECFS storage. We thank Magdalena Balmaseda and Hao Zuo for their support with the ORAS5 data, and we thank two anonymous reviewers for their constructive comments, which helped to improve our manuscript.
Supplemental information related to this paper is available at the Journals Online website: http://dx.doi.org/10.1175/JTECH-D-16-0198.s1.
Current affiliation: Department of Earth and Planetary Sciences, The Johns Hopkins University, Baltimore, Maryland.