A method is developed for estimating the along-isopycnal and vertical mixing coefficients (K and D) and the absolute velocity from time-averaged hydrographic data. The method focuses directly on transports down tracer gradients on isopycnals. When the tracer considered is salinity or an appropriate variable for heat, this downgradient transport constitutes the along-isopycnal component of the thermohaline overturning circulation. In the method, a geostrophic streamfunction is defined that is related on isopycnals by tracer contours and by the thermal wind relationship in the vertical. Volume and tracer conservation constraints are also included. The method is overdetermined and avoids much of the signal-to-noise error associated with differentiating hydrographic data in conventional inverse methods. The method is validated against output of a layered model. It is shown to resolve both K and D, the downgradient isopycnal transport, and the mean flow on isopycnals in the North Pacific and South Atlantic.
Importantly, an understanding is established of both the physics underlying the method and the circumstances necessary for an inverse method to determine the mixing rates and the absolute velocity. If mixing is neglected, the method is the Bernoulli inverse method. At the limit of zero weight on the tracer-contour equations the method is a conventional box inverse method. Comparisons are drawn between each method and their relative merits are discussed. A new closed expression for the absolute velocity is also presented.
How abyssal water masses are returned to shallower regions of the ocean is not fully understood. The two mechanisms identified as possible sources of this “pull” of water back toward the ocean surface are (i) diapycnal upwelling due to vertical mixing, either throughout the deep ocean or at particular hot spots around the globe, and (ii) wind-driven upwelling along sloping isopycnals in the Southern Ocean. Both are elements of ocean dynamics for which our observations are sparse and inconclusive [for a comprehensive review see Kuhlbrodt et al. (2007)]. A detailed global understanding of the diffusive and advective processes that contribute to the modification of water masses and the overturning circulation is needed.
The most established method for diagnosing the general circulation, from hydrographic data, is the box inverse method (Wunsch 1978). Although it is able to accurately resolve major components of the circulation, box inversions of hydrographic data have seldom incorporated vertical mixing and along-isopycnal mixing coefficients as unknowns. The few that do have difficulty resolving such values for areas of the Southern Ocean where isopycnals outcrop either because of the added uncertainty of air–sea fluxes or because interior and mixed layer processes cannot be separated (Ganachaud et al. 2000; Sloyan and Rintoul 2000). The Bernoulli or streamfunction inverse method of Killworth (1986) and the beta-spiral method of Stommel and Schott (1977) rely on the isopycnal gradient of potential vorticity (PV), or other tracers, changing direction in the vertical and, in their classical form, assume that the flow is nondiffusive. Although these methods are applied in different ways, it can be shown that they are based on the same physical principles, particularly the beta-spiral and box inverse methods (Davis 1978). The tracer-contour inverse method presented here draws on aspects of each of the box, Bernoulli, and beta-spiral inverse methods but is designed so that mixing processes are of leading order importance.
In this article, we establish an inverse method that allows us to accurately quantify both the magnitude and distribution of turbulent vertical mixing D, along-isopycnal mixing K, and the geostrophic flow, with particular emphasis on the velocity component down the tracer gradient (usually Θ and S on isopycnals).
This article is structured as follows: section 2 establishes the concept of a tracer contour and presents the central “tracer-contour equations,” which follow from simple combinations of conservation statements. Results from the application of the tracer-contour equations to output of the Hallberg isopycnal model (HIM) (Hallberg 2000) in a region of both the North Pacific and the South Atlantic are given in section 3. We discuss how the tracer-contour equations may be used independently in section 3b and the results are compared with a box inverse method in section 3c. The two approaches are combined to form the tracer-contour inverse method in section 3c. Comparison with the Bernoulli method is made in section 3d. The effect of random error for different relative weights on the box and tracer-contour equations is discussed in section 4. Criteria needed to diagnose the mixing coefficients and theoretical implications are discussed in section 5, including a new closed expression for the absolute velocity. Concluding remarks are made in section 6.
The tracer-contour inverse method is formulated using neutral density γ n as a vertical coordinate and conservative temperature Θ as a variable for “heat.” Here neutral density is considered as it is along neutral tangent planes that parcels of water may be exchanged adiabatically, without restoring buoyancy forces (McDougall 1987). Conservative temperature Θ is used as it is the most conservative temperature variable available. Conservative temperature is proportional to potential enthalpy, which represents the “heat content” per unit mass of seawater (McDougall 2003).
Where this study discusses output of the HIM, potential density σ2 (i.e., potential density referenced to a depth of 2000 m) and potential temperature θ are used. In HIM, unlike in the real ocean, properties are exchanged adiabatically along σ2 surfaces and potential temperature is a perfectly conserved quantity. Note, however, that the distinction between conservative temperature and potential temperature and neutral and potential density is not central to this study. We will frequently refer to Θ or θ as temperature and neutral or potential density surfaces as isopycnals.
2. Tracer contours
Consider the idealized scenario in which no mixing or diffusive processes occur in the ocean. In such an ocean the path taken by a parcel of water of a given tracer concentration C and neutral density is simply the path where the C and γ n are constant (ignoring subtle effects on γ n owing to the changes in pressure and longitude and latitude along the path). Thus, currents in such an ocean would closely follow contours of constant tracer on isopycnals. It is clear then that the amount by which trajectories do not follow such a path, that is, the amount of flow across C contours on isopycnals and vertically through isopycnals, is determined purely by the magnitude of vertical and along-isopycnal mixing processes. Diagnosing the magnitude of such processes is thus best undertaken in a γ n versus C reference frame.
The tracer-contour inverse method, as presented here, considers contours of constant temperature and salinity on isopycnals, as these are the best observed tracers in the ocean. Both temperature and salinity form the same contours on isopycnals. They are effectively passive tracers on isopycnals since mixing them along neutral tangent planes does not affect the stratification. The isopycnal flow across these contours (the cross-contour flow) is the along-isopycnal component of the thermohaline overturning circulation. In principle, the method can be generalized to consider contours of any conservative tracer on isopycnals (appendix C).
In the following h is the temporally averaged “thickness” of a density layer (h = Δρ/ρz for some small Δρ). In principle the thickness h can be infinitely small. The tracers C (e.g., Θ and S) are thickness-weighted temporally averaged quantities. The thickness-weighted velocity is written /h, where v is the lateral velocity vector. The overbar denotes the temporal average. We will later separate the thickness-weighted velocity into an isopycnal-mean component, v (not thickness weighted), and a perturbation component, /h [see Equation (3)].
The thickness-weighted velocity component, directed down the along-isopycnal temperature and salinity gradient, is (/h) · n. The cross-contour direction is defined by n = ∇γΘ/|∇γΘ|, where ∇γ is the epineutral gradient. Along-isopycnal and vertical mixing act downgradient with coefficients K and D, respectively. The equation describing the balance between this flow and mixing process in steady state and, assuming no along-isopycnal gradient of K, is
Equation (1) was first derived by McDougall (1984) in a slightly different form and described as the “water mass transformation” equation. In the current form, (1) was used by Zika et al. (2009) to relate the strength of the Southern Ocean overturning circulation to diffusivities K and D. See appendix A for a full derivation of (1) and appendix C for the analogous form of (1) in terms of the tracer C.
Equation (1) represents a balance between downgradient advection and diffusion through the along-isopycnal and vertical mixing lengths and λ γ, respectively. Along-isopycnal mixing enters largely through the along-isopycnal curvature of either temperature or salinity. More precisely, in the case of a linear equation of state and no thickness gradient (i.e., if the thermal expansion and haline contraction coefficients α and β are constant and ∇γh = 0) the along-isopycnal mixing length simplifies so that . Both linear and nonlinear terms are retained here.
Notably, (1) does not involve the diapycnal velocity component wγ, vertical differences in D, nor second derivatives of tracers in z. Instead, vertical diffusion appears through the Θ–S curvature (noting that Sz3d2Θ/dS2 = SzΘzz − ΘzSzz). The Θ–S curvature is less sensitive to noise in hydrographic data due to heave than Szz and Θzz. Assuming a linear equation of state, the Θ–S curvature is proportional to the curvature of Θ or S with respect to ρ.
In appendix C, we derive a general from of (1) for advection of the tracer C (C3). Again assuming a linear equation of state, vertical mixing influences (C3) through the curvature of C with respect to ρ.
The error term represents the effect of an unsteady or “trend” term, the effect of the gradient of the along-isopycnal mixing coefficient, and any other form of error owing to measurement or sampling (εnoise); that is,
In (2), Θt|γ represents the trend in temperature at constant density and not the variability that is encompassed in the eddy coefficient K. It is straightforward to retain the ∇γK term in (1) if along-isopycnal variations in K are to be resolved.
We split the thickness-weighted velocity into an isopycnal-mean and perturbation component
It is understood that, to a good approximation, the mean circulation on an isopycnal satisfies geostrophy and that the mean velocity can be represented by a geostrophic streamfunction Ψγ (McDougall 1988) such that
In (4), f is the Coriolis frequency and k is the vertical unit vector normal to geopotential surfaces. The perturbation component of (3) may be parameterized as a downgradient flux of thickness (or bolus flux) such that
where υ* is the bolus velocity and KPV is the epineutral thickness or potential vorticity diffusion coefficient. Following McDougall (1991), a beta term, Kβ/f, where β is the meridional gradient of f, may be included on the rhs of (5). However, as HIM has only a thickness diffusion term, the beta term is omitted for the time being. Another way of writing such a parameterization is in terms of the vertical density gradient, as in Gent et al. (1995).
The component of the mean geostrophic velocity that crosses temperature and salinity contours on isopycnals is the cross-contour velocity
Substituting (1) and (5) into (6) and defining a third mixing length for thickness λh, where 1/λh = ∇γ log(h) · n, the cross-contour geostrophic velocity is then purely related to the mixing coefficients K, D, and KPV:
[see appendix A for the complete form of (7)]. Multiplying (7) by f, yields the cross-contour velocity, in terms of both along-isopycnal and vertical mixing coefficients and the gradient of the geostrophic streamfunction,
In (8), we have assumed that KPV = K. However, if desired, KPV could be retained as a third mixing variable.
Integrating along a contour of constant temperature and salinity (or any tracer) on an isopycnal from point (x1, y1) to (x2, y2) and assuming the mixing coefficients K and D are constant along the contour, we have an equation that relates the difference in streamfunction between (x1, y1) and (x2, y2) to the along-isopycnal and vertical mixing coefficients,
Equation (9) states that the total geostrophic flow crossing a Θ contour on an isopycnal, between points at either end of the contour, is the integral of the reciprocal of the scale lengths , λh, and λγ multiplied by the mixing coefficients and the Coriolis frequency. Moving between points (x1, y1) and (x2, y2) on the isopycnal, the gradient of temperature points to the left.
The error term εΨγ may be largely attributed to the integral of f. Here εΨγ also represents the correlations along the contour between coefficients K and D and the scale lengths and λ γ and the inadequacies of the steady geostrophic assumption implicit in using Ψγ.
Both Zika and McDougall (2008) and Zika et al. (2009) exploit temperature contours on isopycnals, relating downgradient advection to mixing. In this study, we consider the entire water column simultaneously and explicitly solve for a geostrophic streamfunction on one reference level.
in which specific volume anomaly δ = 1/ρ(S, Θ, p) − 1/ρ(Sconst, Θconst, p) (conventionally, Θconst = 0°C and Sconst = 35 psu), p0 and δ0 are the pressure and specific volume anomaly on the reference level (γ0), and p̃ and p̃0 are the pressures averaged over isopycnals γ and γ0, respectively.
This formulation may be extended to specifically describe the difference in Ψγ along the tracer contour on the isopycnal γ between points (x1, y1) and (x2, y2) to the difference for the same (x, y) points on the reference level γ0 (Fig. 1) such that
We may now write a single equation for Ψγ0, K, and D for the tracer contour connecting (x1, y1) to (x2, y2) on a particular isopycnal,
Equation (12) may be written for a large number of tracer contours to solve for the reference level streamfunction and individual diffusivities (on each isopycnal or some given spatial function). There are benefits in referencing to a fixed pressure and these are detailed in appendix D. This study retains an isopycnal as the reference level so that the method can be compared to the box inverse method and allows reference level minimizations to be imposed. McDougall and Klocker (2009) discuss different streamfunctions considered on different isopycnals and their corresponding errors. Here we consider isopycnals where the variation in pressure on each isopycnal is small and, hence, difference between the methods has a negligible effect on the result.
The system of equations is represented as
The matrix 𝗔 is made up of the coefficients on the lhs of (12), and b is simply each value for each tracer contour. The vector of unknowns, x, is
In the case of L boundary points with N isopycnals, a solution is sought for Ψγ0 at each L and for K and D on each N. The error ε is minimized using a singular value decomposition solver.
It is likely that, for L streamfunction points, the effective number of tracer contours that add orthogonal information is L/2 (see appendix E for a complete discussion). Thus, the approximate number of effective equations is N × L/2, resulting in an overdetermined system when N × L/2 > 2N + L.
3. Application of the tracer-contour equations
The inverse method is applied to two regions of the ocean using output of the Hallberg isopycnal model. This allows direct comparisons between the inverse method solution and the HIM diffusivities and transports. The inverse method is applied in several ways: using only the tracer-contour equations (12), using box inverse method equations with a reference level minimization, and using a combination of both box and tracer-contour equations without any reference level minimization. The final case forms the tracer-contour inverse method. A discussion of row and column weighting and the effect of random error on the solution is given in section 4.
a. The Hallberg isopycnal model and study regions
There are two main advantages of using the HIM as a means of testing the accuracy of the method. First, by utilizing a density coordinate system the HIM avoids unwanted mixing that can arise from the Veronis effect in z-coordinate GCMs (Veronis 1975). Second, since HIM is already in an isopycnal coordinate, we avoid any error implicit in interpolating on to isopycnals from a Cartesian grid of Θ, S, and p values. The HIM output used is the time average of the final 40 years of a 100-yr coupled experiment at 1° × 1° resolution. The model has a fixed diffusivity for thickness and tracers along isopycnals (K) of 600 m2 s−1. The model employs a Richardson-number-dependent vertical diffusion scheme, which results in high vertical diffusivities, O(10−4 m2 s−1), in the mixed and bottom boundary layers and relaxes to a background profile in the ocean interior. Trends in temperature, salinity, and thickness—spurious or otherwise—are corrected such that the steady-state conservation and geostrophic equations hold exactly.
The inverse method is applied to two regions, both with relatively uncomplicated topography that is everywhere deeper than 2000 m. The first is in the northeast Pacific between 10° and 30°N, 210° and 230°E. This region is characterized by a low velocity at depth (i.e., |v| ≪ 1 cm s−1 at ∼2000 m). The second region is in the southwest Atlantic between 48° and 36°S, 310° and 330°E where the Antarctic Circumpolar Current (ACC) steers northward and interacts with the Malvinas (Falkland) Current in the Scotia Sea. It can be characterized by strong depth-integrated velocity (i.e., |v| > 1 cm s−1 at ∼2000 m). No attempt will be made to compare the model simulation to observations. For simplicity, the experiments presented here do not involve tracer contours that reach the ocean surface or topography; that is, the method is applied to regions distant from continental boundaries and above the sea floor. In principle, tracer contours can be defined for such regions.
b. Contours only
To focus on diffusive processes in the ocean interior we define tracer contours on isopycnals between 500 m and two model grid points above any topographic feature (around 2600 m), the second lowest complete isopycnal being taken as the reference level (in the case of the North Pacific region this isopycnal has the lowest rms velocity). For the North Pacific region there are 56 boundary grid points (L = 56) and 15 isopycnals (N = 15); for the South Atlantic L = 31 and N = 15.
The geostrophic component of the mean velocity crossing Θ contours on isopycnals v · n is −(∇γΨγ × k) · n /f. So, the average value of v · n, along the tracer contour, is approximately −ΔΨγ/f̃Δx, where f̃ is the average f along the contour and Δx is the length of the contour. An important aspect of the tracer-contour method is its ability to determine this downgradient advection. It is instructive to present all the terms in (12), divided by f0Δx, in Θ–S space so that the production of individual water masses can be identified since each contour has a unique temperature and salinity and, thus, represents a specific location on such a diagram (Figs. 2b,c,d and 3b,c,d). In the North Pacific region the strongest cross-contour flow (with respect to the reference level) occurs in the lightest, freshest waters where the subtropical gyre is oriented almost directly southward, bringing cool North Pacific waters into the equatorial zone. This pattern is evident from the terms on the rhs of (12), which are largest in the lightest waters. This suggests that a weak reference level is a reasonable approximation in this region. In the South Atlantic region, the terms on the rhs of (12) are large throughout the whole water column. The patterns in Figs. 3b,c do not coincide with Fig. 3d, suggesting that a weak reference-level velocity is not a good approximation there.
In the possible case that the coefficients and 1/λγ are linearly dependent on each isopycnal, K and D cannot be solved for independently and other information must be included; for example, the ratio of K to D may be set. It is instructive then to plot the mixing terms on the lhs of (12) against one another. We plot versus , where K0 and D0 are equivalent to typical K and D values. These plots show that in the two regions considered here the coefficients are not linearly dependent (Figs. 4 and 5).
In the case of the North Pacific the coefficients vary smoothly with respect to one another, whereas they are much more scattered in the South Atlantic region. For the lightest isopycnals in the South Atlantic region, the coefficient of D (i.e., 1/λ γ) is close to zero because the Θ–S curvature is small relative to the curvature along the isopycnal (Fig. 5). The D diagnosed on these isopycnals for this region is especially sensitive when random error is added to b in the system 𝗔x = b (section 4b).
We define a streamfunction variable at each point around the boundary of the domain on the reference level with the reference level velocity vref(vref = ∇γΨγ0 × k/f ). The streamfunction at the ends of tracer contours that are not at grid points are linearly interpolated to the closest boundary grid points. A simple matrix inversion is applied to solve for the reference level streamfunction on the reference level and the mixing coefficients K and D on each isopycnal with the added constraint that the mean value of Ψγ0 be minimized. This does not minimize the reference level velocity but helps the inversion find a solution for Ψγ0.
The resulting solution for the mixing coefficients is shown in Figs. 6 and 7 (dotted line) and is in agreement with the HIM values (gray solid line). The component of the reference level velocity directed into the domain υref for υref = vref · m, where m is the unit vector perpendicular to the domain boundaries, is also shown (Figs. 8 and 9; dotted line). The reference-level velocity is the same in HIM (gray line in Figs. 8 and 9) to within interpolation error and the definition of the streamfunction.
c. The thickness-weighted mean box inverse method
The principal equations employed in a box inverse method are for the thermal wind and conservation of volume, temperature, and salinity. These are commonly represented in a “layered” framework (neutral or otherwise), assuming steady state and allowing for vertical diffusion D and vertical advection w* or the net diapycnal flux FT of property T. Such methods have been shown to accurately describe the circulation when the row and column weights of the underdetermined problem and the reference level are chosen appropriately (McIntosh and Rintoul 1997). This was done by taking a snapshot of an eddy-permitting ocean model. Model drift or “storage” terms were corrected for in McIntosh and Rintoul as is done here. We formulate the box equations such that the full thickness-weighted mean isopycnal flow is considered; that is, we include the sum of the mean velocity v plus the bolus velocity v* [parameterized as in (5)].
Sloyan and Rintoul (2000) showed that, while net vertical transports may be resolved accurately, diagnosing a diapycnal mixing coefficient is problematic. This is perhaps due to correlations between the vertical velocity and tracer values on the isopycnal faces, . In McIntosh and Rintoul (1997) and Sloyan and Rintoul (2000), snapshots of numerical models were used, in which case the advection rate of the isopycnals in the vertical may have contributed as an error term. The inferred vertical velocity from the inversion, w*, has both a diapycnal component wγ (i.e., through the isopycnals) and a component due to the movement of isopycnals ∂p/∂t|γ. If a mixing coefficient is inferred from the combined velocity w*, an error results. As the HIM data used here are 20-yr means and the steady advective–diffusive equations are constrained to hold exactly, no such error may appear.
We reformulate the box inverse equations somewhat to describe a statistically steady state, assuming the use of temporally averaged hydrographic data (either from an atlas or repeat hydrographic sections). The equations representing the conservation of volume and the tracer C are
Here ∫ … dA is the integral over the isopycnal interface; again, m is directed out of the box, and is the integral over the vertical face of the perimeter of the density layer with thickness h. The conservation equation (16) is written for the general tracer C. Here we will use anomalies to the tracers Θ and S such that C = Θ − Θ̃ and C = S − S̃ in which the tilde represents the mean over the layer.
In (15) we are considering a layer of finite thickness h between isopycnals. To be precise, the layers defined in (15) and (16) are finite volumes “about” the isopycnals on which tracer contours are defined in (9); that is, the tracer contours are written in advective from, whereas the box model equations are in divergence form. In practice, however, h is the same in (1) as in (15) and (16); that is, the layers considered in this section are essentially the same as the isopycnals considered in section 3b, and we will continue to describe them as isopycnals.
If a linear equation of state is assumed, the along-isopycnal diffusion term, , is zero, and, neglecting vertical variations in D, (17) reduces to wγ = Dγzz /γz, as commonly used in studies attempting to infer large-scale dynamics from observations (Karsten and Marshall 2002; Naveira-Garabato et al. 2007). Here we shall retain both the linear and nonlinear components and the vertical gradient of D [see appendix B for a full derivation of (17)].
As in section 3b, we solve for a geostrophic streamfunction on the reference level and a single value for K and D on each isopycnal including one D value above and below the top and bottom isopycnals, respectively. This is so that Dzz may be represented on each isopycnal where the box equations are applied. Later, when tracer-contour and box equations are combined, tracer contours are defined on these upper and lower isopycnals, as Dz, and Dzz do not enter the tracer-contour equations (12).
As the system is underdetermined, an infinite number of solutions exist. To find a unique solution, we impose a weak “level of no motion” constraint. This is done by adding a set of equations that define the reference level velocity to be zero and impose a very small weight on these equations such that a large error is tolerated on them relative to (15) and (16). In practice, the relative weighting used is 10−5 [i.e. when all the terms in 𝗔 and x are weighted to be O(1), the reference level minimization equations are weighted to be O(10−5)]. Changing the relative weight by plus or minus an order of magnitude does not change the resulting solution (not shown).
In the North Pacific region the reference level is chosen as σ2 = 36.89 kg m−3 at a depth of ∼2000 m. Given that the reference level velocity in this region is very small, O(10−5 m s−1), the box inverse method gives a good representation of the vertical and along-isopycnal mixing coefficients (Fig. 6), diverging in the case of K by ±10% or so, on the shallowest isopycnals (Fig. 6). The box inverse method reveals some information about the reference level velocity (Fig. 8).
In the South Atlantic region the reference level is σ2 = 36.96 kg m−3 at a depth of ∼2600 m. In this region depth integrated velocities are large. Using the box inverse method, the reference level velocity and mixing coefficients are much less accurately described by the box model (Fig. 7). The solution would be improved, given a more appropriate a priori value for the reference level velocity.
It is likely that, if more isopycnals were included in the box method, the solution for the reference level velocity would improve. McIntosh and Rintoul (1997) were able to resolve the reference level velocity using approximately 25 levels. They showed that, given a sufficient number of levels, no equation error, and known or zero diapycnal fluxes, in principle, a box inversion can become overdetermined and an exact solution can be found. However, they found that increasing the number of isopycnals makes the solution far more sensitive to errors.
Here, boxes are bound by “sections,” which start and end within an ocean basin, away from continental boundaries. Both McIntosh and Rintoul (1997) and Sloyan and Rintoul (2000) test the box inverse method on boxes bounded by transcontinental sections, which may constrain their solution better. Incorporating the tracer-contour equations into basin-scale inversions is left to future work. Also, box inversions are seldom used as stand-alone tools without a priori oceanographic knowledge included as added constraints. Such constraints may include the net transport across a transcontinental section, velocity from floats or moorings, a reference level minimization based on observations of water masses, or an assumption of the magnitude diapycnal mixing. The box inverse method offers a practical way of including such varied information in a statistically consistent way while constraining the large-scale flow to be mass, heat, and salt conserving. These are all good reasons for combining the box and tracer-contour equations in a flexible manner.
Given that both the box equations and tracer-contour equations are written in terms of the same unknowns, it is trivial to combine the two methods. Solving the full set of equations (without a reference-level minimization) the solution is again well defined to within interpolation error (Figs. 6 and 7). When the system is overdetermined, the reference-level minimization is not needed.
d. The Bernoulli streamfunction inverse method
First proposed by Killworth (1986), the Bernoulli inverse method involves following contours of constant potential vorticity (PV) on isopycnals. In this method, it is assumed that such a contour is a streamline of the geostrophic flow and a reference level streamfunction is diagnosed using the thermal wind. In a more recent application of this method, Cunningham (2000) assumes that streamlines follow contours of constant temperature and salinity. The method is then identical to using the tracer-contour equations for the case in which K and D are assumed to be zero. If the assumption of zero along-isopycnal and vertical mixing is correct, the Bernoulli method should resolve the reference level velocity with the same accuracy as the tracer-contour equations (section 3b). It is clear, however, that the method will give no information about the mixing coefficients, K and D, nor about the downgradient fluxes of tracers, such as temperature and salinity, since these are explicitly minimized in the Bernoulli inverse method.
Given output of HIM, we apply a form of the Bernoulli method such that [from (11)]
In (19), the cross-contour flow ΔΨγ is minimized. The only unknown is the reference level streamfunction Ψγ0.
We apply (19) to the North Pacific and South Atlantic regions, as described in the previous sections. The reference level velocity υref diagnosed using this method is not in good agreement for the North Pacific region. This is because, in the North Pacific region, the velocity on the reference level, the unknown, is small [O(0.1 mm s−1)], whereas the cross-contour velocity on shallower isopycnals [assumed to be zero in (19)] is large [O(5 mm s−1)]. The Bernoulli streamfunction method fares much better in the South Atlantic region where the reference level velocity is large [O(1 cm s−1)] and the flow is largely along the path of the temperature contours of the ACC.
The Bernoulli method may be an appropriate method for diagnosing the along-contour component of the absolute velocity v · (n × k). A recent proposal by Alderson and Killworth (2005) to use Argo data in real time may be one such example. However, no such method can be applied to time-mean hydrographic data to infer any information about the downgradient transports of tracers, because these are explicitly minimized.
4. Effect of random error
The results presented in section 3 leave open some critical questions: what is the sensitivity of the new equations to random error and how should such equations be weighted relative to one another?
a. Weighting the different equations and unknowns
The condition number of a matrix is ||𝗔||∞||𝗔−1||∞, where ||𝗔||∞ is the infinity norm of the matrix 𝗔. The condition number is a key indicator of a matrix’s sensitivity to error (Turing 1948; Tarantola 1987). Following McIntosh and Rintoul (1997), weights are chosen in the tracer-contour method such that the condition number of the matrix 𝗔 is minimized. McIntosh and Rintoul adjust column weights such that the inversion is penalized for error in either the vertical velocity or the reference velocity. As the inclusion of tracer contours makes the inversion overdetermined, column weights have no effect and thus only row weights are important.
We must consider the relative weights of the two types of equations: box property conservation equations (15) and (16) (as well as reference level minimization) and contour equations (12). In practice, the reference level minimization only has a significant effect when its relative weight is equivalent to the weights on the tracer-contour equations. The ratio of box to tracer-contour weight is adjusted to give a minimum condition number. Weights are found to be within an order of magnitude of the initial scaling (K and D coefficients are initially scaled by 103 and 10−5 m2 s−1, respectively, and each row is scaled by its mean value). In the analysis of this section, there is no reference level minimization. Inclusion of a reference level minimization, as in section 3c, makes little or no difference to the result, since the tracer-contour equations are always included and make the inversion overdetermined.
b. Propagation of random error
McIntosh and Rintoul (1997) find that, in cases where interfacial fluxes are known or not present, adding a sufficient number of isopycnals to a box inversion can make the inversion overdetermined and an exact solution can be found. However, they find that adding random errors of ±10% to each equation gives a noisier solution than in the case of less isopycnals.
To test whether adding tracer contours has a similar effect to adding more isopycnals to the inversion, we add 10% randomly distributed errors to b (σ = ±0.1b). This is repeated 100 times, each time with new errors, to find the ensemble average solution and the normal distribution of the solution. We use the full tracer-contour method, with box and tracer-contour equations, optimally weighted to give the lowest condition number. The solution for D and K varies by only 10%–20% in both the North Pacific and South Atlantic regions (Figs. 10 and 11). In the North Pacific region the ensemble mean solution for the reference level is in good agreement with the HIM velocity. However, the magnitude of the error is large relative to the magnitude of the HIM value (Fig. 12). It should be noted, however, that these are small velocities, O(0.1 mm s−1). In the South Atlantic region the full inversion with 10% errors fares far better than the conventional box inverse method, and the reference level velocities are well resolved (Fig. 13).
Here we have added only random independent errors to b. It is possible that the use of a gridded climatology from limited hydrographic observations, may introduce correlated errors. It is uncertain what effect such errors would have on the results shown here. The sensitivity of the method to the way data are gridded and averaged will be the subject of future work.
The terms in (2) are also possible sources of correlated errors. The effect of the unsteady term, Θt|γ /∇γΘ, may be tested using a linear trend, say, over a 50-yr period. Preliminary work has suggested that this term is small, especially in the ocean below the thermocline where changes on isopycnals occur very slowly. The ∇γK term in (2) is likely to be small, distant from the boundary between different mixing regimes: that is, as long as the region of the inverse method is not on either side of a strong and permanent mixing barrier, such as a strong front. The effect of a strong gradient of K could be taken into account, given knowledge of its spatial variability. For example, consider the unknown K on a given isopycnal. One may assume that K has the same spatial variations as a lateral mixing coefficient, observed at the surface, but have a different magnitude. One could then assume K(γ, x, y) = K1(γ)Kobs(x, y), where K1(γ) is an unknown function of density (or depth) and Kobs(x, y) is an assumed function of latitude and longitude. The sensitivity of this method to such choices will be the subject of future work.
c. Value of the box equations
Using the conventional box equations it is not possible to make the inversion overdetermined when interfacial fluxes are included. We have demonstrated in section 3b that using the tracer-contour equation (12) we can develop an overdetermined system that solves for a geostrophic streamfunction Ψγ and the along-isopycnal and vertical mixing coefficients. Here, we demonstrate that there is some skill gained by combining both the box equations [(15) and (16)] with the tracer-contour equation (12) to form the tracer-contour inverse method, as is done in section 3c.
We may plot the condition number as a function of the relative row weight: that is, the relative weight of the tracer-contour equations to box equations (Fig. 14). For low box-to-contour weight (i.e., strongest influence of the tracer-contour equations), the condition number plateaus at O(103). At the other extreme (strongest influence of the box equations), the condition number rises linearly. Notably, near a relative contour to box weighting of 10−1, a dip in the condition number exists, suggesting that the box equations add conditionality to the system. Of course, the condition number does not necessarily correspond to the error in the inverse estimate.
To see what effect the relative weight of each equation type has on the error in the solution for K, D, and υref, we run an ensemble of inversions, as in the previous section, for each relative weighting. We determine the average error for K, D, and υref by subtracting the inverse result from the HIM value for each ensemble run and taking the average over the 100 inversions (Fig. 14). The results clearly show an increase in error for greater weighting on the box inversion equations. For greater weight on the tracer-contour equations, the error plateaus one to three orders of magnitude lower than at the other extreme. There is a slight dip in the error at close to a relative weighting of 10−1. This is especially the case for K. The only case where there is no change in the error with weight is for D in the North Pacific region as the solution is close to perfect at both limits. The minimum error corresponds to the minimum condition number to within an order of magnitude.
The addition of box equations does improve the result of the inversion. The box equations constrain the flow to conserve volume, heat, and salt on isopycnals, whereas the tracer-contour equation does not (although it is derived from those conservation equations). Further reasons for retaining the box equations are given in section 6.
5. How does the tracer-contour method work?
This section shows how the tracer-contour method is able to solve for the mixing coefficients and the absolute velocity. It is shown that, to solve for the mixing coefficients K and D and the absolute velocity, all inverse methods require either nonzero neutral helicity H n or that the direction of the along-isopycnal tracer gradient spirals in the vertical. By neutral helicity, we are not referring to diapycnal upwelling associated nonlinear process but simply
a. Mixing, spiraling, and neutral helicity
For convenience, we split the lateral velocity v into a component perpendicular to tracer contours on isopycnals (i.e., the cross-contour velocity) and a component parallel to tracer-contours on isopycnals υ|| (i.e., the along-contour velocity) such that
Here we consider the along-isopycnal temperature and salinity gradient, ∇γΘ and ∇γS, such that n = ∇γΘ/|∇γΘ| ≡ ∇γS/|∇γS|. However, the analysis here may be generalized for any tracer by replacing Θ with C [see appendix C for a derivation of the analogous form of (7)]. In this section, we consider Θ (which is equivalent to S on an isopycnal) as a connection may be drawn to neutral helicity.
The tracer-contour method considers the component of the lateral velocity, down the along-isopycnal tracer gradient, , given by (7). In steady state, is directly related to mixing (7), whereas υ|| has no direct implication for property conservation or mixing. In this sense, υ|| is truly an “adiabatic” velocity component. In a steady ocean without along-isopycnal or vertical mixing, flow would exactly follow contours of constant tracer on isopycnals.
In the tracer-contour method, velocity is related in the vertical through (11), which, when reduced to a point (laterally and vertically), is simply the thermal wind. Here we write the thermal wind equation, at a point, in terms of the along-isopycnal gradient of pressure p such that
As we are interested in the downgradient component of the flow, , we consider (23) in the downgradient direction n, such that
The rhs of (24) may be expanded using the chain rule into a component involving the vertical derivative of and another involving the vertical derivative of the unit vector n. Hence,
Considering each term in (25), we will show that it describes a balance between neutral helicity, mixing, and spiraling of the along-isopycnal temperature gradient. We will then discuss the physical meaning of each term.
We note that H n is proportional to the gradient of pressure along a temperature contour on an isopycnal (20). Neutral helicity is thus proportional to the component of thermal wind in the cross-contour direction. The lhs of (25) is then
where ϕ is the direction of n, expressed in radians, and ϕz is the rate of spiraling of the isopycnal temperature gradient in the vertical. So, the second term on the rhs of (25) becomes
From (29), we infer that, so as to diagnose mixing, and the along-contour flow υ||, the technique requires that there be some component of the thermal wind vector down the along-isopycnal temperature gradient or that there be spiraling of this gradient. In the first case, there must be a pressure change along a tracer contour on an isopycnal and hence that neutral helicity must be nonzero. Neutral helicity and spiraling give information on the vertical derivative of the cross-contour flow . Equation (29) suggests that the tracer-contour inverse method will be least sensitive to error in determining , and hence K and D, in regions where neutral helicity and spiraling are strong.
b. The absolute velocity and the beta-spiral and Bernoulli methods
Both Needler (1985) and McDougall (1995) follow a similar approach to that presented here, but in terms of PV gradients rather than temperature and salinity gradients. McDougall shows that the beta-spiral method can be collapsed to a point and the absolute velocity can be presented in terms of mixing and the spiraling of the along-isopycnal PV gradient. The Bernoulli method, considered in terms of PV, requires the same spiraling and along-contour pressure gradient as in the beta-spiral method. Using Θ, both the beta-spiral and Bernoulli methods are described by Eq. (29).
Determination of the absolute velocity vector has been the motivation of a long history of work (Stommel and Schott 1977; Needler 1985; Killworth 1986; McDougall 1995). We now develop a new expression for the absolute velocity dependent on neutral helicity, mixing, and spiraling.
Since some degree of spiraling of temperature contours is ubiquitous in the global ocean, (30) is a closed expression for the absolute velocity. The velocity is expressed in terms of H n and mixing, .
Assuming both steady-state and zero mixing (i.e., ), the second and third terms in (30) vanish, in which case,
Equation (31) is a closed expression for the absolute velocity in a mixing-free ocean, where ∇γΘ and ϕz are nonzero—that is, where an along-isopycnal temperature gradient exists and its direction is not uniform in the vertical. The mean velocity would be solely in the along-contour direction k × n with magnitude |v| = |(−N 2/fρgϕz)∇γp × n|. Expression (31) is far less differentiated than those recognized by Needler (1985) or McDougall (1995).
The Bernoulli and beta-spiral methods use the thermal wind equation (23). In their classical formulation, these methods assume no mixing and that the geostrophic flow is along tracer contours (PV or Θ) on isopycnals. This is imposed explicitly in the Bernoulli method and implicitly in the beta-spiral method. Then, for each method, determining the absolute velocity depends on two things: the component of thermal wind down the tracer gradient (neutral helicity) and the rate of spiraling of the tracer gradient (ϕz).
The box model is a special case in that it uses thermal wind and tracer conservation but also has an explicit reference level minimization (i.e., the velocity at some deep level is minimized to an a priori value). This reference level minimization sets the magnitude of v throughout the water column. In this sense, even in the absence of neutral helicity, the along-contour velocity υ|| and the cross-contour velocity may be determined [i.e. the rhs of (29)]. With knowledge of v, the mixing coefficients may also be inferred through property conservation on each isopycnal, (15) and (16).
The dynamical importance of neutral helicity has been associated with the nonlinearity of the equation of state (McDougall and Jackett 1988, 2007) and has implications for diapycnal upwelling (McDougall and Jackett 1988; Klocker and McDougall 2010). Here we have shown its critical role in determining the mean circulation in the ocean, (30). A connection between these two findings has not yet been identified and is left to future work.
6. Discussion and ongoing work
An inverse method has been presented for diagnosing rates of along-isopycnal and vertical mixing and the mean geostrophic velocity from time-averaged hydrographic data. The method involves integrating the “curvatures” of hydrographic data along tracer contours on isopycnals. These curvatures form the coefficients of the diffusivities, D and K, and the sum of these are related to differences in a geostrophic streamfunction along tracer contours. The streamfunction on that isopycnal is then related to a reference level streamfunction, on either a reference isopycnal or on a pressure surface. Assuming that the flow is steady state and that the imposition of a spatial structure on the mixing coefficients is reasonable (e.g., constant on each isopycnal), an overdetermined system is developed, which accurately resolves the mixing coefficients and the mean flow.
The tracer-contour inverse method combines physically consistent equations including conservation of mass, heat, and salt and geostrophic balance. Importantly, it does not require a reference-level velocity minimization. Mixing coefficients are leading order in the inversion, as is the component of the large-scale residual flow down the along-isopycnal tracer gradient. Where the tracer is conservative temperature (equivalently salinity) on an isopycnal, this component of the velocity is effectively the along-isopycnal component of the thermohaline circulation. Tracer contours are combined with the large-scale property budgets of the box inverse method, a framework in which a priori information such as section transports and observed mixing parameters may be included to constrain the results. The influence of the box equations versus the tracer contours can be easily adjusted. Results presented here suggest that the inclusion of tracer contours in large-scale inversions will allow a significant improvement in our understanding of diffusive processes and the general circulation. The method is presently being developed for more general application to time-averaged hydrographic data, including the treatment of coastlines, the mixed and Ekman layers, and multiple boxes.
We thank Dr. Peter McIntosh for his helpful feedback on this work and thank Drs. Carl Wunsch and Andreas Thurnherr for their insightful comments on an earlier version of the manuscript. The authors thank Dr. Robert Hallberg for providing the numerical model output. This work contributes to the CSIRO Climate Change Research Program and has been partially supported by the CSIRO Wealth from Oceans Flagship and the Australian Research Council Centre of Excellence for Mathematics and Statistics of Complex Systems.
Derivation of the Water Mass Equation
The conservation equations for salinity S and conservative temperature Θ are
These equations have been written in the advective form and with respect to neutral density γ n so that wγ is the vertical velocity through isopycnals (i.e., the diapycnal velocity), v is the horizontal velocity, and (/h) is the thickness-weighted horizontal velocity. The thickness is h = Δρ/ρz for some small Δρ. The conservative temperature Θ and salinity S in (A1) and (A2) are the thickness-weighted values. The mixing processes that appear on the right-hand sides are simply along-isopycnal mixing of passive tracers (with coefficient K) along the density surfaces and vertical small-scale turbulent mixing (with coefficient, D). We have not included double-diffusive convection or double-diffusive interleaving. We define the cross-contour direction as n = ∇γS/|∇γS|. Notably, ∇γS/|∇γS| ≡ ∇γΘ/|∇γΘ| as we are considering gradients on a neutral density surface. Multiplying (A1) by Θz and (A2) by Sz gives
Above, we noted that α/β = |∇γS|/|∇γΘ| = (St/Θt)|γ, simplifying the unsteady term in (A5).
The stability ratio is Rρ = αΘz /βSz so that Rρ ≡ |∇γS|Θz /|∇γΘ|Sz. Hence, neither nor 1/λ γ contain singularities if the water column is stably stratified (i.e., if Rρ ≠ 1 ⇒ Θz|∇γS| − Sz|∇γΘ| ≠ 0).
We may remove a thickness gradient term from such that
Note that here ∇γ log(h) · n does not represent the bolus transport. In the case of a linear equation of state ∇γα ≡ ∇γβ ≡ 0, hence simplifies to
Given the equality
where |x,y signifies constant latitude and longitude, 1/λ γ may be written
The thickness-weighted mean velocity /h may be represented as
An additional mixing term, KPV(β/f )j, may be added that relates the bolus velocity v* to gradients of PV rather than thickness h (McDougall 1991).
The cross-contour velocity v · n may now be represented as
where 1/λh = ∇γ log(h) · n.
Derivation of the Density Equation
Again, neither nor 1/ηγ are singular if the water column is stably stratified.
If a linear equation of state is assumed such that ρz ∝ −αΘz + βSz, for α and β constant, and again α/β = |∇γS|/|∇γΘ|, the vertical mixing term 1/ηγ becomes
Also, since ∇γα ≡ ∇γβ ≡ 0, the along-isopycnal mixing term becomes
Hence, for a linear equation of state, the vertical velocity through isopycnals simplifies to wγ = Dρzz /ρz + Dz.
Derivation of an Isopycnal Advective–Diffusive Balance Equation for a Conservative Tracer C
We will now derive an equation that relates advection down the isopycnal gradient of a conservative tracer C to mixing. Consider the steady-state conservation equation for C in an incompressible ocean, distant from sources and sinks:
Clearly, the last terms on the lhs and rhs of (C2) cancel. Expanding (/h) using (A10) grouping the along-isopycnal and diapycnal mixing terms, and dividing through by the absolute tracer gradient |∇γC|, (C2) may thus be written
nC is the unit vector in the direction of the along-isopycnal gradient of C (i.e., nC = ∇γC/|∇γC|).
Assuming a linear equation of state, and . Hence,
Given (A8), the vertical mixing scale length, assuming a linear equation of state, may be written as
So, the isopycnal advective–diffusive balance equation for the tracer C is analogous to the water mass equation (1) or (7). As such, along-isopycnal mixing enters through the ratio of along-isopycnal curvature to along-isopycnal gradient ∇γ2C/|∇γC|, and vertical mixing largely enters through the curvature of the tracer with respect to density d2C/dρ2 (note: if a linear equation of state is assumed, then d2Θ/dS2 can be written in terms of the curvature of either Θ or S with respect to ρ).
Reference Level Streamfunction at a Fixed Pressure
Here we formulate the difference in geostrophic streamfunction on neutral density surfaces ΔΨγ, relative to a difference on a pressure surface. When solving for a streamfunction near topography or outcropping regions, it is cumbersome to use a density surface that is nonexistent at some latitudes and longitudes. As discussed in McDougall and Klocker (2009), errors arise in defining a streamfunction on a density surface over basin to global scales. Although these errors can be reduced by considering appropriate specific volume anomaly and correction terms, it is likely best to reference to a constant pressure, at or near mean sea level.
Choosing the reference level to be a pressure surface, (11) simplifies to
For each equation of type (D1), the choice of Θconst and Sconst and hence the nature of δ is arbitrary. For a specific tracer contour of temperature Θcontour and salinity Scontour, a specific volume anomaly variable δcontour may be defined as
Along such a contour, δcontour = 0; hence, (D1) simplifies further to
Equation (D3) offers a possible increase in accuracy but requires the redefinition of a specific volume anomaly variable δcontour for each individual tracer contour.
Defining Tracer Contours
Here we describe how tracer contours have been defined from fields of the Hallberg isopycnal model (HIM). The North Pacific region (10°–30°N, 210°–230°E) is used as an example.
We start with three-dimensional fields of θ, S, p, and h and the length scales , λ γ, and λh. These variables are defined between the latitudes and longitudes of the region and for all definable isopycnals. The data are spaced by approximately 1° of latitude and longitude. We do not consider isopycnals that interact with the ocean floor or those completely above 500-m depth. The North Pacific region has 19 longitudinal points by 21 latitudinal points by 15 isopycnals (N = 15).
From each point on the boundary, on each isopycnal a temperature contour is followed until it exits the domain. Doing this for the North Pacific region results in 78 contours on each isopycnal with a total of 1170 contours (i.e., 78 × N). Contours are removed if they are less than 2° long, have an average pressure less than 500 db, and/or have a minimum temperature gradient less than 108 K m−1 along them.
In the inversion, contour ends are linearly interpolated onto a set of points around the domain. A streamfunction is solved for at these points on a reference level. If no contour starts or ends close to a given point on the boundary, on any level in the vertical, no solution can be found for the streamfunction difference at that point. An example of this is the northwest corner of the North Pacific region. Since the temperature gradient is to the southwest on each isopycnal, no contour starts or ends near the northwest corner; the corner point is thus isolated from the inversion.
To avoid this problem, we solve for the streamfunction on a set of points around the domain such that each interacts with at least 10 contours (10 is an arbitrary choice). In practice, this is done by moving around the boundary points removing each point that is associated with less than 10 contours. The contours that the original point was associated with are then assigned to the remaining points, and the procedure is repeated. For the North Pacific region, 56 boundary grid points are defined (L = 56). The streamfunction is solved on these points.
The data are at 1° resolution. Thus, tracer contours should be spaced by about 1°. Since the contours are defined at every point on the boundary of the domain, the spacing of tracer contours is about 0.5°. Hence, we expect the effective number of contours to be half those defined: that is, half the number of boundary grid points (28 for the North Pacific region).
Corresponding author address: Jan Zika, Laboratoire des écoulements géophysiques et industriels, BP 53, 38041, Grenoble, CEDEX 9, France. Email: firstname.lastname@example.org