the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
FEOTS v0.0.0: a new offline code for the fast equilibration of tracers in the ocean
Wilbert Weijer
Jiaxu Zhang
In this paper we introduce a new software framework for the offline calculation of tracer transport in the ocean. The Fast Equilibration of Ocean Tracers Software (FEOTS) is an endtoend set of tools to efficiently calculate tracer distributions on a global or regional subdomain using transport operators diagnosed from a comprehensive ocean model. To the best of our knowledge, this is the first application of a transport matrix model to an eddying ocean state. While a Newton–Krylovbased equilibration capability is still under development and not presented here, we demonstrate in this paper the transient modeling capabilities of FEOTS in an application focused on the Argentine Basin, where intense eddy activity and the Zapiola Anticyclone lead to strong mixing of water masses. The demonstration illustrates progress in developing offline passive tracer simulation capabilities, while highlighting the challenges of the impulse response functions approach in capturing tracer transports by a nonlinear advection scheme. Our future work will focus on improving the computational efficiency of the code to reduce timetosolution, using different basis functions to better represent nonlinear advection operators, applying FEOTS to a parent model with unstructured grids (Ocean Model for Prediction Across Scales, MPASOcean), and fully implementing a Newton–Krylov steadystate solver.
Many oceanographic research problems involve the transport and distribution of tracers that do not feed back on the ocean dynamics. Examples of such problems are the diagnostic tracking of water masses using passive tracers (e.g., Dukhovskoy et al., 2016; Zhang et al., 2021), validating the use of isotopes or grain size distributions in marine sedimentary records to infer past ocean circulation changes (e.g., Jahn et al., 2015; Zhang et al., 2017; Gu et al., 2019; Missiaen et al., 2020), assessing anthropogenic carbon uptake by the ocean (e.g., Sarmiento et al., 1992; Khatiwala et al., 2009; Wang et al., 2012), studying the evolution of marine biogeochemical systems (e.g., Séférian et al., 2020), or tracking the fate of microplastics in the ocean (e.g., Mountford and Morales Maqueda, 2019). In many cases, the transport of these tracers is dominated by mesoscale processes like eddies. Typical ocean climate models use grids that are too coarse to explicitly resolve these processes and rely on parameterizations to simulate their impact on tracer fields, but it has become clear that these eddyparameterized models fail to reproduce some critical aspects of the real ocean (e.g., Lozier, 2010). What is more, the ocean also contains dynamical features that rely on mesoscale eddies and which cannot be reproduced by lowresolution models.
A case in point is the Argentine Basin. It is among the most turbulent regions in the world ocean (Fu and Smith, 1996), mostly on account of the confluence of two western boundary currents, the Brazil and Malvinas currents (e.g., Garzoli, 1993). A seamount in the center of the basin is associated with a local minimum in eddy kinetic energy (Fu and Smith, 1996) but is surrounded by a very strong barotropic vortex, the Zapiola Anticyclone (ZA; Saunders and King, 1995; de Miranda et al., 1999). This anticyclone is understood to be driven by the intense eddy field (Dewar, 1998) and has been shown to inhibit exchanges between the interior of the ZA and its surroundings (Weijer et al., 2015, 2020). Strikingly, the Argentine Basin is also a main conduit of water mass exchange between the Atlantic and Southern oceans (e.g., Jullion et al., 2010), so eddydriven mixing of water masses in the Argentine Basin may have implications for the global thermohaline circulation.
It is clear that studying tracer transport and mixing in the Argentine Basin requires an ocean model that resolves the ocean's mesoscale, not only to accurately represent the transports by narrow boundary currents and the turbulent eddy field but also to generate the ZA in the first place. Second, the equilibration of tracers at deeper levels may take many decades or centuries, requiring a model capability that can be run for such long times or that can determine equilibria directly using iterative solvers. Third, this problem is ideally addressed using a representation of the global ocean circulation, not only so that the ocean circulation in the region of interest is fully consistent with largescale oceanographic drivers but also for the practical reason that not every regional problem would require a unique model configuration.
The cost and technical challenges of running global and fully dynamic ocean models make it expensive and difficult at best, but often impossible, to address problems like these. It is obvious that there is a need for simple and efficient tools that solve tracer transport and distribution problems without the need to explicitly simulate ocean dynamics. In the past decade or so, several transport matrix models (TMMs) have been developed that allow for the simulation of passive tracers in a standalone (offline) code, using transport operators that have been diagnosed from comprehensive ocean models (Primeau, 2005; Khatiwala et al., 2005; Khatiwala, 2007; Bardin et al., 2014; Kvale et al., 2017; Zanna et al., 2019; Chamberlain et al., 2019). In particular, Khatiwala et al. (2005) pioneered the approach of empirically estimating an ocean model's transport processes by diagnosing the action of the transport (advection and diffusion) operators on simple basis functions like impulse fields. The resulting impulse response functions (IRFs) can straightforwardly be converted to the desired transport matrices. This process requires minimal intervention in the parent model.
Despite the success of these techniques, TMM models have generally only been applied to lowresolution, noneddying ocean states. This makes them unsuitable to study problems where eddies play a zerothorder role in tracer transport and ocean dynamics, like the Argentine Basin. To the best of our knowledge, the capability introduced in this paper is the first TMM model developed for and applied to eddying ocean states, with transport operators diagnosed from a global ocean model with nominal resolution of 0.3^{∘} and 100 vertical levels (∼10^{8} df's).
Our computational framework is embodied in the Fast Equilibration of Ocean Tracers Software (FEOTS; https://github.com/LANL/FEOTS). FEOTS is based on the methodology of Bardin et al. (2014) but is specifically designed to tackle the large computational problems associated with tracer transport in a global eddying ocean. In particular, FEOTS is written in ObjectOriented Fortran instead of MATLAB, which makes it straightforward to port the code to supercomputers, to parallelize the workflow, and to expose the code to a plethora of existing solver libraries. Also, FEOTS provides an endtoend set of tools that streamlines the process of building and running an offline tracer model. It (i) uses an advanced optimization algorithm to generate an optimal set of impulse functions, given the grid layout and operator stencils for the parent model; (ii) transforms the resulting IRFs obtained from the parent model to transport operators; (iii) sets up regional or global tracer problems with different types of tracers; (iv) runs forward simulations of the offline tracer model; and (v) uses a Newton–Krylov solver to determine steady tracer distributions. To date, the authors have implemented the first four of these features, and their implementation, validation, and verification are presented in this paper. We plan to fully implement a Newton–Krylov steadystate solver in future releases. The framework is applied here to the Parallel Ocean Program (POP; Smith et al., 2010), but the design is flexible enough to be generally applicable – in particular to the new generation of ocean models with unstructured grids, like the Ocean Model for Prediction Across Scales (MPASOcean; Ringler et al., 2013). Other future work aims to improve timetosolution by exposing data parallelism in forwardsimulation components of FEOTS.
In this paper we present validation and verification results for the FEOTS offline tracer solver in a regional forwardsimulation configuration focused on the Argentine Basin. Specifically, we will show that uniform tracer fields are preserved within 0.01 % after 5 years of offline model integration, and we will present a comparison of the offline regional tracer simulation using 5 d averaged transport operators with an online tracer simulation using the parent model. Finally, we comment on the compute costs for running offline simulations with FEOTS.
Our work builds on the methodologies of Bardin et al. (2014) to create FEOTS. FEOTS comes with tools to generate impulse fields, translate impulse response fields to sparse matrices corresponding to advection and lateral diffusion, create vertical diffusion operators from eddy diffusivities reported by the parent model, and execute offline regional transient tracer simulations. In Sect. 2.1 we discuss the parent model POP, which is used for this study. In Sect. 2.2 we present the governing equations for a noninteracting passive tracer system and outline the methodology for capturing transport operators from a comprehensive ocean model in Sect. 2.3. The challenge of using a fluxlimited advection scheme is discussed in Sect. 2.4. The offline forwardstepping algorithm and treatment of vertical mixing is presented in Sect. 2.5. Last, the constant preservation test problem is defined in Sect. 2.6, and the Argentine Basin test problem is defined in Sect. 2.7.
2.1 Parent model
FEOTS is used for offline tracer simulations using transport operators diagnosed from a comprehensive ocean model, the “parent model”. Here we apply FEOTS to the Parallel Ocean Program, in the framework of the E3SMv0HiLAT climate model (Hecht et al., 2019). Our specific configuration is described in Zhang et al. (2019) and is referred to as E3SMv0HiLAT03. E3SMv0HiLAT03 has a tripole grid with nominal 0.3^{∘} spatial resolution. The grid has 1200×800 grid cells and 100 levels in the vertical. The grid has a “seam” in the Arctic that connects the poles in Siberia and Canada. Although technically not eddyresolving in most of the world ocean (Hallberg, 2013), this configuration has a vigorous eddy field (Zhang et al., 2019) and a realistic representation of the eddydriven Zapiola Anticyclone in the Argentine Basin (Weijer et al., 2020).
With O(10^{8}) df's, and an eddying field that requires transport operators at high temporal frequency (e.g., 5daily averages), the data volume of the diagnosed IRFs can become unmanageable quickly. It is therefore critical to keep the number of IRFs to an absolute minimum. This means that we need to choose advection and diffusion treatments that have the most compact stencils. Typically highresolution ocean models use a biharmonic mixing scheme to dampen out the dispersive errors caused by the advection operator (e.g., Hecht et al., 2008). Biharmonic mixing is clearly undesirable for our application, given its huge stencil. Of the three advection schemes currently implemented in POP (Smith et al., 2010), the centered and thirdorder upwind schemes require explicit diffusion to manage the dispersion error. The fluxlimited Lax–Wendroff scheme does not, making this the only reasonable choice for our application, even though its stencil is quite large (27 grid points, $\mathrm{3}\times \mathrm{3}\times \mathrm{3}$) compared with the other two (7 grid points for centered, 13 for the thirdorder upwind scheme). This choice requires 53 impulse functions to capture the advection operator (compared to 34 impulse functions for the thirdorder upwind scheme) but eliminates the need for explicit diffusion.
The model is forced by the normalyear Coordinated OceanIce Reference Experiments version 2 (COREII; Griffies et al., 2012) climatology, which has been a widely used framework to force ocean and/or sea ice models for hindcast simulations. With a time step of 7 min the model typically yields maximum CFL values of O(10^{−1}) or smaller. Although the model was run for 186 years, we diagnosed the transport operators for the 5year period starting at simulation year 64. Even though 63 years of spinup is not sufficient to fully equilibrate the stratification in the deep ocean, the main circulation features (e.g., boundary currents, the eddy field, the Zapiola Anticyclone) are well established by then, making this an appropriate data set to demonstrate the capability of FEOTS. We refer to Weijer et al. (2020) for evaluation of the hydrography and circulation in the Argentine Basin in a companion simulation.
2.2 Governing equations
We model a passive dye tracer as a concentration field that is subjected to advection and diffusion,
where v is the fluid volume anomaly, c is the tracer concentration, t is time, u is the ocean velocity field, and 𝕂 is the diffusivity tensor that models unresolved eddy activity. The fluid volume anomaly is a unitless quantity that is a measure of the relative change of the fluid volume due to movement of the free surface. This formulation is chosen so that the total tracer is conserved and so that the offline model is consistent with the parent model (Smith et al., 2010).
In practice, the fluid volume anomaly is defined as the ratio of the freesurface height to the uppermost grid cell; only in the uppermost grid cell
where η is the freesurface height, dz_{1} is the grid cell thickness in the uppermost grid cell, and δ_{k,1} is the Kronecker delta function.
The initial and boundary conditions are set to be
and
where z is depth (measured positive downward), θ is longitude, and ϕ is latitude. The semidiscrete form of Eq. (1) can be written as
where c is a vector of the discrete values of the tracer, A is the advection matrix, D_{h} is the horizontal diffusion matrix, and D_{v} is the vertical diffusion matrix. In the examples presented in this paper, the advection matrix corresponds to the fluxlimited Lax–Wendroff advection scheme on an Arakawa B grid; this scheme is chosen in the parent model configuration. Since we do not enable explicit lateral tracer diffusion in the parent model in this study, all elements of D_{h} are zero. To create the vertical diffusion matrix, we diagnose diffusivity coefficients, which are generated from the KPP parameterization, and create a vertical mixing operator offline using centered differencing with noflux conditions applied at the ocean free surface and at the seafloor.
2.3 Graphcoloring approach to operator diagnosis
The purpose of capturing the impulse response functions is to diagnose sparse matrices that are consistent with the advection discretization in the parent model. To illustrate this procedure, suppose that the parent model has N grid cells and that N impulse fields are set as the Kronecker delta functions:
In this setup, impulse function i is zero at all grid cells except for grid cell i where the impulse function has a value of one. Application of the transport matrix to each of the c^{(i)} returns column i of A,
While using a set of Kronecker delta functions will completely diagnose all of the elements of the transport matrix, this strategy is computationally expensive. For each time step this strategy requires computing the advective tendency for N tracer fields, where N is the number of grid cells. For example, coarseresolution models at O(1^{∘}) resolution have roughly 10^{6} grid cells. Storing 10^{6} impulse and impulse response functions would require approximately 3 TB of memory at double precision.
To reduce the number of required impulse functions to fully diagnose the transport matrices, we can take advantage of the fact that the advection scheme results in a sparse matrix. Equivalently, the domain of influence of the advection operator is limited to nearby grid cells. The parent model employed in Bardin et al. (2014) used a thirdorder upwind scheme, where the impulse response is guaranteed to extend no further than two grid cells in each spatial dimension, giving a $\mathrm{5}\times \mathrm{5}\times \mathrm{5}$ brick for the domain of influence. Because of this, the authors used a set of 125 tracer fields,
FEOTS offers a unique capability to generate a minimal set of impulse functions by posing the problem as a graphcoloring problem. A graph G(V,E) is defined by a set of vertices V and edges E that connect the vertices. Two vertices connected by an edge are said to be adjacent. A valid graph coloring of G(V,E) assigns colors to each vertex so that no two adjacent vertices have the same color. To calculate impulse functions that can be used to diagnose transport operators, FEOTS offers functionality to express a POP mesh and an advection stencil into an equivalent graph that is colored with a greedy algorithm. This formulation has the benefit that it can be generalized to parent models based on unstructured grids, and it takes into account irregular boundaries from variable bathymetry.
In FEOTS, graph vertices V correspond to each ocean grid cell, centered on tracer points, in the POP mesh. Two vertices are adjacent if their impulse response functions overlap. Because a valid coloring results in adjacent vertices having distinct colors, vertices with the same color can safely be assigned to the same impulse function. Consequently, the chromatic number of the graph corresponds to the number of impulse functions used for model diagnosis. For this work the parent model uses a 0.3^{∘} periodic tripole mesh and the thirdorder fluxlimited Lax–Wendroff advection scheme. This approach results in 53 impulse functions required to uniquely diagnose the transport operators.
The transport operators in Eq. (5) are diagnosed empirically from the parent model, using the methodology used by Bardin et al. (2014) and pioneered by Khatiwala et al. (2005). This process involves diagnosing and time averaging the impulse response functions corresponding to each of the 53 impulse fields in the passive tracer equations in POP. Modifications are made in POP to initialize the passive tracers to the impulse functions at the beginning of each time step. The impulse response function is set equal to the advective tendency as diagnosed in POP. After diagnosing the IRF in each time step, the tracer field is reset to the impulse field so that the IRF at the next time step is saved. The IRFs are timeaveraged over a configurable averaging period and written out to file.
For each timeaveraged IRF, we create a sparse matrix representation of the advection operator in POP. Modeling the advection operator as a matrix–vector multiplication, Eq. (7) shows that passing an impulse at grid point i through the advection operator returns column i of the matrix that corresponds to the advection. We use this to construct a sparse matrix in compressedrowstorage format that corresponds to the advection operator.
For our test problems, we diagnosed the 5 d averaged IRFs and vertical diffusivities for the 5year analysis period of the parent model. We repeated the simulation for 105 d, diagnosing 1 d averaged IRFs. With this methodology and the 7 min time step, the 1 d averaged operators are each an average of 1440 IRF snapshots, and the 5 d averaged operators are each an average of 7200 IRF snapshots. The data volume of the global parent model's 5 years' worth of 5 d averaged operators (365 IRFs and diffusivities) is about 9 TB. Once transformed to transport operators, the data volume is 4 TB.
2.4 Fluxlimited advection
The parent model uses a nonlinear fluxlimited Lax–Wendroff advection scheme; the flux limiter is equivalent to the ULTIMATE flux limiter described in Leonard (1991) and Hundsdorfer and Trompert (1994). Because FEOTS treats offline advection as a linear operation, online advection with the parent model, in this case, can not be equivalent. In what follows, we illustrate how the impulse functions result in the diagnosis of a more diffusive advection scheme than Lax–Wendroff.
In POP, the tracer equations are discretized using a finite volume approximation on an Arakawa B grid. The nonlinear fluxlimited Lax–Wendroff scheme is applied in three dimensions using a dimensional splitting technique, the details of which can be found in Chap. 6 of Smith et al. (2010). For the purpose of this discussion, it is sufficient to consider the onedimensional problem of calculating the advective flux at a tracer cell face. The advective flux at a tracer cell face can be written as
where F_{f} is the advective flux at a tracer cell face, u_{f} is the normal velocity at a tracer cell face, c_{u} corresponds to the tracer concentration in the upwind direction, c_{d} corresponds to the tracer concentration in the downwind direction, and Ψ(r) is the limiter function where r is a measure of the monotonicity of c. The monotonicity in this method is defined as
where c_{u−1} is the tracer concentration one cell further in the upstream direction. When r<0, the tracer concentration has a local maximum or minimum centered around c_{u}. For the ULTIMATE flux limiter, Ψ(r)=0 when r<0, in which case the flux evaluates to the upwind flux.
The impulse functions used in this study all have a local maximum centered at the impulse locations. At the impulse locations r<0, and therefore the diagnosed flux is the more diffusive upwind flux. This behavior will be apparent and quantified when comparing the online and offline simulations in Sect. 3.2.
2.5 Time integration
Forward integration of the offline tracer model uses a backward Euler method for vertical mixing and can use forward Euler, Adams–Bashforth second order, or Adams–Bashforth thirdorder methods for transport. As in Bardin et al. (2014), we forward step an equation for the volume anomaly using a forward Euler method. Volume anomalies arise due to divergence in the transport field at the uppermost z level that are associated with fluctuations of the free surface.
In general, the time integration scheme can be written as
where v is the volume anomaly, i is a constant vector whose elements are all set to one, c^{*} depends on the time integration scheme that is used (Table 1), and 𝕍^{n+1} is a diagonal matrix whose diagonal elements are the volume anomalies.
In ocean models, advection and horizontal diffusion operators have a compact stencil, enabling the use of the IRF approach described in the previous section to “capture” these advection and diffusion operators. However, vertical diffusion is usually treated differently, by solving a tridiagonal system that touches the entire water column. The reason is that high values of vertical diffusivity are applied where the water column is unstable, and these values easily render any explicit scheme unstable. Consequently, as the region of influence of the vertical diffusion operator is the entire water column, the IRF approach would demand a separate IRF field for each vertical level of the model grid, which is prohibitive for finely resolved grids. Instead, FEOTS treats the vertical solve similarly as the parent model, so rather than IRFs, the vertical diffusivities are diagnosed, saved, and used to recreate the vertical diffusion operators offline.
Forwardstepping FEOTS requires inverting a tridiagonal system of equations, given by Eq. (12), for the tracer concentration in order to incorporate vertical mixing. To solve this system, we use the preconditioned conjugate gradient algorithm (Shewchuk, 1994) with a diagonal preconditioner. The initial solution guess for the vertical mixing solver is set as the tracer concentration that is predicted without vertical mixing,
For the results presented in this paper, we use the thirdorder Adams–Bashforth time integrator, and the conjugate gradient solver is stopped when the residual magnitude, relative to the initial solution guess magnitude, is less than 10^{−6}. We use a 15 min time step, and a typical maximum CFL value, obtained by eigenvalue analysis of the transport operators, is O(0.1).
2.6 Constant preservation test problem
We now show that constant preservation is expected in the discrete system. This provides the foundation for a test case that we use to verify the FEOTS implementation. From Eqs. (11) and (12), we start by assuming that there is no lateral or vertical diffusion, and the initial tracer field is a constant value of 1 (c^{0}=i). The discrete system at the initial time step is then
Note that we can write these equations using indicial notation, which makes it easier to obtain the result for c^{1}:
Substituting Eq. (15a) into (15b) gives
The only solution to Eq. (16) is ${c}_{j}^{\mathrm{1}}=\mathrm{1}$ for each j, implying that the discrete solution remains a constant. If we progress to the next time step using a forward Euler method, the same result is obtained. If instead we switch to a secondorder Adams–Bashforth method for the next time step, we have
Using this for c^{*} results in
which again implies that c^{2}=i. This result also holds for the thirdorder Adams–Bashforth method in this case, since
when ${\mathit{c}}^{\mathrm{2}}={\mathit{c}}^{\mathrm{1}}={\mathit{c}}^{\mathrm{0}}=\mathit{i}$.
When vertical mixing is introduced, the solution obtained in Eq. (16) serves as an initial guess to the conjugate gradient solver. If the initial guess for the solver is a constant tracer field,
This result in Eq. (20) is due to the use of a finite volume discretization for the vertical diffusion with noflux conditions at the ocean surface and at the seafloor. In exact arithmetic, we would therefore expect the tracer field to remain a constant.
This analysis shows that a constant tracer field remains a constant tracer field under the discretizations employed in FEOTS. Finite precision arithmetic, however, can produce slight deviations from a constant field, and the results presented in Sect. 3.1 characterize the behavior of roundoff errors for a constant tracer field, both with and without vertical mixing.
2.7 The Argentine Basin test problem
As discussed in the introduction, the Argentine Basin is an ideal region for testing a tracer transport capability in an eddying ocean model. In this paper we compare online and offline simulations of dye tracers that are initialized at the boundaries of the Argentine Basin, here chosen as (52.18^{∘} S, 28.06^{∘} S) × (70.25^{∘} W, 24.90^{∘} W). In all simulations, the only source of tracers comes from the model boundary conditions applied along the southern, eastern, and northern boundaries.
We want to distinguish between water masses that originate from each of the domain boundaries and above and below 1000 m depth. This is accomplished by simulating six passive tracer fields D_{i} with boundary conditions:
where H(x) is the Heaviside step function. With this configuration, tracers D_{1}, D_{3}, and D_{5} are released in the upper 1000 m on the southern, eastern, and northern boundaries, respectively; while tracers D_{2}, D_{4}, and D_{6} are released below 1000 m. Note that maximum mixedlayer depths in the Argentine Basin in winter are around 500 m in this model, so deep convection should not play a role in the transport of these tracers across the 1000 m depth horizon.
In this section, we present results of an offline regional study focused on the Argentine Basin. We first present verification of the constant preservation property (Sect. 3.1). Then we present a comparison of the tracer simulation results in the parent model (online) and the offline model using 5 d averaged transport operators (Sect. 3.2), followed by a comparison between 1 and 5 d averaged transport operators (Sect. 3.3). Finally, we discuss the computational performance of the code (Sect. 3.4) on the systems where our simulations were conducted.
3.1 Constant preservation
The parent model uses a finite volume discretization that guarantees the preservation of constant tracer fields. To verify that FEOTS accurately diagnoses transport operators that are representative of the parent model, our first simulation involves verifying that a constant tracer field remains constant under the action of the diagnosed transport operators.
In all of our simulations, we have opted to use single precision arithmetic and have enabled aggressive compiler optimizations (compiler option Ofast with GCC 9.2.0). These choices were made to minimize data storage costs for the transport operators and postprocessing output and to optimize the timetosolution for the offline simulations. Although analytically we expect that the FEOTS algorithm should preserve constant tracer fields, errors from floating point arithmetic are expected to be the main source of constant preservation errors.
To better understand the spatial distribution and the relative impact of roundoff errors, we configure a simulation where the initial tracer field is set to c_{0}=1, and there is no external source or sink of tracer. Since ∇c_{0}=0 and $\mathrm{\nabla}\cdot \mathit{u}=\mathrm{0}$ under the discretizations used in the parent model, we expect that c_{t}=0 and $c={c}_{\mathrm{0}}=\mathrm{1}$ for all time, as demonstrated in Sect. 2.6.
Figure 1 shows the maximum volume anomaly in the domain as a function of depth after 5 d and 5 years of integration. Analytically, the volume anomaly is expected to be zero in all cells except the topmost layer. At the surface layer, fluctuations in the freesurface height are associated with nonzero fluid divergences that contribute to changes in the fluid volume. Beneath the surface layer, the fluid velocity field is expected to be divergence free. In general, larger errors in the volume anomaly are observed above 1000 m depth. After 5 d, errors in the deep ocean are O(10^{−5}), and after 5 years, the deep ocean volume anomaly errors have grown by an order of magnitude to O(10^{−4}). Larger errors are observed above 1000 m, reaching O(10^{−3}) after 5 years. Note that the volume anomaly field is identical for all choices of the time integrator for the dye tracer and is independent of vertical mixing.
Errors in the volume anomaly lead to spurious values for predicted tracer concentrations. For this simulation, any deviation of the tracer concentration from its initial uniform value is erroneous. Figure 2 shows the max error in the dye tracer as a function of depth after 5 d and 5 years of integration, with and without mixing. After 5 years of integration, the maximum relative error with mixing is about 0.05 %, and without mixing it is about 0.01 %. At depth, the errors in the tracer are comparable with and without mixing. However, above 1000 m, particularly in the mixed layer, the inclusion of mixing results in an accumulation of roundoff errors.
3.2 The Argentine Basin test case: offline vs. online comparison
In addition to quantifying errors for the constant tracer scenario, a practical concern is in the comparison between the online and offline tracer simulations. Do the tracer distributions simulated by FEOTS using transport operators diagnosed from the parent model E3SMv0HiLAT03 (offline) faithfully represent the tracer distributions simulated by the parent model itself (online)? The time averaging of the transport operators will introduce differences between the online and offline simulations, as will the fluxcorrected advection operator. Our aim is to quantify and qualitatively describe the differences between a tracer simulation conducted directly in the online parent model and the offline model. We do this for our example problem, namely determining the source waters of the Zapiola Anticyclone.
Figure 3 compares the tracer fields obtained with the online and offline method, at the end of the 5year analysis period. The dye tracers shown are those that are released in the upper 1000 m at each of the three domain boundaries. Visual inspection shows that the online simulation has sharper gradients compared to the offline simulation, even though the same dynamical features are visible. This is a result of the nonlinear advection operator that is more diffusive when applied to impulse functions than to a smooth tracer field, as explained in Sect. 2.4.
The advective errors also affect integrated quantities of dye tracers within the Zapiola Anticyclone. Figure 4 shows the vertical distribution of tracer concentrations averaged over the Zapiola Anticyclone, while Fig. 5 shows the total tracer stock.
Dye tracer D_{1}, sourced at the southern boundary in the upper 1000 m, reaches the ZA after about 45 d in the online simulation (Fig. 4, upper row). Concentrations in the upper 1000 m rise steadily in the first 2 years, then become mostly steady, indicating saturation is reached (Fig. 5a, black dashed). Just below the surface, concentrations in the final 3 years average about 0.76, suggesting that 76 % of the surface waters within the ZA may be derived from the Southern Ocean through the Malvinas Current. At 984 m this fraction is 0.45. Tracer stock below 1000 m increases slowly and contains 17 % of the column stock after 5 years (Fig. 5a, black dotted). The offline simulation reproduces this behavior qualitatively but with significant quantitative differences. D_{1} reaches the ZA about 5 d earlier and increases faster than the online simulation, but upperlayer stock plateaus at a slightly lower level (Fig. 5a, gray dashed). Tracer fraction reaches 0.74 just below the surface, close to the online simulation, but the saturation value of 0.37 at 1000 m depth is significantly lower. Figure 4 (upper row) shows that this is due to a significant vertical redistribution of tracers by vertical diffusion that depletes tracers from the upper 1000 m and increases concentrations below. Indeed, the inventory below 1000 m accounts for 30 % of the column stock after 5 years.
The story is similar for dye tracer D_{2}, which is sourced at the southern boundary below 1000 m. It arrives in the ZA after 87 d in the online simulation. The stock below 1000 m (Fig. 5b, black dotted) rises quickly through the first 2 years, followed by a more gradual (linear) rise after that. Tracer concentrations at 1049 m level out at 0.20 after 3 years. The total contribution of Southern Ocean waters (D_{1}+D_{2}) is about 59 % at 1000 m depth. The offline simulation reproduces this behavior quite well but displays higher inventories in the upper 1000 m that increase the overall column stock by 9 % after 5 years. The vertical profiles again clearly show the impact of vertical diffusion that depletes tracers in the 1000–3000 m depth range and increases concentrations below and above (Fig. 4, second row). The D_{2} concentration at 1049 m depth levels out at 0.21, with a total Southern Ocean contribution at 1000 m of 56 %.
The nextlargest contribution to the ZA tracer inventory is coming from the north through dye tracers D_{5} (upper 1000 m) and D_{6} (below 1000 m). It takes about 110 d for D_{5} to arrive at the ZA, and the surface concentration saturates at about 0.18 after 550 d (as does the upperlayer tracer stock; Fig. 5e, black dashed, overlain by solid). This suggests that the Brazil Current may contribute about 20 % of the surface waters in the ZA. At 984 m, however, this fraction is still only 0.015 after 5 years, and rising, probably reflecting the strongly sheared character of the Brazil Current and the long transit time from the northern domain boundary to the ZA. Notable concentrations of D_{6} reach the ZA after 2 years, but trace quantities already arrive after about 270 d, having been mixed upward into the upper layer and transported southward in the Brazil Current. The offline simulations display qualitatively similar behavior, but D_{5} inventories are significantly higher (+40 %), with again significant sequestration below 1000 m. D_{6} stock is slightly lower (−7 %) than in the online simulations, despite higher values above 1000 m.
Dye tracers D_{3} and D_{4}, released at the eastern domain boundary, take much longer to reach the ZA due to very low westward flow velocities in the interior part of the basin. Offline stock of D_{3} in the upper 1000 m is about 50 % smaller than in the online simulation, a deficiency that can only be partly explained by sequestration below 1000 m. D_{4} stock is 65 % too high, with enhanced stock both below and above 1000 m depth.
Based on the propagation speed of the diffusion front of D_{1} in the offline simulation (Fig. 4a, center column), we can make a rough estimate of the artificial vertical diffusivity that is introduced by the advection issue. We are able to model the depth of the diffusion front below ${z}_{\mathrm{0}}=\mathrm{1000}$ m and after t_{0}=132 d as $z={z}_{\mathrm{0}}\sqrt{\mathrm{4}D(t{t}_{\mathrm{0}})}$, when $D=\mathrm{1.74}\times {\mathrm{10}}^{\mathrm{2}}$ m^{2} s^{−1}. This is 2 to 3 orders of magnitude larger than typical values for background diffusivity used in ocean models.
3.3 The Argentine Basin test case: the role of temporal averaging
The parent model is capable of producing velocity fields that have a wide range of scales of spatial and temporal variability. The shortest temporal periods are on the order of a few time steps, and the longest period is the duration of the simulation. In general, higherresolution models introduce more variability on shorter length and time scales, and some consideration is needed when selecting an averaging period for the transport operator diagnosis. For storage reasons, it is not practical to store snapshots of the transport operators at every time step. Conversely, representing the ocean transport with long time averages may exclude the effects of important variability. The choice of the timeaveraging period for the transport operators can impact the evolution of tracers calculated in FEOTS, and an appropriate balance of practicality and accuracy should be struck.
To our knowledge, there is currently no definitive guidance on choosing the averaging period for velocity fields in an offline tracer simulation; the averaging period should be chosen so that variability in the underlying advection field is well resolved. Here, the impact of temporal averaging is investigated by comparing simulations with 1 and 5 d averaged operators. To that end, we ran another online simulation for which we saved 1 d averaged IRFs. This simulation was run for 105 d, producing 105 operator sets. The offline simulations were run for 365 d hence cycling through the operator set almost 3.5 times.
The online simulations that produced the 1 and 5 d averaged operators were not bitforbit identical, as we inadvertently specified different processor counts. A positive consequence of this oversight is that we have two different realizations of the same chaotic system, allowing us to assess how the tracer advection error compares to trajectory divergence, in terms of their impact on tracer stock. Figure 6 shows that averaging 1 d averaged operators to 5 d averages has a very small impact on the tracer stock (light gray); the impact on spatial distributions is similarly small (not shown). This shows that 5 d averaged operators are sufficient to accurately reproduce tracer distributions. Comparison of the black and dark gray lines shows that the impact of trajectory divergence on tracer stock (dark gray) is of similar magnitude as the impact of the advection error (black), at least for the first year of simulations.
3.4 Computational performance
One goal of FEOTS is to perform tracer calculations at a lower computational cost than the parent model. Additionally, FEOTS allows researchers to take advantage of transport operators produced by stateoftheart climate simulations to conduct regional offline simulations. This provides flexibility in studying ocean transport phenomena and increases the value of online produced model data while considerably reducing the computational expense for researchers solely interested in studies involving passive tracers. Here we evaluate the computational performance of a regional FEOTS configuration and compare it with the global parent model.
The total cost of using FEOTS is associated with the following steps:

impulse functions are generated from the model grid,

an online simulation is done with POP that reads in the impulse functions and outputs IRFs averaged over the desired averaging interval (e.g., 1, 5 d),

the diagnosed IRFs are translated from gridded output to a sparse matrix format, and

offline passive tracer simulations are run.
The first three steps are onetime costs that are necessary to generate the transport operator database. In our experience, impulse function generation introduces a negligible cost, requiring only a few minutes to run in serial. Simulation of the passive tracers with the parent model to generate the impulse response functions requires about a factor of six more CPU hours than when running without tracers. This is the most significant upfront cost in generating the transport operator database. Diagnosis of the sparse matrices introduces a small cost; for the parent model presented here, about 15 min of wall time on a single core is needed per transport operator. The expense of the offline passive tracer simulation depends on the specific use case. Below, we provide an example based on the Argentine Basin test problem and a simple global simulation.
Table 2 summarizes a comparison of the computational expense and compute platform size requirements between POP and FEOTS. We ran the 0.3^{∘} ocean/sea ice configuration of E3SMHiLAT on the Los Alamos National Laboratory's Institutional Computing clusters. A typical simulation with six dye tracers costs 9020 CPU h per simulated year, using 2432 cores on 76 nodes. The throughput is 6.5 simulated years per wallclock day. The simulation with 53 IRFs (in addition to the six dye tracers) typically costs 45 000 CPU h, with a throughput of 1.60 simulated years per day when using 2912 cores.
In contrast, a 1year offline FEOTS simulation of the Argentine Basin problem with six tracers takes roughly 47 CPU h and a throughput of 3 simulated years per day on six cores. Table 3 shows the five most expensive routines in FEOTS and the percentage of the total wall time spent executing those routines. The iterative treatment of vertical mixing is most expensive; a 1year simulation without vertical mixing costs 6.8 CPU h and provides a throughput of 20 simulated years per day on six cores. All of these routines execute sparse matrix–vector multiplication in order to compute advective and diffusive tendencies, suggesting future improvements to sparse matrix–vector multiplication in FEOTS would be beneficial.
Currently, offline global simulations at eddyresolving resolutions are slow and require a large amount of memory per core. For our global model, about 15 GB of memory is needed for the offline simulation at single precision and 30 GB at double precision. At the time of this study, our focus was on proving the transient simulation capabilities, which could be done at a low computational cost in regional simulations; we did not obtain direct measurements of FEOTS' runtime for global transient simulations. However, we are interested in understanding whether this methodology is potentially computationally competitive in comparison to online tracer simulations in the parent model.
To develop an estimate for a global offline simulation with FEOTS, we assume that the CPU hours scale linearly with the number of grid cells. This is reasonable, since the majority of FEOTS' runtime is spent in sparse matrix–vector multiplication, which scales linearly with the number of rows. The Argentine Basin simulation has 1.02×10^{6} grid cells, whereas a global configuration has about 45 times more grid cells. Thus, an offline simulation with FEOTS is expected to cost about 2115 CPU h per simulation year in the current version, 77 % less CPU hours than the parent model. However, the expected throughput for an offline global simulation is estimated as 1 model year for every 15 d. Although FEOTS is currently expected to take longer to produce global transient simulations, it is probable that exposing parallelism in FEOTS will provide comparable runtimes to the parent model at a reduced computational expense. Though this is not a definitive result, it is a reasonable estimate that motivates future work in exposing parallelism in FEOTS.
Based on this work, we posit that offline tracer simulations provide a viable modeling capability at a significantly reduced computational expense compared to online models. With the reduction in computational expense, however, we have shown that the offline simulations can produce tracer distributions consistent with a more diffusive ocean circulation than online models. A potential approach to mitigate this problem, when using nonlinear fluxlimiting advection schemes, is to replace the impulse fields with smoother basis functions, as discussed by Khatiwala et al. (2005).
When using impulse fields with sharp discontinuities, the thirdorder fluxlimited Lax–Wendroff advection scheme reduces effectively to a firstorder upwind scheme. By leveraging smoother basis functions for the impulse fields with this scheme, a less diffusive advection operator can be diagnosed. Using a smoother basis function, however, is expected to introduce additional complications for diagnosing transport operators. In general, the process of diagnosing the transport operators can be thought of as a matrix projection problem,
where 𝔸 is the transport operator we want to diagnose, 𝔽 is a matrix whose columns are the impulse fields, and ℝ is a matrix whose columns are the impulse response fields. The advection operator is obtained by multiplying Eq. (22) by the right inverse of 𝔽:
Ideally, the basis function we choose should be smooth enough to retain higherorder terms in the thirdorder fluxlimited Lax–Wendroff advection scheme. For computational purposes, the basis functions would ideally be mutually orthonormal so that the inverse is easy to calculate:
This approach, along with experimentation with other advection schemes, is planned for future work.
With regard to performance, FEOTS can provide offline tracer modeling capabilities at a significantly reduced computational expense (CPU hours) and with practical runtimes for regional simulations. Projected computational resource requirements (CPU hours) were shown to be about 77 % less than the online parent model. Estimated runtimes for global simulations, however, indicate that parallelism in FEOTS needs to be exposed before it is practical for this use case. Addressing the slow runtime for global offline simulations is critical for working towards a framework that allows for the calculation of steadystate solutions within a practical amount of time. Currently, our plan for reducing runtime is to leverage opensource toolkits, like PSBLAS (Filippone and Colajanni, 2000), to handle sparse matrix operations in parallel.
In this paper we introduced the Fast Equilibration of Ocean Tracers Software (FEOTS), which is an endtoend set of tools to efficiently calculate tracer distributions on a global or regional subdomain. Key features currently implemented in FEOTS and discussed in this paper include impulse field generation through graph coloring, impulse response function translation into a sparse matrix representation of transport operators (advection), regional subdomain transport operator generation from global transport operators, and offline forward simulation with the diagnosed operators. These are key components for calculating equilibrium tracer fields, which we aim to implement in future releases. In our example problem we diagnose transport operators from a global eddypermitting configuration of the POP ocean model and calculate the distribution of passive dye tracers in the Argentine Basin for a 5year period. The offline simulations are shown to produce more diffuse tracer distributions than the online parent model simulations. This demonstration shows the feasibility of this approach, while at the same time highlighting the challenges of the impulse response functions approach.
The current version of FEOTS (v0.0.0) is available from the project website: https://github.com/FluidNumerics/FEOTS/ under the 3Clause BSD License. The exact version of the model used to produce the results used in this paper is archived on Zenodo (https://doi.org/10.5281/zenodo.5576912, Schoonover et al., 2021). The input data and scripts to run the model and produce the plots for all the simulations presented in this paper are also archived on Zenodo (https://doi.org/10.5281/zenodo.6250938, Schoonover et al., 2022). A codelab tutorial to walk through running the Argentine Basin simulations is available online at https://fluidnumerics.github.io/FEOTS/codelabs/feotsongooglecloud/#0 (last access: 22 May 2023).
JS and WW wrote the manuscript. JS conceived, developed and wrote the FEOTS code, and performed the error analysis; WW performed online and offline tracer simulations and analysis; and JZ contributed to the development and testing of FEOTS.
The contact author has declared that none of the authors has any competing interests.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
We are grateful to Keith Lindsay (NCAR) as well as Ann Bardin and François Primeau (UCI) for sharing code and for useful discussions. This research was supported by the Los Alamos National Laboratory (LANL) through its Center for Space and Earth Science (CSES). CSES is funded by LANL's Laboratory Directed Research and Development (LDRD) program under project number 20210528CR. Joseph Schoonover and Jiaxu Zhang acknowledge support from LANL's Center for Nonlinear Studies (CNLS). Computations were performed on LANL's Institutional Computing platforms and on the Darwin cluster. This research was supported by the Regional and Global Model Analysis (RGMA) component of the Earth and Environmental System Modeling (EESM) program of the U.S. Department of Energy's Office of Science, as a contribution to the HiLATRASM project. We thank Ann Bardin and an anonymous reviewer for constructive reviews of this paper.
This research was supported by the Los Alamos National Laboratory (LANL) through its Center for Space and Earth Science (CSES). CSES is funded by LANL's Laboratory Directed Research and Development (LDRD) program under project number 20210528CR.
This paper was edited by Olivier Marti and reviewed by Ann Bardin and one anonymous referee.
Bardin, A., Primeau, F., and Lindsay, K.: An offline implicit solver for simulating prebomb radiocarbon, Ocean Model., 73, 45–58, 2014. a, b, c, d, e, f
Chamberlain, M. A., Matear, R. J., Holzer, M., Bi, D., and Marsland, S. J.: Transport matrices from standard oceanmodel output and quantifying circulation response to climate change, Ocean Model., 135, 1–13, 2019. a
de Miranda, A. P., Barnier, B., and Dewar, W. K.: On the dynamics of the Zapiola Anticyclone, J. Geophys. Res.Oceans, 104, 21137–21149, 1999. a
Dewar, W. K.: Topography and barotropic transport control by bottom friction, J. Mar. Res., 56, 295–328, 1998. a
Dukhovskoy, D. S., Myers, P. G., Platov, G., Timmermans, M.L., Curry, B., Proshutinsky, A., Bamber, J. L., Chassignet, E., Hu, X., Lee, C. M., and Somavilla, R.: Greenland freshwater pathways in the subArctic Seas from model experiments with passive tracers, J. Geophys. Res.Oceans, 121, 877–907, 2016. a
Filippone, S. and Colajanni, M.: PSBLAS: A Library for Parallel Linear Algebra Computation on Sparse Matrices, ACM T. Math. Software, 26, 527–550, https://doi.org/10.1145/365723.365732, 2000. a
Fu, L.L. and Smith, R. D.: Global ocean circulation from satellite altimetry and highresolution computer simulation, B. Am. Meteorol. Soc., 77, 2625–2636, 1996. a, b
Garzoli, S. L.: Geostrophic velocity and transport variability in the BrazilMalvinas Confluence, DeepSea Res. Pt. I, 40, 1379–1403, 1993. a
Griffies, S. M., Winton, M., Samuels, B., Danabasoglu, G., Yeager, S., Marsland, S., Drange, H., and Bentsen, M.: Datasets and protocol for the CLIVAR WGOMD Coordinated Oceansea ice Reference Experiments (COREs), WCRP Report, 21, 1–21, 2012. a
Gu, S., Liu, Z., Jahn, A., Rempfer, J., Zhang, J., and Joos, F.: Modeling neodymium isotopes in the ocean component of the community earth system model (CESM1), J. Adv. Model. Earth Sy., 11, 624–640, 2019. a
Hallberg, R.: Using a resolution function to regulate parameterizations of oceanic mesoscale eddy effects, Ocean Model., 72, 92–103, 2013. a
Hecht, M., Veneziani, M., Weijer, W., Kravitz, B., Burrows, S., Comeau, D., Hunke, E., Jeffery, N., UrregoBlanco, J., Wang, H., Wang, S., Zhang, J., Bailey, D., Mills, C., Rasch, P., and Urban, N.: E3SMv0HiLAT: A modified climate system model targeted for the study of highlatitude processes, J. Adv. Model. Earth Sy., 11, 2814–2843, 2019. a
Hecht, M. W., Hunke, E., Maltrud, M., Petersen, M. R., and Wingate, B. A.: Lateral mixing in the eddying regime and a new broadranging formulation, Eddy resolving ocean models, Geophys. Monogr., 177, 339–352, 2008. a
Hundsdorfer, W. and Trompert, R.: Method of lines and direct discretization: a comparison for linear advection, Appl. Numer. Math., 13, 469–490, https://doi.org/10.1016/01689274(94)900094, 1994. a
Jahn, A., Lindsay, K., Giraud, X., Gruber, N., OttoBliesner, B. L., Liu, Z., and Brady, E. C.: Carbon isotopes in the ocean model of the Community Earth System Model (CESM1), Geosci. Model Dev., 8, 2419–2434, https://doi.org/10.5194/gmd824192015, 2015. a
Jullion, L., Heywood, K. J., Naveira Garabato, A. C., and Stevens, D. P.: Circulation and water mass modification in the Brazil–Malvinas Confluence, J. Phys. Oceanogr., 40, 845–864, 2010. a
Khatiwala, S.: A computational framework for simulation of biogeochemical tracers in the ocean, Global Biogeochem. Cy., 21, GB3001, https://doi.org/10.1029/2007GB002923, 2007. a
Khatiwala, S., Visbeck, M., and Cane, M. A.: Accelerated simulation of passive tracers in ocean circulation models, Ocean Model., 9, 51–69, 2005. a, b, c, d
Khatiwala, S., Primeau, F., and Hall, T.: Reconstruction of the history of anthropogenic CO_{2} concentrations in the ocean, Nature, 462, 346–349, 2009. a
Kvale, K. F., Khatiwala, S., Dietze, H., Kriest, I., and Oschlies, A.: Evaluation of the transport matrix method for simulation of ocean biogeochemical tracers, Geosci. Model Dev., 10, 2425–2445, https://doi.org/10.5194/gmd1024252017, 2017. a
Leonard, B.: The ULTIMATE conservative difference scheme applied to unsteady onedimensional advection, Comput. Method. Appl. M., 88, 17–74, https://doi.org/10.1016/00457825(91)90232U, 1991. a
Lozier, M. S.: Deconstructing the conveyor belt, Science, 328, 1507–1511, 2010. a
Missiaen, L., Menviel, L. C., Meissner, K. J., Roche, D. M., Dutay, J.C., Bouttes, N., Lhardy, F., Quiquet, A., Pichat, S., and Waelbroeck, C.: Modelling the impact of biogenic particle flux intensity and composition on sedimentary Pa/Th, Quaternary Sci. Rev., 240, 106394, https://doi.org/10.1016/j.quascirev.2020.106394, 2020. a
Mountford, A. and Morales Maqueda, M.: Eulerian modeling of the threedimensional distribution of seven popular microplastic types in the global ocean, J. Geophys. Res.Oceans, 124, 8558–8573, 2019. a
Primeau, F.: Characterizing transport between the surface mixed layer and the ocean interior with a forward and adjoint global ocean transport model, J. Phys. Oceanogr., 35, 545–564, 2005. a
Ringler, T., Petersen, M., Higdon, R. L., Jacobsen, D., Jones, P. W., and Maltrud, M.: A multiresolution approach to global ocean modeling, Ocean Model., 69, 211–232, 2013. a
Sarmiento, J. L., Orr, J. C., and Siegenthaler, U.: A perturbation simulation of CO_{2} uptake in an ocean general circulation model, J. Geophys. Res.Oceans, 97, 3621–3645, 1992. a
Saunders, P. M. and King, B. A.: Bottom currents derived from a shipborne ADCP on WOCE cruise A11 in the South Atlantic, J. Phys. Oceanogr., 25, 329–347, 1995. a
Schoonover, J., Zhang, J., and Weijer, W.: FluidNumerics/FEOTS: v0.0.0 (v0.0.0), Zenodo [code], https://doi.org/10.5281/zenodo.5576912, 2021. a
Schoonover, J., Weijer, W., and Zhang, J.: FEOTS Argentine Basin Transport Operators (5day Average), Zenodo [data set], https://doi.org/10.5281/zenodo.6250938, 2022. a
Séférian, R., Berthet, S., Yool, A., Palmieri, J., Bopp, L., Tagliabue, A., Kwiatkowski, L., Aumont, O., Christian, J., Dunne, J., and Gehlen, M.: Tracking improvement in simulated marine biogeochemistry between CMIP5 and CMIP6, Current Climate Change Reports, 1–25, https://doi.org/10.1007/s40641020001600, 2020. a
Shewchuk, J. R.: An Introduction to the Conjugate Gradient Method Without the Agonizing Pain, Tech. rep., School of Computer Science, Carnegie Mellon University, http://www.cs.cmu.edu/~quakepapers/painlessconjugategradient.pdf (last access: 17 May 2023), 1994. a
Smith, R., Jones, P., Briegleb, B., Bryan, F., Danabasoglu, G., Dennis, J., Dukowicz, J., Eden, C., FoxKemper, B., Gent, P., Hecht, M., Jayne, S., Jochum, M., Large, W., Lindsay, K., Maltrud, M., Norton, N., Peacock, S., Vertenstein, M., and Yeager, S.: The Parallel Ocean Program (POP) reference manual: Ocean component of the Community Climate System Model (CCSM), Tech. rep., Los Alamos National Laboratory, http://nldr.library.ucar.edu/repository/collections/OSGC000000000954 (last access: 17 May 2023), 2010. a, b, c, d
Wang, S., Moore, J. K., Primeau, F. W., and Khatiwala, S.: Simulation of anthropogenic CO_{2} uptake in the CCSM3.1 ocean circulationbiogeochemical model: comparison with databased estimates, Biogeosciences, 9, 1321–1336, https://doi.org/10.5194/bg913212012, 2012. a
Weijer, W., Maltrud, M. E., Homoky, W. B., Polzin, K. L., and Maas, L. R.: Eddydriven sediment transport in the Argentine Basin: Is the height of the Zapiola Rise hydrodynamically controlled?, J. Geophys. Res.Oceans, 120, 2096–2111, 2015. a
Weijer, W., Barthel, A., Veneziani, M., and Steiner, H.: The Zapiola Anticyclone: A Lagrangian study of its kinematics in an eddypermitting ocean model, DeepSea Res. Pt. I, 164, 103308, https://doi.org/10.1016/j.dsr.2020.103308, 2020. a, b, c
Zanna, L., Khatiwala, S., Gregory, J. M., Ison, J., and Heimbach, P.: Global reconstruction of historical ocean heat storage and transport, P. Natl. Acad. Sci. USA, 116, 1126–1131, 2019. a
Zhang, J., Liu, Z., Brady, E. C., Oppo, D. W., Clark, P. U., Jahn, A., Marcott, S. A., and Lindsay, K.: Asynchronous warming and δ^{18}O evolution of deep Atlantic water masses during the last deglaciation, P. Natl. Acad. Sci. USA, 114, 11075–11080, 2017. a
Zhang, J., Weijer, W., Maltrud, M. E., Veneziani, C., Jeffery, N., Hunke, E. C., Urrego Blanco, J. R., and Wolfe, J. D.: An eddypermitting oceansea ice general circulation model (E3SMv0HiLAT03): Description and evaluation, Tech. rep., Los Alamos National Lab.(LANL), Los Alamos, NM (United States), https://doi.org/10.2172/1542803, 2019. a, b
Zhang, J., Weijer, W., Steele, M., Cheng, W., Verma, T., and Veneziani, M.: Labrador Sea freshening linked to Beaufort Gyre freshwater release, Nat. Commun., 12, 1229, https://doi.org/10.1038/s41467021214703, 2021. a