Approximating Isoneutral Ocean Transport via the Temporal Residual Mean

Ocean volume and tracer transports are commonly computed on density surfaces because doing so approximates the semi-Lagrangian mean advective transport. The resulting density-averaged transport can be related approximately to Eulerian-averaged quantities via the Temporal Residual Mean (TRM), valid in the limit of small isopycnal height fluctuations. This article builds on a formulation of the TRM for volume fluxes within Neutral Density surfaces, (the “NDTRM”), selected because Neutral Density surfaces are constructed to be as neutral as possible while still forming well-defined surfaces. This article derives a TRM, referred to as the “Neutral TRM” (NTRM), that approximates volume fluxes within surfaces whose vertical fluctuations are defined directly by the neutral relation. The purpose of the NTRM is to more closely approximate the semi-Lagrangian mean transport than the NDTRM, because the latter introduces errors associated with differences between the instantaneous state of the modeled/observed ocean and the reference climatology used to assign the Neutral Density variable. It is shown that the NDTRM collapses to the NTRM in the limiting case of a Neutral Density variable defined with reference to the Eulerian-mean salinity, potential temperature and pressure, rather than an external reference climatology, and therefore that the NTRM approximately advects this density variable. This prediction is verified directly using output from an idealized eddy-resolving numerical model. The NTRM therefore offers an efficient and accurate estimate of modeled semi-Lagrangian mean transports without reference to an external reference climatology, but requires that a Neutral Density variable be computed once from the model’s time-mean state in order to estimate isopycnal and diapycnal components of the transport.


Introduction
Accurate quantification of the ocean's meridional overturning circulation (MOC) is required to infer oceanic transport of heat and other tracers around the globe [1,2]. It has been common practice for decades to quantify the MOC in some form of density coordinate system (e.g., [3][4][5][6][7][8]). Because flow in the ocean interior is quasi-adiabatic, this approach approximately corresponds to taking a semi-Lagrangian mean, which in turn approximates the full Lagrangian-mean transport [9]. In comparison to an MOC based on mass fluxes averaged at constant depth, the semi-Lagrangian mean transport velocity contains an additional "eddy" transport, essentially a generalized Stokes drift [10]. This eddy transport is particularly important in the Southern Ocean [11][12][13], and therefore forms a key component of recent conceptual models of the global overturning circulation [14][15][16][17][18][19].
The separation of the semi-Lagrangian mean transport into mean and eddy components can be made explicit via the Temporal Residual Mean (TRM) formulation, which is valid in the asymptotic limit of small isopycnal height fluctuations [9,20,21]. In this article we use "TRM" to refer to the TRM-II [21], but it is straightforward to extend our results to the TRM-I; the distinction between the two is discussed by [22]. The TRM relates the eddy component of the transport to correlations between the lateral velocity (u) and density (γ) fields, and thereby allows an estimate of the semi-Lagrangian mean transport streamfunction to be computed from Eulerian-averaged quantities. The resulting TRM streamfunction Ψ advects the density variable γ in the following sense: in the absence of non-conservative sources and sinks of γ, the density γ is materially conserved following the transport velocity defined by Ψ, up to a consistent asymptotic order in the amplitude of the isopycnal height fluctuations [22].
This article builds upon the Neutral Density Temporal Residual Mean (NDTRM) formulation ( [23], [hereafter ST15]). In ST15's formulation the Neutral Density variable of [24], referred to henceforth as γ JM97 , was the advected density variable. However, the transport streamfunction Ψ was calculated from correlations between the lateral velocity, salinity, potential temperature and pressure. This approach bears the following advantages: (i) The NDTRM can be applied in regions where potential density referenced to any level becomes vertically non-monotonic somewhere in the water column (e.g., parts of the high-latitude oceans). (ii) The NDTRM circumvents the computational expense of calculating Neutral Density globally at every model time step, as would be required to compute the semi-Lagrangian mean transport velocity, either exactly or via the TRM. (iii) More fundamentally, isopycnals of Neutral Density are more closely aligned with the local neutral tangent plane, in a global sense, than any potential density variable [24]. ST15 compared various formulations of the TRM, and verified that higher-order approximations to the NDTRM most closely approximated the volume fluxes computed exactly within γ JM97 surfaces from numerical model output. The NDTRM therefore offers a computationally efficient estimate of isopycnal volume fluxes within surfaces of a stably stratified density variable that has been used widely in quantifications of the ocean's global overturning circulation (e.g., [5,25]).
In this article we pursue the related but distinct goal of approximating the semi-Lagrangian mean transport using vertical fluctuations defined directly by the neutral relationship [26], which we refer to as local Neutral Surfaces. In pursing this goal we make the assumption that it is these local Neutral Surfaces that should heave vertically under perfectly adiabatic motions [27,28], and therefore that they are best suited for calculation of the semi-Lagrangian mean transport [9]. The validity of these assumptions is discussed further in Section 5. In Section 2 we summarize relevant background literature, specifically the formulation of the TRM [20,21] and the derivation of the NDTRM by ST15. In Section 3 we derive a form of the TRM, which we refer to as the "Neutral TRM" (NTRM), that approximates volume fluxes within local Neutral Surfaces. However, Neutral Surfaces are globally ill-defined [29]. This means that isoneutral volume fluxes cannot be associated with a globally-defined, stably-stratified density variable [30], which is desireable for quantifying global circulation and water mass transformation [5,25,31]. This motivates us to return to the NDTRM, which is associated with a globally-defined density variable, and consider the special case of a Neutral Density constructed with reference to the Eulerian-mean salinity, potential temperature and pressure. We show that in this case, the NDTRM reduces to the NTRM to a consistent asymptotic order in the amplitude of the isopycnal height fluctuations. This implies that the NTRM advects a Neutral Density variable defined with respect to the Eulerian-mean ocean state, again up to a consistent asymptotic order in the amplitude of the isopycnal height fluctuations. In Section 4 we evaluate our results using idealized numerical simulations, and verify that the NTRM and NDTRM more closely approximate volume fluxes within Neutral Density surfaces constructed based on the model's mean state and based on an independent reference dataset, respectively. Finally, in Section 5 we summarize and provide concluding remarks.

Approximating Volume Fluxes within Neutral Density Surfaces
There is a substantial existing literature on the subjects of Neutral Surfaces and the TRM approximation (e.g., [20,21,24,26]). To achieve a self-contained presentation, this section summarizes relevant concepts from this literature, and from the derivation of the NDTRM by ST15.

The Temporal Residual Mean (TRM) Streamfunction
The TRM transport is fundamentally related to fluctuations of material or quasi-material surfaces, typically defined in terms of a density/buoyancy variable [20,21]. While it is impractical to quantify the Lagrangian-mean transport of fluid parcels in the ocean in the most general sense [32], a close approximation can be obtained by exploiting vertically stratified, quasi-material surfaces, commonly defined via a suitable density or buoyancy variable [22]. Consider a series of vertically monotonic, temporally-evolving, quasi-material surfaces identified by a label γ 0 , currently not assigned any physical meaning. We define a streamfunction Ψ at any horizontal point (x 0 , y 0 ) on surface γ 0 as the lateral mass flux above that surface, where u = (u, v) is the horizontal velocity vector and we have assumed a flat ocean surface (z = 0) for simplicity. The overline indicates a time average, but could be interpreted more broadly as an average over an ensemble of realizations of the flow [33]. The TRM approximates the mass flux defined in (1) under the assumption of small fluctuations in the quasi-material surfaces, i.e., z 0 = z 0 + z 0 with z 0 = O(ε) assumed to be asymptotically small.
Here ε may be interpreted as the magnitude of the surface fluctuation relative to some vertical length scale such as the ocean depth H, so ε = z 0 /H 1. Using Taylor expansions of (1), it may be shown (e.g., [21,23]) that Ψ may be approximated in Eulerian coordinates as a sum of mean and "eddy" components, defined as Here, and henceforth in this manuscript, the rules of Reynolds averaging are assumed to apply unless otherwise stated, i.e., • = • and • = 0 for any •. Note that Ψ and Ψ are streamfunctions in the sense that their curl is equal to a non-divergent velocity field. To write (2) in terms of only mean quantities, and therefore to create a useful approximation to the ocean's Lagrangian-mean velocities, we require a suitable approximation to the surface fluctuations z 0 .
Correlations of quantities such as u z 0 and z 0 2 are not routinely provided as output from ocean models formulated in Eulerian vertical coordinate systems, partly due to the ambiguity in the definition of z 0 , so instead z 0 is typically related to Eulerian fluctuations of thermodynamic variables (e.g., [22,23]). For a vertically monotonic label γ, we can further Taylor-expand our continuum of labels γ to obtain an O(ε 2 ) approximation to z 0 [21,22], allowing us to express the TRM "eddy" streamfunction in terms of purely Eulerian variables, Note that all variables on the right-hand side are evaluated at the mean surface elevation z 0 . Additionally, note that the O(ε 3 ) error terms on the right-hand side of (6) differ from those on the right-hand side of (4); additional O(ε 3 ) terms have been included in the error term following the substitution for γ using (5). In the interst of brevity, throughout this article we will write the error terms in terms of their lowest-order dependence on the relevant asymptotically small parameters without specifically identifying when the mathematical definition of the error terms has changed. In the form of equation (6), the TRM streamfunction is in principle straightforward to calculate from ocean model output if we associate γ with, for example, potential density or Neutral Density ( [24], discussed in Section 2.2). However, potential density surfaces are not particularly neutral, while the Neutral Density variable of [24] is too computationally expensive for practical use in ocean models [23].

The Neutral Density Temporal Residual Mean (NDTRM)
ST15 derived alternative expressions to (5) for the vertical fluctuations (z 0 ) of isopycnals of a Neutral Density variable (γ). These expressions contain no explicit dependence on γ itself, and so circumvent the need to explicitly compute a Neutral Density variable in order to construct the TRM. Here we briefly summarize the portions of this derivation and the resulting expressions for z 0 required in Section 3.
The Neutral Density variable of [24] is assigned based on geographically distributed reference casts. Each reference cast is pre-labeled with Neutral Density γ as a function of depth, such that the cast potential temperature, salinity and pressure may be written as θ c (γ), S c (γ) and p c (γ), respectively. To assign a Neutral Density label to a given point in space and time, the potential temperature (θ 0 ), salinity (S 0 ) and pressure (p 0 ) at that point are compared with a nearby reference cast. Specifically, the assigned value of γ is that which satisfies the "discrete neutral relation", which ST15 wrote as This relation states that water parcels with properties (S 0 , θ 0 , p 0 ) and (S c (γ), θ c (γ), p c (γ)) would have equal densities if both were moved adiabatically and isentropically to the mid-point pressure between the two parcels. Here β m and α m are short-hands for i.e., the saline contraction and thermal expansion coefficients evaluated at the thermodynamic mid-point (S m , θ m , p m ), which is defined as The small parameter ∆ measures deviations of the properties (S 0 , θ 0 , p 0 ) from the mid-values (S m , θ m , p m ). See [23] for further details on the definition of ∆ and the asymptotic derivation of (7). To derive an expression for the vertical fluctuations of a Neutral Density surface, ST15 considered the semi-Lagrangian evolution of the thermodynamic properties on such a surface, i.e., S 0 = S(z 0 (t), t), θ 0 = θ(z 0 (t), t), and p 0 = p(z 0 (t), t). They then posed Taylor expansions of (7) in terms of ε and ∆, e.g., where S (z 0 , t) = S(z 0 , t) − S(z 0 , t) is the Eulerian fluctuation of the salinity at z = z 0 . In this article we restrict our attention to the Boussinesq case for simplicity. This effectively corresponds to replacing p 0 , p c and p m by z 0 , z c and z m , where z c is the depth of the Neutral Density label γ on the reference cast and z m = 1 2 (z 0 + z c ). In this case, the highest-order NDTRM approximation to z 0 (the "NDTRM2" of ST15) is Here we use shorthand notation for the Eulerian-mean salinity and potential temperature at constant height, and for the haline contraction and thermal expansion coefficients evaluated as functions of Eulerian-mean properties, Equation (13) approximately relates the Neutral Surface height fluctuations to Eulerian fluctuations of S and θ, in analogy with (5). By excluding terms from (13), ST15 derived lower-order approximations to z 0 , referred to as the "NDTRM1", and "NDTRM0" where we define in analogy with (15). For brevity we omit complete expressions for the NDTRM streamfunctions, derived by substituting (13)- (17) into (2), and a discussion of the evaluation of the "cast" terms S c and θ c ; we refer the reader to ST15 for information on such specifics of the NDTRM. These expressions differ from those given by ST15 via the introduction of an additional small parameter,∆, that measures deviations of the Eulerian mean properties, S and θ, from the cast properties, S c and θ c . Thus, conceptually, the three small parameters in (13)-(17) are related via ∆ ∼∆ + ε.

Approximating Volume Fluxes within Local Neutral Surfaces
The Neutral Density variable of [24] (γ JM97 ) was formulated with the aim of obtaining a density-like variable whose surfaces lie as parallel as possible to the local neutral tangent plane everywhere in the ocean. However, due to the ocean's equation of state, is not possible to construct such a variable exactly [29,30]. Thus in general the vertical fluctuations of Neutral Density surfaces differ from from fluctuations of local Neutral Surfaces [34]. In Section 1 we posited that local Neutral Surfaces are best suited to estimating semi-Lagrangian mean transport, and therefore also to formulating the TRM, because adiabatic flows are constrained to follow these surfaces. The validity of these assumptions is discussed further in Section 5.
In this section we aim to formulate a TRM that approximates fluxes along local Neutral Surfaces as closely as possible (the NTRM). We first derive an approximation to the vertical fluctuations of local Neutral Surfaces, denoted z NS , for comparison with the NDTRM approximations to z 0 given in Section 2. This approximation alone yields an incomplete TRM because the resulting streamfunction (2) is not associated with a globally-defined set of surfaces (i.e., γ). To circumvent this issue we consider a Neutral Density variable γ mean whose reference dataset (i.e., S c , θ c and p c ) is chosen to be the local Eulerian-mean thermodynamic properties (i.e., S, θ and z). For such a Neutral Density variable the NDTRM reduces to the NTRM, up to the same asymptotic order in ε, which implies that the NTRM does in fact advect a globally-defined density variable, γ mean . We also show that the NTRM coincides with a TRM derived based on fluctuations of surfaces of locally-referenced potential density, up to the same asymptotic order in ε.

Vertical Fluctuations of Local Neutral Surfaces
Consider a local Neutral Surface with instantaneous height z = z 0 (t) at a fixed horizontal location (x 0 , y 0 ). To simplify our presentation we will drop all dependence on horizontal location in what follows. We define the time-evolution of z 0 (t) via the semi-Lagrangian neutral relationship [26], where δS and δθ denote infinitesimal changes in the salinity and potential temperature on the surface. Both β and α are themselves evaluated at the local salinity S(z 0 (t), t), potential temperature θ(z 0 (t), t) and pressure p(z 0 (t), t). Dividing by an infinitesimal unit of time δt, Equation (19) may be expanded to describe the time evolution of z 0 , We now use (20) to derive a relationship between the Neutral Surface height fluctuations z 0 = z 0 (t) − z 0 and the Eulerian fluctuations of S and θ. Note that fluctuations of S and θ may not only be associated with vertical heaving of Neutral Surfaces; they may also result from lateral stirring of property gradients along Neutral Surfaces. To simplify the presentation we assume that the magnitudes of both vertical heaving-induced and lateral stirring-induced property fluctuations are small and characterized by the same small parameter ε. A caveat to this assumption is that even outside of the ocean's surface and bottom boundary layers there may be a non-negligible influence of surface buoyancy fluxes [35]. For further discussion of these points and a rationalization of the small parameter ε, the reader is referred to ST15.
By posing Taylor expansions of β and α as e.g., we can approximate (20) as Finally, integrating (22) with respect to t yields an expression for the local Neutral Surface height fluctuation, which we denote as z NS , in terms of S and θ , Note that the constant of integration vanishes between (22) and (23) because by definition, z 0 = S = θ = 0. Substituting (23) into (2)-(4) yields the following expression for the NTRM streamfunction, Note that in applying (24) we are free to select z 0 , for example to coincide with a model vertical grid level. Equation (23) then approximates the vertical fluctuations of a local Neutral Surface z 0 (t) whose mean elevation is z 0 , and (24) approximates the transport above z = z 0 .

Connection to the Neutral Density Temporal Residual Mean (NDTRM)
While the NTRM streamfunction (24) approximates the semi-Lagrangian mean transport, in isolation it does not quantify isopycnal/diapycnal fluxes because we have not yet associated it with a globally-defined set of surfaces (i.e., γ). It therefore does not advect a density variable, and the fluxes cannot be partitioned into adiabatic and diabatic components (e.g., [22]). We therefore now consider whether the NDTRM, which is associated with a globally-defined density variable (e.g., γ JM97 ), can be adapted to more closely approximate the fluxes within local Neutral Surfaces.
We first note that there is a close similarity between our equations for the vertical fluctuations of local Neutral Surfaces, (23), and the vertical fluctuations of Neutral Density surfaces, (13)- (17). Indeed, (23) is identical to the lowest-order form of the NDTRM, the NDTRM0 (17), except for the differing asymptotic order of the error terms, In other words, there are O(∆ 3 ), O(∆ 2 ) and O(ε∆) terms that appear explicitly on the right-hand sides of (13) and (16), and these terms allow z NDTRM2 and z NDTRM1 to more closely approximate Neutral Density surface fluctuations than z NDTRM0 . However, by treating these O(∆ 3 ), O(∆ 2 ) and O(ε∆) terms as error terms in (17) we actually obtain a closer approximation of local Neutral Surface fluctuations (i.e., of (23)). Motivated by this observation, we now seek a hypothetical Neutral Density variable that eliminates the O(∆ 2 ), O(∆ 2 ) and O(ε∆) terms in the NDTRM. This is achieved by constructing a Neutral Density variable, referred to as γ mean , using the Eulerian-mean thermodynamic properties in place of the reference dataset, Here we have implicitly assumed that the local water column has been pre-labeled with γ based on the Eulerian-mean thermodynamic properties, and the reference elevation z r (γ) is the elevation in this water column that has been labeled with Neutral Density γ. The isopycnal fluctuation z = z 0 (t) is therefore defined by the discrete neutral relation (7) in the form Here β m and α m are defined as in (8), but with the thermodynamic mid-point defined as z m (t, γ) ≡ (z 0 (t) + z r (γ))/2.
We have now established that the transformation described by (26) and (28)-(30) leaves the NDTRM formulae (13), (16) and (17) largely unchanged, but modifies the interpretation of the O(∆ 2 ), O(∆ 2 ) and O(ε∆) error terms. To approach the NTRM, we now demonstrate that This implies that the difference (∆) between the "parcel" properties on the heaving Neutral Density surface and the properties of the reference "cast" (consisting of the Eulerian-mean thermodynamic properties) necessarily scales with the amplitude (ε) of the isopycnal height fluctuations in the limit ε → 0. To derive (32), we first note that the reference elevation must approach z 0 in the limit of vanishingly small ε, i.e., z r (γ) → z 0 as ε → 0, because in this limit z 0 (t) → z 0 , and the (stationary) isopycnal must satisfy the neutral relation with itself. We therefore write z r = z 0 + δz r , where δz r = O(∆) by definition, and seek a relationship between δz r and ε by posing a Taylor expansion of (27) in terms of ε and∆. First, we note that the midpoint thermodynamic properties can be written as where we have moved all terms in S , θ and z 0 into the error terms on the right-hand sides of these equations. We can write β m and α m as where β and α are defined as in (18). Substituting (33) and (34) into (27), Taylor expanding S(z r (γ)) and θ(z r (γ)) in terms of δz r , and Taylor expanding S(z 0 (t), t) and θ(z 0 (t), t) as in (12), we obtain Noting again that the left-hand side of (35) is O(∆) by definition, we conclude that δz r = O(ε) and that (32) holds. Finally, we combine the transformation defined by (26) and (28)-(30) with the asymptotic relationships (33), (34) and (32) to obtain the following transformation of the variables in (13), (16) and (17): Substituting (26)-(36c) into (13)- (17), we find that all forms of the NDTRM collapse to the NTRM, to a consistent asymptotic order in ε, This implies that, to a consistent order in ε, the NTRM advects a Neutral Density variable γ mean constructed via the discrete neutral relation (7), using the Eulerian-mean thermodynamic properties (S, θ, z 0 ) as the reference dataset.

Equivalence to Fluctuations of Locally-Referenced Potential Density Surfaces
We now show that our Equation (23) for the vertical fluctuations of Neutral Surfaces also coincides with the fluctuations of surfaces of locally-referenced potential density, to the same order of accuracy in ε. Ref. [36] used a TRM based on locally-referenced potential density fluctuations, but did not demonstrate its equivalence to fluctuations of local Neutral Surfaces.
We define locally-referenced potential density σ z 0 as the density referenced to the mean height z 0 , Here ρ is the in situ density and z 0 (t) is the instantaneous height of the σ z 0 surface whose mean depth is z 0 . Below we identify this σ z 0 surface to a consistent order of asymptotic approximation in the amplitude of the surface height fluctuation, ε.
Posing Taylor expansions of S and θ following (12), we may expand (40) as Taking the mean of (41) yields an expression for the semi-Lagrangian-mean locally-referenced potential density As z 0 is defined by a surface of constant locally-referenced potential density, it follows that Substituting (43) and (42) into (41) eliminates the first two terms in (41) to a consistent order of approximation in ε, and thus (41) reduces to where z LRPD denotes vertical fluctuations of the locally-referenced potential density surface. Therefore, to the same asymptotic order of approximation in ε, the NTRM estimate of isopycnal height fluctuations also coincides with vertical fluctuations of locally-referenced potential density surfaces.

Comparison and Assessment Using an Idealized Numerical Model
The key result of Section 3 is that the NTRM streamfunction (24) approximates volume fluxes within local Neutral Surfaces, and advects a Neutral Density variable defined using the Eulerian-mean thermodynamic properties (S, θ, z 0 ) as a reference dataset. An implication of this result is that while the NDTRM should more accurately approximate volume fluxes within surfaces of γ JM97 , the NTRM should more accurately approximate volume fluxes within surfaces of γ mean . We now test this prediction using output from an idealized numerical model of the Antarctic continental shelf and slope. We use this model because its geometry simplifies the task of computing γ mean , and because nonlinearities in the seawater equation of state lead to relatively pronounced deviations of potential density and Neutral Density surfaces in the simulated region [23,37,38].

Model Configuration
The modeling approach has been described in detail by [23,[39][40][41][42], so here we present only salient aspects of its configuration and refer the reader to previously published articles. The MIT general circulation model [43,44] was configured in a zonally re-entrant channel of length L x = 400 km and widh L y = 450 km, with its depth ranging from 500 m over the continental shelf at the southern boundary to 3000 m at the northern boundary. The model was forced via a combination of a steady westward surface wind stress, linear drag at the sea floor, a two-equation representation of surface sea ice heat and salt exchanges (see [42]), a constant input of salt over the southernmost 50 km of the model domain, and restoring toward profiles of θ and S derived from hydrography within a sponge layer covering the northernmost 50 km of the domain. The model's horizontal grid spacing is around 1 km, sufficient to resolve the eddy field [41], and the vertical grid consists of 53 vertical levels with spacings increasing from around 13 m at the surface to around 100 m at the sea floor. The analysis conducted here uses 5 years of daily model snapshots from a period in which the simulation had reached statistically steady state.

Constructing Neutral Density from the Model's Mean State
To test the NTRM we must first construct γ mean , i.e., for each time t n and model grid location (x i , y j , z k ) we must assign Neutral Densities γ mean (x i , y j , z k , t n ), where i = 1 . . . N x , j = 1 . . . N y , k = 1 . . . N z , and n = 1 . . . N t are gridpoint indices. Similar to the procedure laid out by [24], we require a dataset to which the instantaneous model stratification can be referred in order to assign a γ mean value to each grid point. In [24], this dataset was the hydrographic climatology of [45] and a corresponding set of neutral density labels on a coarse global grid. Similarly, for the purpose of constructing γ mean the reference dataset is the model's time-mean hydrography, S(y j , z k ) and θ(y j , z k ) (see Figure 1a,b), and a corresponding set of neutral density labels γ mean ref (y j , z k ). Note that the reference dataset is approximately independent of x due to the zonally symmetric model geometry and forcing. The approximate two-dimensionality of the time-mean hydrography also implies that the neutral helicity of the reference dataset is zero [29]. Taken together with the simply connected model domain, this implies that Neutral Surfaces are well defined [46].
We now outline the procedure for labeling the reference dataset with γ mean ref and assigning γ mean values to instantaneous model output. This procedure closely follows that described by [24], to which the reader is referred for more specific details of the algorithm. We first labeled the northernmost model water column, γ mean ref (y N y , z k ), by setting it equal to γ JM97 (y N y , z k ). Here γ JM97 was calculated from the hydrographic measurements used to restore the model stratification in the offshore sponge layer, taking the model domain to lie along a slice between (61W,67S) at y = 0 and (50.6W,67S) at y = 450 km [23]. We then iterated southward through the model water columns, assigning the j th column, γ mean ref (y j , z k ) | k = 1 . . . N z , from the (j + 1) th column using the form of the discrete neutral relation (7) given by [24] and linear vertical interpolation. For grid points in the j th column that were denser or lighter than any grid point in the (j + 1) th water column, we linearly extrapolated γ mean ref (y j , z k ) vertically. Note that in our vertical density extrapolation we assigned the vertical gradient of γ mean ref to be equal to the vertical gradient of locally-referenced potential density, whereas in this region the vertical gradient of γ JM97 was typically assigned to be twice the vertical gradient of locally-referenced potential density (the "b-factor" of [24]). This choice of extrapolation was made for simplicity and should not influence our results, though a consequence is that γ mean spans a smaller range of density values than γ JM97 in this model domain. Figure 1. Diagnostics of the mean model state, where averages are taken in time (over 5 years of daily snapshots) and in the along-slope (zonal) direction (400 km distance at ∼1 km horizontal grid spacing. (a) Potential temperature θ, (b) practical salinity S, (c) Neutral Density constructed based on the model's mean state, which we denote as γ mean (see Section 4), (d) difference between γ JM97 , where γ JM97 is the Neutral Density calculated using the reference dataset and algorithm of [24], and γ mean .
To calculate γ mean (x i , y j , z k , t n ) for a given model gridpoint with salinity S(x i , y j , z k , t n ) and potential temperature θ(x i , y j , z k , t n ), we construct continuous functions of z, S(y j , z) and θ(y j , z), via linear interpolation. We then find a depth z r such that (S(y j , z r ), θ(y j , z r ), z r ) and (S(x i , y j , z k , t n ), θ(x i , y j , z k , t n ), z k ) satisfy the discrete neutral relation. We then assigned γ mean (x i , y j , z k , t n ) by linear vertical interpolation of γ mean ref (y j , z k ) to the point z ref . Note that our algorithm for assigning γ mean differs slightly from the algorithm of [24] for assigning γ JM97 : to assign γ mean we refer each gridpoint to a horizontally local "cast" consisting of the Eulerian-mean model properties, whereas the algorithm of [24] refers each gridpoint to a four geographically local reference casts. Additionally, we use linear vertical interpolation to assign a γ JM97 value, whereas the algorithm of [24] uses quadratic interpolation.
In Figure 1a,b we plot the time-and zonal-mean potential temperature and salinity, θ and S, that serve as the reference dataset in the construction of γ mean . The mean state is qualitatively similar to hydrographic measurements from the western Weddell Sea [41], and exhibits a layer of relatively cold, fresh, shelf-sourced bottom water that descends down the continental slope. Figure 1c shows the time and zonal mean of the Neutral Density constructed from the model's mean state, i.e., γ mean , while Figure 1d shows the difference between γ JM97 and γ mean . As expected, γ JM97 takes lower values than γ mean in waters lighter than the lightest water at the northern boundary, and higher values than γ mean in waters denser than the densest water at the northern boundary, due to the procedure for extrapolating γ mean ref . Some variations are visible across the continental slope due to the presence of alternating along-slope jets that slowly drift offshore, and which do not perfectly cancel under integration over the 5-year analysis period [42].
We anticipate that γ mean should be more closely aligned with local Neutral Surfaces than γ JM97 . To illustrate this, in Figure 2 we plot the γ mean , γ JM97 , θ and S at the sea floor as a function of cross-slope distance, i.e., within the dense outflowing layer of shelf-sourced bottom water. Figure 2 shows that γ JM97 exhibits spurious fluctuations and local minima, on the order of 0.01 kg m −3 , that are absent in γ mean . For example, around y = 200 km there are negligible variations in θ and S, indicating that the sea floor should coincide with a neutral surface. This is reflected in γ mean , but not in γ JM97 , which exhibits a spurious local minimum at this point. We conclude that iso-surfaces of γ mean are more closely aligned with local Neutral Surfaces than γ JM97 , though the spurious variations in γ JM97 are nonetheless quite small compared to the domain-scale ranges of γ mean and γ JM97 .  [42], as a function of distance from the shore. This plot illustrates that the use of an independent reference dataset leads to a spurious local minima, e.g., around y = 200 km, in γ JM97 . These are absent in γ mean , indicating that it more closely tracks Neutral Surfaces. The rapid density variations around y = 100 km are due to "steps" in the model's geopotential-coordinate representation of the continental slope.

Comparing Diagnosed and TRM Transports
Having constructed γ mean , we now compare the transports within γ mean and γ JM97 surfaces, as estimated by the TRM and diagnosed directly from the model simulations. Based on Section 3, we anticipate that the NTRM should more closely approximate the diagnosed fluxes within γ mean surfaces, while the NDTRM should more closely approximate the diagnosed fluxes within γ JM97 surfaces.
We diagnose fluxes within density surfaces via direct calculation of (1). The procedure is summarized here to achieve a self-contained presentation; for further details the reader is referred to ST15. For each model snapshot, the model-reported meridional volume fluxes and Neutral Density (γ mean or γ JM97 ) on geopotential surfaces are interpolated to a vertical grid that is 10 times finer than the model's vertical grid (see [23,41,42]). The volume fluxes are then assigned to a discrete set of density bins, averaged in time to compute mean volume fluxes within density surfaces, and integrated vertically in density space to obtain the streamfunction. To map the streamfunction back to z-space, we also average the height z 0 of each density surface, and map the streamfunction at density γ 0 to height z 0 . This allows for direct comparison with the TRM, which approximates the transport above the isopycnal whose mean elevation is z 0 . Note that due to the zonal symmetry of the model domain, we only calculate the meridional component of the streamfunction, which quantifies circulation in the y/z plane. We denote the streamfunction corresponding to volume fluxes in γ mean surfaces as ψ mean , and the streamfunction corresponding to volume fluxes in γ JM97 surfaces as ψ JM97 .
We compare ψ mean and ψ JM97 with the NTRM and NDTRM estimates of the volume fluxes. The NTRM is evaluated as the meridional component of (24), with the mean operator defined as a time and zonal average. We denote the resulting streamfunction as ψ NTRM . The NDTRM streamfunction is defined following equation (B.4) of ST15, and is denoted as ψ NDTRM . Throughout this section we use "NDTRM" to refer to the NDTRM2 of ST15, which was shown to exhibit closer agreement with diagnosed fluxes in γ JM97 surfaces than lower-order formulations of the NDTRM.
In Figure 3 we compare ψ mean and ψ JM97 , ψ NTRM and ψ NDTRM directly. Figure 3a shows the overturning circulation diagnosed from volume fluxes within γ mean surfaces. The circulation is characterized by wind-driven shoreward transport within the surface mixed layer and shoreward interior eddy transport across the continental slope [42]. These shoreward transports are returned offshore via dense shelf water export down the continental slope. Figure 3 quantifies the difference in the computed transports within γ mean and γ JM97 surfaces. The most substantial difference is that the amplitude of ψ JM97 is enhanced relative to ψ mean over the upper continental slope and the continental shelf, exceeding the ∼0.23 Sv shelf overturning magnitude of ψ mean by around 0.08 Sv. This suggests that a non-negligible fraction of the diagnosed shelf overturning in ψ JM97 surfaces is due to errors introduced by the independent reference dataset of [24]. Figure 3c,d show the differences between ψ mean and ψ NTRM , and between ψ JM97 and ψ NDTRM , respectively. Both TRM streamfunctions somewhat underestimate the amplitude of the overturning circulation over the continental shelf, and also exhibit errors close to the ocean surface and floor, where the vertical stratification is weak. The latter might be ameliorated by adopting a generalized form of the NTRM that utilizes vertical tracer fluxes in weakly-stratified regions, analogous to the methods proposed by [47,48]. The errors are typically smaller than 0.05 Sv in magnitude, which is smaller than the difference between ψ mean and ψ JM97 over the shelf, suggesting that ψ NTRM more closely approximates ψ mean and that ψ NDTRM more closely approximates ψ JM97 . (a) Diagnosed transport within γ mean surfaces (ψ mean ), (b) difference between diagnosed transports within γ JM97 surfaces and γ mean surfaces (ψ mean − ψ JM97 ), (c) difference between NTRM-estimated transport within γ mean surfaces and diagnosed transport within γ mean surfaces (ψ NTRM − ψ mean ), (d) difference between NDTRM-estimated transport within γ n surfaces and diagnosed transports within γ JM97 surfaces (ψ JM97 − ψ NDTRM ). The diagnosed streamfunctions, ψ mean and ψ JM97 , have been mapped from density to z coordinates using the time-mean depths of γ mean and γ JM97 surfaces, respectively.
To quantify the agreement between the TRM and diagnosed transports, in Figure 4a we compare the net southward transport reported by each streamfunction, as a function of latitude. This calculation follows ST15, and only includes southward transport at any given latitude if that transport has been supplied from the northern boundary of the domain, i.e., only counting streamlines that connect that latitude to the northern boundary. Figure 4b shows a similar comparison, but only accounting for southward interior transport, i.e., excluding transport in the surface mixed layer. In each case, ψ NTRM and ψ NDTRM more closely tracks ψ mean and ψ JM97 , respectively, over the continental slope (y 150 km). Over the continental slope (150 km y 250 km) the transports are difficult to distinguish, as suggested by Figure 3. In the deep ocean (y 250 km) the TRM and diagnosed transports diverge due to large differences in the streamfunctions close to the sea floor (see Figure 3c,d). Thus, our net southward transport metrics confirm the result suggested by the overturning streamfunction comparison in Figure 3: the NTRM more closely approximates volume fluxes in ψ mean surfaces, while the NDTRM more closely approximates fluxes in ψ JM97 surfaces.

Advective Tracer Transport
Finally, we note that that in addition to estimating isopycnal/diapycnal transports, TRM streamfunctions may also be used to estimate the advective flux of tracers. Importantly, no globally-defined density variable is needed to perform such an estimate, so, e.g., it would not be necessary to construct a Neutral Density variable like γ mean in conjunction with the NTRM for this purpose. We now compare the advective meridional fluxes of θ and S estimated by the NTRM and NDTRM with those diagnosed from the model simulation. We focus on the continental shelf and open ocean regions (y 150 km), in which the heat and salt transports have been shown to be primarily advective, rather than diffusive [42]. Figure 5a,b shows the northward eddy fluxes of heat and salt, diagnosed directly from the daily model output via Here χ = θ or χ = S, z = −h(y) is the sea floor elevation, the mean operator • denotes a time and zonal average, and primes denote deviations from that average. The heat and salt fluxes have been rescaled to take units of TW and Gg/ s, respectively, for ease of interpretation. We calculate NTRM and NDTRM estimates of the advective fluxes via Here we set ψ TRM = ψ NTRM and ψ TRM = ψ NDTRM to define the flux estimates F χ NTRM and F χ NDTRM , respectively. In Figure 5c,d we plot the differences between the diagnosed and TRM-estimated heat and salt fluxes, which typically differ by around 10-20%. There is very little difference between the NTRM and NDTRM salt fluxes, whereas the NTRM yields errors in the heat flux that are up to 40% smaller than the errors incurred by using the NDTRM. Fluxes are plotted over a subset of the model domain in which the eddy transport was found to be primarily advective, rather than diffusive [42]. (a,b) Diagnosed offshore heat and salt fluxes, respectively. (c,d) Differences between diagnosed and NTRM-/NDTRM-estimated offshore heat and salt fluxes, respectively.

Summary and Conclusions
The Temporal Residual Mean (TRM) transport requires an Eulerian approximation for vertical fluctuations of a set of quasi-material surfaces (e.g., [20][21][22]). The TRM approximates the transport within these surfaces, which in turn approximates the semi-Lagrangian mean and Lagrangian-mean transports [9]. This article builds on a recent study by [23] (hereafter ST15), who used the Neutral Density variable γ JM97 of [24] to construct a Neutral Density TRM (NDTRM) that can be applied anywhere in the ocean, circumventing the limitations of e.g., potential density at high latitudes. A key result of this article is that the vertical fluctuations of local Neutral Surfaces, which are defined directly via the continuous form of the neutral relation, can be approximated via (23). The corresponding TRM streamfunction (24), referred to as the "Neutral TRM" (NTRM) should more closely approximate the semi-Lagrangian mean transport than the NDTRM, and therefore also better approximate the generalized Lagrangian mean transport [9]. The NTRM also coincides with a TRM based on vertical fluctuations of locally-referenced potential density surfaces, again to the same asymptotic order in the amplitude of the isopycnal height fluctuations.
In isolation the NTRM streamfunction does not directly quantify isopycnal and diapycnal ocean transports (see e.g., [31]) because local Neutral Surfaces are, by definition, not globally well-defined. However, in Section 3 we showed for the special case of a Neutral Density variable (γ mean ) constructed using the Eulerian-mean ocean state, the NDTRM reduces to the NTRM, up to the same asymptotic order in the amplitude of the isopycnal height fluctuations. Thus the NTRM approximately quantifies adiabatic and diabatic volume fluxes within and across γ mean surfaces, again up to a consistent asymptotic order in the amplitude of the isopycnal height fluctuations.
In Section 4 we tested this theoretical prediction explicitly using an idealized eddy-resolving simulation of the Antarctic continental shelf and slope. Using the model's time-mean state as a reference dataset for γ mean , we labeled the model output with both γ mean and γ JM97 and computed volume fluxes within surfaces of both Neutral Density variables. We showed that γ mean iso-surfaces were more closely aligned with local Neutral Surfaces than γ JM97 , indicating that fluxes within γ mean surfaces should be more closely aligned with the local neutral tangent plane. Consistent with Section 3, we found that the NTRM more closely approximated volume fluxes within γ mean surfaces, while the NDTRM more closely approximated volume fluxes within γ JM97 surfaces. Finally, we compared the advective fluxes of heat and salt estimated by the NTRM and NDTRM, with those diagnosed from the model. Both formulations of the TRM yielded flux estimates in close agreement with the diagnosed fluxes, though the NTRM yielded a somewhat closer approximation to the heat fluxes.
The results summarized above suggest that for most oceanographic purposes, the NTRM is preferable to the NDTRM for the purposes of estimating isopycnal and diapycnal volume fluxes, and of estimating advective tracer fluxes. It also bears the advantage of being more convenient to compute than the NDTRM, as it makes no explicit reference to globally-defined Neutral Density variable or to the set of reference casts used to defined such a variable. A practical consideration is that the NTRM advects a Neutral Density variable constructed from the model's time-mean state, i.e., γ mean [22]. Though algorithms for constructing global, approximately neutral density variables have been proposed (e.g., [24,49,50]), translating these algorithms to the time-mean state of a three-dimensional model with an arbitrary ocean geometry is non-trivial. More recent algorithms have been proposed to construct global, approximately neutral surfaces [46,51], but have not yet been adapted to the construction of three-dimensional Neutral Density variables. Furthermore, a Neutral Density variable such as γ mean would need to be constructed with great care to achieve consistency with γ JM97 , which has been used to map the global overturning circulation [5,25].
Throughout this manuscript we have assumed that local Neutral Surfaces coincide with purely adiabatic flows [26], and are therefore optimal for the construction of the TRM. Neutral helicity precludes the construction of a continuous density variable that lies exactly parallel to these surfaces everywhere in the ocean [29]. However, on the basis of two-parcel energetic considerations, Ref. [52] argued that that adiabatic parcel displacements should not, in fact, follow the neutral tangent plane. Variance of potential temperature and salinity has been found to be larger on Neutral Density surfaces than on surfaces of potential density referenced to 2000 decibars [26,53], though this does not explicitly contradict the notion that adiabatic parcel displacements follow neutral surfaces. Although this topic remains under debate [54][55][56], an implication is that there may yet be a more accurate choice than local Neutral Surfaces to define adiabatic heaving, and therefore to approximate the semi-Lagrangian mean and generalized Lagrangian mean transports. Further work is required to compare estimates of thickness-weighted average transports within different formulations of quasi-neutral surfaces, such as ω-surfaces [51], thermodynamic neutral density (γ T ) surfaces [50], and orthobaric density surfaces [49].
Previous studies have used the NTRM0, and thus implicitly the NTRM (see Section 3), to estimate volume fluxes within γ JM97 surfaces (e.g., [42,[57][58][59]). In principle, this incurs an O(∆) error in the approximation, where∆ quantifies differences between the model's mean state and the independent reference dataset. On the other hand, our comparison in Section 4 shows that the NTRM approximates transports within γ JM97 surfaces almost as closely as it approximates transports within γ mean surfaces, so in practice such errors may not be of concern. However, further work is required to evaluate errors associated with neutral helicity when quantifying transports within γ mean surfaces using the NTRM.