The exact geostrophic streamfunction for neutral surfaces

TL;DR

This paper derives the exact geostrophic streamfunction for neutral surfaces using a path integral approach, introduces a practical approximation called the topobaric geostrophic streamfunction, and demonstrates its superior accuracy in estimating oceanic geostrophic velocities.

Contribution

The paper provides the first explicit form of the exact geostrophic streamfunction for neutral surfaces and introduces the topobaric approximation for practical use in oceanography.

Findings

Topobaric geostrophic streamfunction reduces velocity estimation error by an order of magnitude.

The exact streamfunction involves a path integral over a neutral trajectory and is multivalued due to surface holes.

The orthobaric Montgomery potential offers a computationally efficient approximation with improved accuracy.

Abstract

McDougall (1989) proved that neutral surfaces possess an exact geostrophic streamfunction, but its form has remained unknown. The exact geostrophic streamfunction for neutral surfaces is derived here. It involves a path integral of the specific volume along a neutral trajectory. On a neutral surface, the specific volume is a multivalued function of the pressure on the surface, $\underset{\sim}{p}$. By decomposing the neutral surface into regions where the specific volume is a single-valued function of $\underset{\sim}{p}$, the path integral is simply a sum of integrals of these single-valued functions. The regions are determined by the Reeb graph of $\underset{\sim}{p}$, and the neutral trajectory is encoded by a walk on this graph. Islands, and other holes in the neutral surface, can create cycles in the Reeb graph, causing the exact geostrophic streamfunction on a neutral surface to be multivalued. Moreover, neutral surfaces are ill-defined in the real ocean. Hence, the topobaric geostrophic streamfunction is presented: a single-valued approximation of the exact geostrophic streamfunction for neutral surfaces, for use on any well-defined, approximately neutral surface. Numerical tests on several approximately neutral surfaces reveal that the topobaric geostrophic streamfunction estimates the geostrophic velocity with an error that is about an order of magnitude smaller than that for any previously known geostrophic streamfunction. Also, the Montgomery potential is generalized, leading to an alternate form of the exact geostrophic streamfunction for neutral surfaces. This form is used to construct the orthobaric Montgomery potential, an easily computable geostrophic streamfunction that estimates the geostrophic velocity more accurately than any previously known geostrophic streamfunction, but less so than the topobaric geostrophic streamfunction.

Figures

Tables

Table 1. Table 1: The ℓ 1 subscript ℓ 1 \ell_{1} (a) and ℓ 2 subscript ℓ 2 \ell_{2} (b) norms of the error ϵ italic-ϵ \epsilon [ m ​ s − 1 m superscript s 1 \mathrm{m}\,\mathrm{s}^{-1} ] in estimating the “true” geostrophic velocity, for various velocity estimates (rows) and surfaces (columns). In the top row, the velocity estimate is the full velocity, so ϵ italic-ϵ \epsilon is the ageostrophic velocity.

(a)	$\underset{~}{z}\left({𝒙}_0\right)=-997.99m$			$\underset{~}{z}\left({𝒙}_0\right)=-1988.60m$
${\Vert ϵ\Vert }_1$	${\sigma }_1$	${\tau }^{\prime }$	$\omega$	${\sigma }_2$	${\tau }^{\prime }$	$\omega$
$\underset{~}{𝐮}$	6.0e-03	6.0e-03	5.9e-03	4.4e-03	4.4e-03	4.3e-03
${\underset{~}{{𝐮}_g}}^{\left(ZH\right)}$	1.5e-03	4.0e-04	3.3e-04	3.6e-03	3.7e-04	3.0e-04
${\underset{~}{{𝐮}_g}}^{\left(Cu\right)}$	6.2e-04	6.0e-04	4.0e-04	1.7e-03	1.3e-03	7.9e-04
${\underset{~}{{𝐮}_g}}^{\left(MK\right)}$	5.1e-04	1.4e-04	9.3e-05	1.4e-03	2.3e-04	1.4e-04
${\underset{~}{{𝐮}_g}}^{\left(oM\right)}$	2.3e-04	9.3e-05	8.9e-05	2.1e-04	3.1e-05	6.4e-05
${\underset{~}{{𝐮}_g}}^{\left(tb\right)}$	3.8e-05	1.0e-06	1.7e-05	5.9e-05	1.3e-06	1.4e-05
${\underset{~}{{𝐮}_g}}^{\left({\nabla }_{z}p\right)}$	1.7e-05	1.4e-05	1.2e-05	2.0e-05	1.0e-05	1.0e-05

Table 2. Table 2: The percentage of the surface area, on different surfaces (columns), for which the error in estimating the “true” geostrophic velocity by methods (b)–(g) (rows) exceeds the ageostrophic velocity, in magnitude.

	$\underset{~}{z}\left({𝒙}_0\right)=-997.99m$			$\underset{~}{z}\left({𝒙}_0\right)=-1988.60m$
	${\sigma }_1$	${\tau }^{\prime }$	$\omega$	${\sigma }_2$	${\tau }^{\prime }$	$\omega$
${\underset{~}{{𝐮}_g}}^{\left(ZH\right)}$	6.90	3.30	2.50	7.91	5.05	4.25
${\underset{~}{{𝐮}_g}}^{\left(Cu\right)}$	3.41	3.39	2.04	4.31	4.84	3.21
${\underset{~}{{𝐮}_g}}^{\left(MK\right)}$	3.07	1.32	0.71	4.24	2.34	1.40
${\underset{~}{{𝐮}_g}}^{\left(oM\right)}$	2.00	1.08	1.09	3.11	0.48	1.02
${\underset{~}{{𝐮}_g}}^{\left(tb\right)}$	0.28	0.00	0.14	0.47	0.01	0.15
${\underset{~}{{𝐮}_g}}^{\left({\nabla }_{z}p\right)}$	0.12	0.11	0.10	0.19	0.15	0.15

Equations

$$-\hat{\mathbf{k}}\times f{\mathbf{u}}_g=\nu {\nabla }_{z}p,$$

$${\nabla }_r\Psi =\nu {\nabla }_{z}p$$

$${\nabla }_r\Psi =\nu {\nabla }_{r}p+{\nabla }_r(gz).$$

$$\nabla \underset{˜}{\Psi }=\underset{˜}{\nu }\nabla \underset{˜}{p}+\nabla (g\underset{˜}{z}).$$

$$\nabla \underset{˜}{\nu }=\underset{˜}{{\nu }_p}\nabla \underset{˜}{p}.$$

$$0=\underset{˜}{{\nu }_S}\nabla \underset{˜}{S}+\underset{˜}{{\nu }_{\theta }}\nabla \underset{˜}{\theta },$$

$$\underset{˜}{\nu }=\hat{\nu }(\underset{˜}{p})$$

$$\boxed{{\underset{˜}{\Psi }}^{(n)}=g\underset{˜}{z}+{\int }_{p_0}^{\underset{˜}{p}}\hat{\nu }(p^{\prime })\mathrm{d}{p}^{\prime }}$$

$$\underset{˜}{z}=\eta -\frac{1}{g}{\int }_{p^{(\eta )}}^{\underset{˜}{p}}\nu \mathrm{d}p,$$

$${\underset{˜}{\Psi }}^{(p)}=g\eta -{\int }_{p^{(\eta )}}^{\underset{˜}{p}}\nu \mathrm{d}p$$

$${\underset{˜}{\Psi }}^{(n)}(\mathbf{x})={\underset{˜}{\Psi }}^{(n)}({\mathbf{x}}_0)+{\int }_{\mathcal{P}}\nabla {\underset{˜}{\Psi }}^{(n)}\cdot \mathrm{d}\mathbf{r}.$$

$${\underset{˜}{\Psi }}^{(n)}=g\underset{˜}{z}+{\int }_{\mathcal{P}}\hat{\nu }(\underset{˜}{p})\mathrm{d}p.$$

$${\underset{˜}{\Psi }}^{(n)}=g\underset{˜}{z}$$

$$+$$

$$0=\sum _{j=1}^{\mathcal{J}-1}{\int }_{p_{n_j}}^{p_{n_{j+1}}}{\hat{\nu }}_{a_j}(p^{\prime })\mathrm{d}{p}^{\prime }.$$

$$\displaystyle 0=\sum_{j=1}^{\mathcal{J}-1}\Biggl{(}$$

$$\displaystyle+\int_{p_{n_{j}}}^{p_{n_{j+1}}}\int_{p_{n_{j}}}^{p^{\prime}}\hat{\pi}_{a_{j}}(p^{\prime\prime})\,\mathrm{d}p^{\prime\prime}\,\mathrm{d}p^{\prime}\Biggr{)},$$

$${\hat{\nu }}_a(p^{\prime })=K_a+L_a(p^{\prime }-p_{{\ell }_a})+\mathcal{V}(S_0,{\theta }_0,p^{\prime }),$$

$$\delta =\nu -{\nu }_0,$$

$$\nabla {\underset{˜}{\Psi }}^{(\delta )}=\nabla \left(\underset{˜}{\delta }\underset{˜}{p}+{\int }_P^{\underset{˜}{p}}{\nu }_0\mathrm{d}p+g\underset{˜}{z}\right)-\underset{˜}{p}\nabla \underset{˜}{\delta },$$

$${\int }_P^{\underset{˜}{p}}{\nu }_0\mathrm{d}p+g\underset{˜}{z}=-{\int }_{p^{(\eta )}}^{\underset{˜}{p}}\delta \mathrm{d}p+g\eta .$$

$$\displaystyle\nabla\utilde{\Psi}^{(ZH)}=\$$

$$-(\underset{˜}{p}-p_0)\nabla \underset{˜}{\delta },$$

$$\displaystyle\utilde{\Psi}^{(MK)}=\$$

$$-\frac{1}{12}\frac{T_b}{\underset{˜}{\rho }}(\underset{˜}{\theta }-{\theta }_0)(\underset{˜}{p}-p_0{)}^2,$$

$${\underset{˜}{\Psi }}^{(Cu)}={\int }_P^{\underset{˜}{p}}\mathcal{V}(\underset{˜}{S},\underset{˜}{\theta },p^{\prime })\mathrm{d}{p}^{\prime }+g\underset{˜}{z}.$$

$$\nabla {\underset{˜}{\Psi }}^{(Cu)}=\underset{˜}{\nu }\nabla \underset{˜}{p}+\nabla (g\underset{˜}{z})+{\int }_P^{\underset{˜}{p}}{\nabla }_{p^{\prime }}\mathcal{V}(\underset{˜}{S},\underset{˜}{\theta },p^{\prime })\mathrm{d}{p}^{\prime }.$$

$$\displaystyle\nabla\utilde{\Psi}^{(oM)}=\$$

$$\displaystyle-\big{(}\utilde{p}-\hat{p}(\utilde{\delta})\big{)}\nabla\utilde{\delta},$$

$${\underset{˜}{\Psi }}^{(n)}$$

$$=g\underset{˜}{z}$$

$$={\underset{˜}{\Psi }}^{(oM)}$$

$$-\hat{\mathbf{k}}\times f\underset{˜}{{\mathbf{u}}_g}={\rho }_B^{-1}\underset{˜}{{\nabla }_{z}p}={\rho }_B^{-1}\nabla \underset{˜}{p}+g{\rho }_B^{-1}\underset{˜}{\rho }\nabla \underset{˜}{z}=\nabla \underset{˜}{\Psi }.$$

Full text

††journal: Ocean Modelling

1 Introduction

Outside the mixed layer and bottom boundary layer and on scales larger than about $10\text{\,}\mathrm{m}$, the structure of oceanic flow is predominantly two-dimensional, largely confined to isosurfaces of a quasi-conservative density. This enabled Reid and Lynn [1971], for example, to draw conclusions about how the deep tropical oceans are connected to the shallow polar oceans. This method of analysis is even more powerful if, rather than just knowing the flow is within such a surface, we can describe the flow within that surface, for example by means of a geostrophic streamfunction (GSF).

A GSF is the potential, in a specified surface, for the acceleration by the horizontal pressure gradient. (Hence, this quantity is sometimes called an “acceleration potential”; nomenclature is discussed at the end of this section.) A GSF is not, technically, a streamfunction of the geostrophic velocity: none exists, because meridional variation of the Coriolis parameter implies the geostrophic velocity is divergent. Still, the geostrophic velocity (a 2D vector field) is exactly determined by the GSF (a scalar field) and the known Coriolis parameter. Moreover, the geostrophic velocity is everywhere (except the equator) tangent to contours of the GSF, which are therefore streamlines of the geostrophic velocity.

GSFs possess many theoretical graces, not least being their deep connection with potential vorticity. If a GSF exists on isosurfaces of a materially conserved 3D variable, then this variable gives rise to a materially conserved Ertel potential vorticity: when the momentum equations are recast with this variable as the vertical coordinate, the gradient of the GSF exactly represents the horizontal pressure gradient acceleration, and no spurious acceleration terms are produced [de Szoeke, 2000]. Such 3D variables and their GSFs are thus of prime theoretical and practical importance in layered models. Also, because of the divergence-free flow (namely, the geostrophic velocity multiplied by the Coriolis parameter) they provide, GSFs serve as building blocks for many theoretical models [Stommel, 1948, Stommel and Arons, 1959, Welander, 1971, Rhines and Young, 1982, Luyten et al., 1983, Marshall and Radko, 2003] and inverse models [Killworth, 1986, Cunningham, 2000, Zika et al., 2010], and are useful to analyse quantities in recirculating flow, such as streamwise budgets and averages [Gille, 1997, Shuckburgh et al., 2009, Chapman and Sallée, 2017].

Mapping GSFs was first done on specific volume anomaly surfaces [Montgomery, 1938, Reid, 1965], on which the Montgomery [1937] potential is the GSF. Later, potential density became the preferred quasi-conserved density variable, upon whose surfaces many authors analysed the geostrophic circulation by means of the Montgomery potential [e.g. Bower et al., 1985, Lozier et al., 1995, Aksenov et al., 2011]. However, the Montgomery potential is not a GSF for a potential density surface, and indeed none exists[McDougall, 1989]. Its use on a potential density surface can predict a geostrophic velocity that differs substantially from the true geostrophic velocity, as shown by Zhang and Hogg [1992], who then upgraded the Montgomery potential to minimize these errors. We say the Montgomery potential is the exact GSF on specific volume anomaly surfaces, but it is an inexact, or approximate, GSF on potential density surfaces.

Since McDougall [1987] defined neutral surfaces and also highlighted the problems of using potential density far from its reference pressure, another shift has occurred, towards using surfaces that are more closely aligned with the neutral tangent plane — the plane in which fluid parcels can move adiabatically and infinitesimally without experiencing a buoyant restoring force. Neutral trajectories and neutral surfaces are, respectively, paths and surfaces that are everywhere parallel to the local neutral tangent plane [McDougall, 1987]. Unfortunately, non-linearity in the equation of state for seawater makes neutral trajectories path-dependent, and so neutral surfaces are formally ill-defined [McDougall and Jackett, 1988]. However, in practice this path-dependence is small enough [McDougall and Jackett, 1988, 2007] that we can usefully craft approximately neutral surfaces, which are well-defined and are approximately parallel with the neutral tangent plane; examples include neutral density surfaces [Jackett and McDougall, 1997], orthobaric density surfaces [de Szoeke et al., 2000], $\omega$-surfaces [Klocker et al., 2009], and topobaric surfaces [Stanley, 2019, hereafter S19]. Of course, specific volume anomaly surfaces and potential density surfaces are also approximately neutral surfaces, but “approximately” carries heavier emphasis.

An opportunity now presents itself: use the exact GSF for neutral surfaces, rather than the Montgomery potential or its variants, on these approximately neutral surfaces. McDougall [1989] proved that there is an exact GSF for neutral surfaces, at least in a theoretical ocean where neutral surfaces are well-defined. However, the analytic form of this GSF was not given, and has remained elusive ever since.

The key theoretical result of this paper is to derive the exact GSF for neutral surfaces. It is found to also suffer from path-dependency, this time caused by islands and other holes in the neutral surface. Its ill-defined nature is overcome by the key practical result of this paper: the topobaric GSF, which approximates the exact GSF for neutral surfaces but is well-defined, and can be used on any approximately neutral surface. These key results are essentially corollaries to a theory of neutral surfaces developed in a companion paper [S19]. The root of this theory is a multivalued functional relationship between the specific volume and the pressure on neutral surfaces, and its characterisation by the Reeb [1946] graph. The necessary ideas of that theory are briefly described here, but the reader is encouraged to read the companion paper first.

Complimenting these results, the way Zhang and Hogg [1992] upgraded the Montgomery [1937] potential is generalized. This generalization turns out to be none other than the exact GSF for neutral surfaces, but presenting a different functional form. This alternate form is useful in creating an easily calculable yet highly accurate GSF called the orthobaric Montgomery potential — the second key practical result of this paper.

Section 2 reviews GSFs and neutral surfaces, then derives the exact GSF for neutral surfaces. The topobaric GSF is discussed in Section 3. The Montgomery potential and its variants are reviewed in Section 4, then generalized in Section 5 to create the orthobaric Montgomery potential. Section 6 discusses the model data to be used and the Boussinesq approximation. Numerical comparisons of various GSFs are presented in Section 7, before summarizing in Section 8. A provides formulas for various GSFs in a Boussinesq ocean, and B details their numerical computation.

A final note on nomenclature is warranted, a subject of some confusion at least since a joint letter by Wexler and Montgomery [1941] suggested both the names “stream function” and “acceleration potential”, respectively, for the quantity here called a GSF. There are two issues with the name “acceleration potential”. First, it has a specific, older definition in potential flow theory. Second, it neglects other acceleration terms; technically it should be called the “horizontal pressure gradient acceleration potential” in a particular surface. Later, the name “Montgomery potential” and arose, but this is best reserved for the specific function defined by Montgomery [1937] for use in a specific volume anomaly surface. “Montgomery function” has also been used, but this name bears no direct relation to the horizontal pressure gradient acceleration, the geostrophic velocity, or surfaces other than specific volume anomaly surfaces. This manuscript adopts the now widely-used modification of Wexler’s suggestion: “geostrophic streamfunction”, to be thought of as a compound term with a specific definition — not as shorthand for “streamfunction of the geostrophic velocity”, as no such streamfunction exists.

2 The exact geostrophic streamfunction on a neutral surface

2.1 Preliminary definitions

The salinity, potential temperature, pressure, in-situ density, and specific volume are 3D scalar fields denoted by $S$, $\theta$, $p$, $\rho$, and $\nu$ respectively.222The theory presented here could equivalently use Absolute Salinity and Conservative Temperature. Here, practical salinity and potential temperature are used to match what is used by the equation of state in the ocean model whose data shall be analysed. They are related by $\nu=\rho^{-1}=\mathcal{V}(S,\theta,p)$ where $\mathcal{V}$ is the (inverse of the) equation of state. Also let $\nu_{p}=\partial_{p}\mathcal{V}(S,\theta,p)$ be another 3D scalar field, akin to the compressibility. Let $\bm{\hat{i}}$, $\bm{\hat{j}}$, and $\bm{\hat{k}}$ be the eastward, northward, and upward unit vectors, respectively. The depth $z$ is another 3D scalar field, giving the distance in the $\bm{\hat{k}}$ direction from each point to a reference geopotential (near the sea surface), and signed so that $z<0$ below the reference geopotential. The sea surface height $\eta$ is a 2D scalar field, similar to $z$ but for the sea surface. Let $g$ be the gravitational acceleration.

2.2 Background of geostrophic streamfunctions

The geostrophic velocity $\bm{u}_{g}$ is defined by the geostrophic approximation of the horizontal momentum equation,

[TABLE]

where $f$ is the Coriolis parameter, and $\nabla_{z}p\equiv\partial_{x}p|_{z}\,\bm{\hat{i}}+\partial_{y}p|_{z}\,\bm{\hat{j}}$ is the gradient of pressure at constant depth $z$. A GSF $\Psi$ in a generic $r$-surface (e.g. an isosurface of a 3D scalar field $r$) is defined as that which satisfies333For some surfaces, it is useful to instead require $\nabla_{r}\Psi=\nabla_{z}p$, which gives a streamfunction for $\rho f\bm{u}_{g}$ rather than $f\bm{u}_{g}$ [McDougall, 1989]. Also, the geopotential of the sea surface is sometimes included in the definition, i.e. requiring $\nabla_{r}\Psi=\nu\nabla_{z}p-\nabla(g\eta)$ [McDougall and Klocker, 2010]. As $g\eta$ can be absorbed into $\Psi$, the two definitions are essentially interchangeable.

[TABLE]

where $\nabla_{r}$ is the “projected non-orthogonal gradient” in the $r$-surface, first introduced by Starr [1945]. Specifically, given a tracer $C$ (a 3D scalar field), $\nabla_{r}C\equiv\partial_{x}C|_{r}\,\bm{\hat{i}}+\partial_{y}C|_{r}\,\bm{\hat{j}}$, where the partial derivatives are taken “in the $r$-surface”, sampling $C$ from the $r$-surface even though the radial (vertical) position is ignored when measuring distance [McDougall et al., 2014].

The first step to determine $\Psi$ is usually to transform the gradient of $p$ from one at constant depth to one in the surface, according to $\nabla_{r}p=\nabla_{z}p+\partial_{z}p\,\nabla_{r}z$ [Starr, 1945]. This is combined with hydrostatic balance, $\partial_{z}p=-g\nu^{-1}$, transforming (2) into

[TABLE]

treating $g$ as constant, for simplicity of presentation.444More generally, the gravitational potential $\Phi$ can be used as the vertical coordinate instead of $z$, in which case hydrostatic balance is $\mathrm{d}p/\mathrm{d}\Phi=-\nu^{-1}$, and (3) becomes $\nabla_{r}\Psi=\nu\nabla_{r}p+\nabla_{r}\Phi$.

We use an alternative notation that captures the precise details of the projected non-orthogonal gradient, but uses a standard gradient on a transformed variable [S19]. Specifically, let $\utilde{C}$ be the projection onto a perfect sphere (the centre and radius of which match those of Earth) of the restriction of a 3D field $C$ to the surface in question.555For simplicity, assume the surface in question exists at a unique depth in each water column, or not at all where it has grounded or outcropped. Were it not, simply restrict attention to a local region where this is true, and work region by region. With $C$ a scalar field, $\nabla\utilde{C}$ can be evaluated using the standard gradient in spherical coordinates. The only difference between $\nabla\utilde{C}$ and $\nabla_{r}C$ is that the former deals with a single 2D surface under consideration, whereas the latter deals with the whole 3D ocean. Note $C$ could instead be a vector field, such as $\bm{u_{g}}$.

Considering a single $r$-surface with this notation, (3) becomes

[TABLE]

The goal is to move $\utilde{\nu}$ inside the gradient, expressing the RHS as the gradient of some quantity, namely $\utilde{\Psi}$. The vector field $f\utilde{\bm{u}_{g}}$ can then be determined entirely from the scalar field $\utilde{\Psi}$.

2.3 Background of neutral surfaces

As discussed in Section 1, neutral surfaces are formally ill-defined. Only under special conditions are neutral surfaces well-defined. A necessary condition is that the neutral helicity is everywhere zero [McDougall and Jackett, 1988]. To make theoretical progress, let us assume neutral surfaces are well-defined. Later, we will translate the theoretical results into the realistic case where this is not so.

A neutral surface is one in which specific volume variations in the surface are solely due to pressure variations in the surface. Mathematically, this is usually expressed as $\nabla_{n}\nu=\nu_{p}\nabla_{n}p$, using the projected non-orthogonal gradient in a neutral surface denoted $n$ [McDougall, 1987]. With the under-tilde notation, a neutral surface is one that satisfies

[TABLE]

An equivalent and perhaps more familiar condition is

[TABLE]

where $\nu_{S}=\partial_{S}\mathcal{V}(S,\theta,p)$ and $\nu_{\theta}=\partial_{\theta}\mathcal{V}(S,\theta,p)$. This is all we need to proceed, but for a fuller review of neutral surfaces, see S19.

2.4 The exact geostrophic streamfunction on a neutral surface, locally

The neutral condition (5) implies the gradients, and thus the contours, of $\utilde{\nu}$ and $\utilde{p}$ are always aligned. Therefore, there is a functional relationship between $\utilde{\nu}$ and $\utilde{p}$ :

[TABLE]

for some function $\hat{\nu}$. In fact, $\hat{\nu}$ is a multivalued function: $\utilde{\nu}$ is constant on a closed contour of $\utilde{p}$, but there can be different values of $\utilde{\nu}$ on different contours of $\utilde{p}$ at the same pressure value.

For now, consider a region where $\hat{\nu}$ is single-valued. In this region, the exact GSF for a neutral surface is

[TABLE]

for some constant pressure $p_{0}$ in the domain of $\hat{\nu}$. Note $p_{0}$ provides an arbitrary additive constant to $\utilde{\Psi}^{(n)}$. That $\utilde{\Psi}^{(n)}$ exactly satisfies (4) is easily confirmed by taking its gradient using the Leibniz integral rule, then using (7).

2.5 Discussion

In (8), $\utilde{z}$ is obtained by inverting hydrostatic balance to get $\partial_{p}z=-g^{-1}\nu$, then integrating in the local water column to get

[TABLE]

where $p^{(\eta)}$ is the pressure at $z=\eta$, the sea surface. (If $p^{(\eta)}$ is constant, it is common, but not necessary, to set $p_{0}=p^{(\eta)}$.) This transformation is standard when studying GSFs. For example, multiplying (9) by $g$ and taking the gradient in an isobaric surface (on which $\utilde{p}$ is constant) yields $\nabla(g\utilde{z})$ on the LHS, which is the entire RHS of (4) in this case. Thus

[TABLE]

is the exact GSF in an isobaric surface; this is the well-known “dynamic height”.666Numerical calculations of the dynamic height can be made more accurate by using the specific volume anomaly in the integrand of (10); this is then called the “dynamic height anomaly”. With double precision computing, this is no longer necessary [McDougall and Klocker, 2010]. Though $\utilde{\Psi}^{(p)}$ is really just $g\utilde{z}$, the form (10) is essential because it is computable from hydrographic casts which measure the salinity and potential temperature as functions of pressure: denoting these $\check{S}(p)$ and $\check{\theta}(p)$ respectively at a given water column, the integrand in (9) and (10) is a function of $p$ alone, namely $\mathcal{V}(\check{S}(p),\check{\theta}(p),p)$. With model data, the actual depth of an isobaric or other such surface may be known, but this is not so for real ocean measurements.

One might argue that (8) is not a closed form expression for $\utilde{\Psi}^{(n)}$, since $\hat{\nu}$ is unspecified. Or is it? We are accustomed to think there is ambiguity about neutral surfaces, given their ill-defined nature in the real ocean with non-zero neutral helicity, and this ambiguity transfers to $\hat{\nu}$. But if we actually — hypothetically — have a well-defined neutral surface, then we do have the function $\hat{\nu}$. It is determined by the data $(\utilde{p},\utilde{\nu})$: a scatter plot of these would show a functional relationship around which there is precisely zero scatter. Also, consider that $\utilde{\Psi}^{(p)}$ has an analogous definition: it too is defined using a function of specific volume, and this function is determined from data. In this case, the data is taken vertically above or below the location in question, which is easy enough to fathom. For $\utilde{\Psi}^{(n)}$, the problem lies in determining the local neutral surface, from which to get the data $(\utilde{p},\utilde{\nu})$. Assuming the neutral surface is known, $\utilde{\Psi}^{(n)}$ is the exact GSF on that surface.

In fact, $\utilde{\Psi}^{(n)}$ from (8) is the exact GSF on any surface for which $\utilde{\nu}$ is a function of $\utilde{p}$. One example is an isosurface of $\nu$, but this is very nearly a geopotential, far from neutral, and so not terribly useful. Dividing (5) by $\utilde{\nu_{p}}$, the entire argument can be reversed, and any surface where $\utilde{p}$ is a function of $\utilde{\nu}$ also possesses an exact GSF — such as an isobaric surface, again far from neutral. More generally, isosurfaces of any function $\hat{\sigma}$ that is pycnotropic (a function only of $p$ and $\nu$) possess an exact GSF [de Szoeke, 2000]. This is because an isosurface of $\sigma=\hat{\sigma}(p,\nu)$ has $0=\nabla\utilde{\sigma}=\utilde{\partial_{\nu}\sigma}\nabla\utilde{\nu}+\utilde{\partial_{p}\sigma}\nabla\utilde{p}$; thus $\nabla\utilde{\nu}$ and $\nabla\utilde{p}$ are parallel, as in (5) but replacing $\utilde{\nu_{p}}$ by $-\utilde{\partial_{p}\sigma}\,/\,\utilde{\partial_{\nu}\sigma}$, and the derivation leading to (8) proceeds similarly. By empirically determining a pycnotropic $\hat{\sigma}$ from oceanic data, de Szoeke et al. [2000] created orthobaric density, isosurfaces of which are approximately neutral. Orthobaric density surfaces possess an exact GSF, and its form is identical to $\utilde{\Psi}^{(n)}$ in (8) — see Eq. 2.32 of de Szoeke et al. [2000]. What is new about (8) is the realization that this form applies to neutral surfaces — or rather to small regions of neutral surfaces where $\hat{\nu}$ is single-valued. For neutral surfaces more broadly, $\hat{\nu}$ is multivalued, and we must generalize (8) to handle this.

2.6 The global structure of $\utilde{\Psi}^{(n)}$

The view so far has been local, with $\hat{\nu}$ a single-valued function. From a global perspective, $\hat{\nu}$ is multivalued: $\utilde{\nu}$ is constant on a contour of $\utilde{p}$, but may differ between disjoint contours of $\utilde{p}$ at the same pressure value. How this multivalued nature enters (8) must be made explicit.

Consider an arbitrary path $\mathcal{P}$ in a neutral surface from $\bm{x}_{0}$ to $\bm{x}$. (This is a neutral trajectory.) For simplicity of presentation, suppose the neutral surface is path-connected, so that all points $\bm{x}$ in the neutral surface are reachable from $\bm{x}_{0}$ by such a path. (If this is false, simply apply the following theory to each path-connected component separately.) Using the gradient theorem, the global structure of $\utilde{\Psi}^{(n)}$ is

[TABLE]

Using (4) and (7) and choosing $\utilde{\Psi}^{(n)}(\bm{x}_{0})=g\utilde{z}(\bm{x}_{0})$ for convenience, (11) becomes

[TABLE]

It is by the neutral trajectory $\mathcal{P}$ that the multivalued nature of $\hat{\nu}$ enters $\utilde{\Psi}^{(n)}$ in (12), whereas this was less clear from (8). Of course, to use the gradient theorem we assumed $\utilde{\Psi}^{(n)}$ is well-defined. But is $\utilde{\Psi}^{(n)}(\bm{x})$ independent of the choice of the path $\mathcal{P}$ from $\bm{x}_{0}$ to $\bm{x}$? To answer this, the Reeb graph is introduced.

A contour of $\utilde{p}$ is a connected component of the level set of $\utilde{p}$ at a specified pressure value. A level set can be the disjoint union of multiple contours. The Reeb graph of $\utilde{p}$ contracts each contour of $\utilde{p}$ to a single point. The essence of the resulting object can be represented as a graph, a collection of $N$ nodes and $A$ arcs between pairs of nodes. Each node $n$ corresponds to a critical point of $\utilde{p}$, denoted $\bm{x}_{n}$, at which point the critical value is denoted $p_{n}$; that is, $p_{n}=\utilde{p}(\bm{x}_{n})$. Leaf nodes correspond to extrema of $\utilde{p}$, while internal nodes correspond to saddle points of $\utilde{p}$. Each arc $a$ corresponds to a geographic region $\mathcal{R}_{a}$ within which there is precisely one $\utilde{p}$ contour per pressure value. In this way, the Reeb graph partitions space into geographic regions within each of which the multivalued function $\hat{\nu}$ is actually single-valued. Restricting the multivalued $\hat{\nu}$ to $\mathcal{R}_{a}$ gives a single-valued function, denoted $\hat{\nu}_{a}$ and called a branch of $\hat{\nu}$. An example is shown in the left and centre panels of Fig. 1. For oceanographic data, many of the associated regions are the shape of mesoscale eddies, which form closed $\utilde{p}$ contours.

Now re-consider the neutral trajectory $\mathcal{P}$. Each point on $\mathcal{P}$ is part of a contour of $\utilde{p}$ that is contracted to some point on the Reeb graph. In this way, $\mathcal{P}$ corresponds to a walk777In graph theory, a walk is an alternating sequence of nodes and arcs where each node incident upon each arc adjacent to it in the sequence. Moreover, the sequence must start and end with nodes, but this is relaxed here. through the Reeb graph, alternately moving along arcs and passing through nodes, in a sequence $a_{0}$, $n_{1}$, $a_{1}$, …, $n_{\mathcal{J}}$, $a_{\mathcal{J}}$. Moreover, $\bm{x}_{0}$ is in $\mathcal{R}_{a_{0}}$, so $\utilde{\nu}(\bm{x}_{0})=\hat{\nu}_{a_{0}}(p_{0})$, where $p_{0}=\utilde{p}(\bm{x}_{0})$. Similarly, $\bm{x}$ is in $\mathcal{R}_{a_{\mathcal{J}}}$, so \utilde{\nu}(\bm{x})=\hat{\nu}_{a_{\mathcal{J}}}\big{(}\utilde{p}(\bm{x})\big{)}. Thus, the path integral in (12) becomes a “graph integral”,

[TABLE]

This reduces to (8) when $\mathcal{J}=0$, as happens when $\bm{x}$ and $\bm{x}_{0}$ are in the same region $\mathcal{R}_{a}$, and $\mathcal{P}$ stays within this region. Note that (2.6) is analogous to Eq. 20 of S19, in which $\hat{\nu}$ is obtained by integrating a multivalued function $\hat{\pi}$ that satisfies $\utilde{\nu_{p}}=\hat{\pi}(\utilde{p})$.888S19 presents the theory of topobaric surfaces using in-situ density rather than specific volume, but it is trivial to use the latter instead.

What are the circumstances that guarantee $\utilde{\Psi}^{(n)}$ is well-defined? If all paths from $\bm{x}_{0}$ to $\bm{x}$ are equivalent to the same walk through the Reeb graph, then (2.6) guarantees path-independence of $\utilde{\Psi}^{(n)}$. Of course, a path could backtrack on itself, so that its walk has repeated nodes and differs from the walk of another path between the same end points. Fortunately, such backtracking is perfectly cancelled by the integrals in (2.6). Thus, $\utilde{\Psi}^{(n)}$ is well-defined if, for every pair of nodes in the Reeb graph with a walk between them, there is a unique walk between them having no repeated nodes.

In other words, $\utilde{\Psi}^{(n)}$ is well-defined if the Reeb graph of $\utilde{p}$ is a tree, i.e. it contains no cycles. (If the neutral surface is not path-connected, $\utilde{\Psi}^{(n)}$ is well-defined if the Reeb graph is a forest, i.e. a collection of trees.) A cycle is a walk in the graph that starts and ends at the same node, and contains no repeated nodes or arcs, aside from the first and last node. Cycles can arise in the Reeb graph of $\utilde{p}$ when the surface has holes, such as made by islands, seamounts, and other places where the surface grounds or outcrops. However, not every hole in the surface produces a cycle in the Reeb graph. For example, if a $\utilde{p}$ contour encloses a single hole, that hole does not produce a cycle [S19]. The small island in Fig. 1 does not produce a cycle, whereas the big island does. Thus, one cannot assess whether $\utilde{\Psi}^{(n)}$ is well-defined or not based purely on the existence of holes in the neutral surface: the Reeb graph is needed.

There is a second, special circumstance in which $\utilde{\Psi}^{(n)}$ is well-defined, even when the Reeb graph possesses cycles. Let $n_{1}$, $a_{1}$, $n_{2}$, $a_{2}$, $...$, $a_{\mathcal{J}-1}$, $n_{\mathcal{J}}=n_{1}$ be a cycle (all $a_{j}$ are distinct), and let $\bm{x}_{0}=\bm{x}=\bm{x}_{n_{1}}$ without loss of generality. If $\utilde{\Psi}^{(n)}$ is well-defined, then (2.6) along this cycle reduces to

[TABLE]

Noting this cycle was arbitrary, (14) must hold for every such cycle. If it does, then all paths between a given pair of endpoints yield the same value for $\utilde{\Psi}^{(n)}$, despite the paths potentially taking distinct walks. Thus, $\utilde{\Psi}^{(n)}$ is well-defined if and only if (14) holds for every cycle in the Reeb graph of $\utilde{p}$ (which is trivially true when the graph is a tree).999In fact, this need only hold on a subset of cycles, namely a cycle basis; see S19 for details.

In general, the Reeb graph has cycles and there is no reason for (14) to hold, so $\utilde{\Psi}^{(n)}$ is ill-defined. Indeed, $\utilde{\Psi}^{(n)}$ is a multivalued function of geographic location, because it can be evaluated using different path integrals that loop, any number of times, around holes in the neutral surface. Pictured as a surface with radial (vertical) coordinate given by its value, $\utilde{\Psi}^{(n)}$ qualitatively resembles a multistorey car park, with interior ramps around holes in the neutral surface (Fig. 2 of S19 illustrates this, while simultaneously illustrating a similar phenomenon about neutral surfaces themselves). The “pitch” of these interior ramps is the (generally non-zero) RHS of (14). A map of $\utilde{\Psi}^{(n)}$ — that chooses a single value of $\utilde{\Psi}^{(n)}$ for each geographic location — will exhibit discontinuities that emanate from islands and other holes in the neutral surface, as illustrated in Fig. 1. The proof by McDougall [1989] for the existence of $\utilde{\Psi}^{(n)}$ was local in scope, so this issue could not have been foreseen.

The multivalued nature of $\utilde{\Psi}^{(n)}$ does not, of course, mean the geostrophic velocity is multivalued or discontinuous. The geostrophic velocity is determined not by $\utilde{\Psi}^{(n)}$, but by its gradient; the path-ambiguity of $\utilde{\Psi}^{(n)}$ adds a constant offset to each “storey” of $\utilde{\Psi}^{(n)}$ and so does not change its gradient (evaluated using a consistent “storey”). Indeed, $\nabla\utilde{\Psi}^{(n)}$ is determined by the gradient of (12) or (2.6), which just reverts to the unique value given by (4) upon using (7). So, if one desires only the geostrophic velocity, simply use (4). But often, more advanced analyses (as discussed in Section 1) require the GSF itself — and require it to be well-defined, without discontinuities emanating from holes. Remedying these discontinuities is the goal for the next section.

3 The topobaric geostrophic streamfunction

From the previous section’s theoretical results, a practical idea emerges for a well-defined GSF that approximates $\utilde{\Psi}^{(n)}$ and is useful on any approximately neutral surface. This approximation is called the topobaric GSF and denoted $\utilde{\Psi}^{(tb)}$ because it is built upon a topological analysis of the pressure on an approximately neutral surface.

Given an approximately neutral surface on which the pressure is $\utilde{p}$, $\utilde{\Psi}^{(tb)}$ is constructed as follows. First, calculate the Reeb graph of $\utilde{p}$. Then, obtain a multivalued function $\hat{\nu}$ that satisfies (14) for each cycle in the Reeb graph, and whose branches $\hat{\nu}_{a}$ approximately satisfy (7) in each region, i.e. \utilde{\nu}(\bm{x})\approx\hat{\nu}_{a}\big{(}\utilde{p}(\bm{x})\big{)} for all $\bm{x}$ in $\mathcal{R}_{a}$. Then, obtain $\utilde{\Psi}^{(tb)}$ by integrating $\hat{\nu}$ according to the RHS of (2.6). Because $\hat{\nu}$ satisfies the cycle constraints (14), $\utilde{\Psi}^{(tb)}$ is well-defined.

When the surface in question is a topobaric surface, $\hat{\nu}$ is given a priori. Otherwise, $\hat{\nu}$ must be empirically determined. These two cases are discussed next.

3.1 Use on topobaric surfaces

It should be emphasized now that $\utilde{\Psi}^{(tb)}$ is not the exact GSF on a topobaric surface. A topobaric surface [S19] is a well-defined surface that satisfies (7) not just approximately, but exactly, and $\hat{\nu}$ is obtained by integrating another multivalued function $\hat{\pi}$ which approximately satisfies $\utilde{\nu_{p}}\approx\hat{\pi}(\utilde{p})$. (To ensure topobaric surfaces are well-defined, $\hat{\pi}$ satisfies cycle constraints identical to (14) but the integrands use $\hat{\pi}$ rather than $\hat{\nu}$ — see Eq. 21 of S19.) The resulting $\hat{\nu}$ does not, in general, satisfy the cycle constraints (14). Like true neutral surfaces, topobaric surfaces do possess an exact GSF, but it is the multivalued $\utilde{\Psi}^{(n)}$, not the well-defined $\utilde{\Psi}^{(tb)}$.

However, $\utilde{\Psi}^{(tb)}$ is the exact GSF on a modified topobaric surface, which is a topobaric surface in which $\hat{\nu}$ satisfies the cycle constraints (14). The methods of S19 to create topobaric surfaces are easily extended to create modified topobaric surfaces, simply by adding the cycle constraints (14) when empirically fitting $\hat{\pi}$. With $\hat{\nu}$ expressed in terms of $\hat{\pi}$ (according to Eq. 20 of S19, analogous to (2.6) here), some algebra reduces (14) to

[TABLE]

which must be satisfied for each cycle $n_{1}$, $a_{1}$, $...$, $a_{\mathcal{J}-1}$, $n_{\mathcal{J}}=n_{1}$ in the Reeb graph of $\utilde{p}$. These additional constraints tend to make modified topobaric surfaces slightly less neutral than regular topobaric surfaces, but their exact GSF, $\utilde{\Psi}^{(tb)}$, is well-defined.

A goal for Section 7 is to demonstrate, numerically, that $\utilde{\Psi}^{(tb)}$ is exact on modified topobaric surfaces.

3.2 Use on other approximately neutral surfaces

Now consider approximately neutral surfaces more generally. No exact GSF exists for these surfaces in general, but $\utilde{\Psi}^{(tb)}$ can be constructed as an approximate GSF. Now, $\hat{\nu}$ is not given a priori, and (7) is only approximate. Nonetheless, $\hat{\nu}$ can be empirically fit, then integrated to obtain $\utilde{\Psi}^{(tb)}$. This is exactly analogous to how, for topobaric surfaces, $\hat{\pi}$ is empirically fit, then integrated to obtain $\hat{\nu}$. Hence, the method is only summarized here; see S19 for further details.

First, a reference location $\bm{x}_{0}$ is chosen: following S19, $\bm{x}_{0}=(180^{\circ}\mathrm{E},0^{\circ}\mathrm{N})$ by default. This is used to select a reference salinity $S_{0}=\utilde{S}(\bm{x}_{0})$ and a reference potential temperature $\theta_{0}=\utilde{\theta}(\bm{x}_{0})$.

Second, the Reeb graph of $\utilde{p}$ is calculated by the algorithm of Doraiswamy and Natarajan [2013].

Third, the branches of $\hat{\nu}$ are empirically fit in each region. Given an arc $a$, denote the two nodes incident to $a$ as $\ell_{a}$ and $h_{a}$, with the critical values of $\utilde{p}$ at these nodes satisfying $p_{\ell_{a}}<p_{h_{a}}$. Using a simple functional form

[TABLE]

the unknown constants $K_{a}$ and $L_{a}$ are determined by fitting $K_{a}+L_{a}(\utilde{p}-p_{\ell_{a}})$ to the specific volume anomaly on the surface, $\utilde{\delta}=\utilde{\nu}-\mathcal{V}(S_{0},\theta_{0},\utilde{p})$, using ordinary least squares and restricting data to the region $\mathcal{R}_{a}$. Including $\mathcal{V}(S_{0},\theta_{0},p^{\prime})$ in (16) helps to capture some of the non-linearity in the equation of state. Each branch is fit independently, except for those whose arcs exist on cycles of the Reeb graph: these branches are also fit by ordinary least squares, but as a coupled problem subject to the cycle constraints (14). For truly neutral surfaces and (modified or regular) topobaric surfaces, $\hat{\nu}$ meets continuously at the pressure saddles [S19]; however, we are using $\utilde{\Psi}^{(tb)}$ on other surfaces, so $\hat{\nu}$ is allowed to meet discontinuously at the pressure saddles. This is analogous to how topobaric surfaces allow the branches of $\hat{\pi}$ to meet discontinuously at the pressure saddles.

Reassuringly, any such discontinuities are eliminated in the fourth step, which obtains $\utilde{\Psi}^{(tb)}$ by integrating $\hat{\nu}$ according to the RHS of (2.6). Practically, this is done by a breadth-first search in the Reeb graph: at each step a node $m$ is discovered that is adjacent to a previously discovered node $n$ (to initialize the search, one chosen node is marked as discovered), and the arc $a$ that is incident to both $m$ and $n$ is determined; then $\utilde{\Psi}^{(tb)}$ in the region $\mathcal{R}_{a}$ is obtained from (8) using the local branch $\hat{\nu}_{a}$ and using $p_{0}=p_{n}$, so that $\utilde{\Psi}^{(tb)}$ matches continuously at $\bm{x}_{n}$. At the end of the breadth-first search, every arc has been processed except one arc $a$ per cycle. Obtain $\utilde{\Psi}^{(tb)}$ in each $\mathcal{R}_{a}$ as above, using $p_{0}=p_{l_{a}}$. The cycle constraints (14) ensure the result is identical had we instead used $p_{0}=p_{h_{a}}$.

4 The Montgomery potential and variants

We now shift gears, reviewing the Montgomery [1937] potential and its variants, as preparation for deriving another new GSF in Section 5. The Montgomery potential is the exact GSF on specific volume anomaly surfaces, i.e. isosurfaces of the specific volume anomaly

[TABLE]

where $\nu_{0}=\mathcal{V}(S_{0},\theta_{0},p)$ is the specific volume of a fluid parcel with reference salinity $S_{0}$ and reference potential temperature $\theta_{0}$, at the local pressure.

To derive the Montgomery potential, add and subtract $\utilde{\nu_{0}}\nabla\utilde{p}$ from the RHS of (4). The subtracted term combines with $\utilde{\nu}\nabla\utilde{p}$ to give $\utilde{\delta}\nabla\utilde{p}$, upon which the product rule is used; the added term has the Leibniz integral rule applied to it. The result is

[TABLE]

for some constant pressure $P$. On an isosurface of $\delta$, the last term is zero, so $\utilde{\Psi}^{(\delta)}$ is just the expression in parentheses. As before, the arbitrary $P$ provides an arbitrary additive constant to $\utilde{\Psi}^{(\delta)}$, and $\utilde{z}$ is given by (9). If $p^{(\eta)}$ in (9) is constant, we may choose $P=p^{(\eta)}$, so that

[TABLE]

Substituting (19) into (18) gives a more familiar expression for $\utilde{\Psi}^{(\delta)}$. We shall continue with the more general form for non-constant $p^{(\eta)}$, but if $p^{(\eta)}$ is constant, one can also substitute (19) in any of (4), (4), (22), and (5) below.

On surfaces other than specific volume anomaly surfaces, $\utilde{\delta}$ is non-constant, so $\utilde{p}\nabla\utilde{\delta}$ in (18) is non-zero; this creates an error in estimating the geostrophic velocity by use of $\utilde{\Psi}^{(\delta)}$. To minimize this error, Zhang and Hogg [1992] subtracted a constant pressure $p_{0}$ from $\utilde{p}$ in (4), and percolated $p_{0}$ through the above derivation. Equivalently, add and subtract $\nabla(\utilde{\delta}\,p_{0})$ from (18), to obtain

[TABLE]

having used $\utilde{\delta}\nabla p_{0}=\bm{0}$. When the surface in question is not an isosurface of $\delta$, the error in using $\utilde{\Psi}^{(ZH)}$ can be made less than that of $\utilde{\Psi}^{(\delta)}$ by choosing $p_{0}$ to reduce the prefactor $(\utilde{p}-p_{0})$. Zhang and Hogg [1992] took $p_{0}$ as the mean of $\utilde{p}$, which minimizes the root-mean-square of $(\utilde{p}-p_{0})$. Traditionally, $(S_{0},\theta_{0})$ is taken as ($35\text{\,}\mathrm{p}\mathrm{s}\mathrm{u}$, $0^{\circ}\mathrm{C}$), following Montgomery [1937].

McDougall and Klocker [2010] augmented $\utilde{\Psi}^{(ZH)}$ (which is designed for $\delta$-surfaces) with additional information appropriate to neutral surfaces, obtaining

[TABLE]

which is their Eq. 62 re-expressed using (19) and with $g\eta$ added. To derive (4), $T_{b}/\utilde{\rho}$ is treated as constant, namely $2.7\text{\times}{10}^{-15}\text{\,}\mathrm{K}^{-1}\,\mathrm{P}\mathrm{a}^{-2}\,\mathrm{m}^{2}\,\mathrm{s}^{-2}$.

The choice of reference values $p_{0}$, $S_{0}$, and $\theta_{0}$ does affect the errors. For $\utilde{\Psi}^{(ZH)}$, one would ideally choose $p_{0}$ in tandem with $S_{0}$ and $\theta_{0}$ to minimize $|(\utilde{p}-p_{0})\nabla\utilde{\delta}|$. Such optimization efforts will not be pursued here, as there are greater gains to be had by using other GSFs. Instead, we follow the suggestion by McDougall and Klocker [2010] to take $p_{0}=\utilde{p}(\bm{x}_{0})$, $S_{0}=\utilde{S}(\bm{x}_{0})$, and $\theta_{0}=\utilde{\theta}(\bm{x}_{0})$ with $\bm{x}_{0}$ in the equatorial Pacific; for consistency with Section 3, we take $\bm{x}_{0}=(180^{\circ}\mathrm{E},0^{\circ}\mathrm{N})$. (This will indeed give $p_{0}$ quite close to the mean of $\utilde{p}$ over the whole surface.)

Finally, the Cunningham [2000] GSF is [following Eq. 25 of McDougall and Klocker, 2010]

[TABLE]

Note that $\utilde{\Psi}^{(Cu)}$ is similar to $\utilde{\Psi}^{(\delta)}$ but uses water-column specific “reference” values, namely $\utilde{S}$ and $\utilde{\theta}$. This re-defines $\delta$ as $\nu-\mathcal{V}(\utilde{S},\utilde{\theta},p)$, leading to $\utilde{\delta}=\utilde{\nu}-\mathcal{V}(\utilde{S},\utilde{\theta},\utilde{p})=0$. However, a new type of error is created due to gradients of the “reference” values: by the Leibniz integral rule,

[TABLE]

The last term is the error, which can be minimized by careful choice of $P$. For consistency, we again use $\bm{x}_{0}=(180^{\circ}\mathrm{E},0^{\circ}\mathrm{N})$ and set $P=\utilde{p}(\bm{x}_{0})$.

5 The orthobaric Montgomery potential

We now generalize the progress that Zhang and Hogg [1992] made on the Montgomery [1937] GSF. Let $\hat{p}$ be a function of $\delta$. Add and subtract $\hat{p}(\utilde{\delta})\nabla\utilde{\delta}$ from (18), using Leibniz’s integral rule on the subtracted term, to obtain

[TABLE]

for some constant $\Delta$. The expression in large parentheses is a generalization of the Montgomery potential.

In a surprising twist, $\utilde{\Psi}^{(oM)}$ has merely re-expressed $\utilde{\Psi}^{(n)}$ in the form of a Montgomery potential, as follows:

[TABLE]

In the second equality, $\hat{\delta}$ is defined by $\hat{\nu}$ and the reference profile $\nu_{0}$: specifically, $\hat{\delta}(p^{\prime})=\hat{\nu}(p^{\prime})-\mathcal{V}(S_{0},\theta_{0},p^{\prime})$. The third equality uses Laisant’s inverse integral rule, with $\hat{p}$ the inverse function of $\hat{\delta}$. (If the latter is not invertible, first replace the integral of $\hat{\delta}$ with a sum of integrals over restricted domains on which $\hat{\delta}$ is invertible.) Finally, $P\Delta$ is an arbitrary constant, which can be ignored.

Perhaps this equivalence between $\utilde{\Psi}^{(oM)}$ and $\utilde{\Psi}^{(n)}$ should not be a surprise. The error in using $\utilde{\Psi}^{(oM)}$ to estimate the geostrophic velocity can be made small by choosing $\hat{p}$ so that $\hat{p}(\utilde{\delta})$ stays close to $\utilde{p}$. This is eminently possible if the surface in question is a neutral surface: then $\utilde{p}$ and $\utilde{\nu}$ are perfectly functionally related, hence so too are $\utilde{p}$ and $\utilde{\delta}=\utilde{\nu}-\mathcal{V}(S_{0},\theta_{0},\utilde{p})$. This gives $\hat{p}$ such that $\utilde{p}-\hat{p}(\utilde{\delta})=0$ identically, making $\utilde{\Psi}^{(oM)}$ exact on neutral surfaces. Of course, this requires $\hat{p}$ to be a multivalued function of $\delta$, with branches defined on the Reeb graph of $\utilde{\delta}$. In this case, $\utilde{\Psi}^{(oM)}$ re-expresses $\utilde{\Psi}^{(n)}$ from (2.6) in the form of the Montgomery potential, and we would call it the topobaric Montgomery potential. However, this is little different from the topobaric GSF since the Reeb graph ensures all branches are single-valued.

Instead, the goal here is simplicity, so we take $\hat{p}$ as a single-valued function of $\delta$, and call $\utilde{\Psi}^{(oM)}$ the orthobaric Montgomery potential. This re-expresses the orthobaric geostrophic streamfunction — $\utilde{\Psi}^{(n)}$ from (8) with a single-valued $\hat{\nu}$ — in the form of a Montgomery potential. The function $\hat{p}$ is empirically fit over an entire (connected part of the) surface; as such, $\utilde{\Psi}^{(oM)}$ is only an approximate GSF on surfaces other than specific volume anomaly surfaces.

Numerical tests have persuaded the author that the orthobaric Montgomery potential is typically more accurate than the orthobaric GSF, despite their mathematical equivalence. This is because $\utilde{\nu}$ (and $\utilde{\delta}$) can differ substantially between the Arctic and the Southern Ocean, say, even at the same pressure: it is important to handle $\hat{\nu}$ and $\hat{\delta}$ as multivalued. On the other hand, $\hat{p}$ can more closely approximate $\utilde{p}$ as a single-valued function of $\utilde{\delta}$, provided the reference values $S_{0}$ and $\theta_{0}$ are well-chosen. (Think of $\utilde{p}$ as a quadratic function of $\utilde{\delta}$, nicely single-valued, whereas $\utilde{\delta}$ is like a square-root of $\utilde{p}$, multivalued. Data illustrating this point will be shown in Section 7.2.) To separate $\utilde{\delta}$ in the Arctic from $\utilde{\delta}$ in the Southern Ocean, we will choose $S_{0}=\utilde{S}(\bm{x}_{\mathrm{deep}})$ and $\theta_{0}=\utilde{\theta}(\bm{x}_{\mathrm{deep}})$, where $\bm{x}_{\mathrm{deep}}$ is such that $\utilde{p}(\bm{x}_{\mathrm{deep}})$ is a maximum. Other values for $(S_{0},\theta_{0})$ may be chosen, but this simple method seems to yield good results on a variety of surfaces. Then, $\hat{p}$ will be obtained by fitting $\utilde{p}$ to $\utilde{\delta}$ as a cubic spline with 12 pieces. Again, this form seems to robustly yield good results, but other choices are possible.

6 Gridded data and Boussinesq ocean models

Numeric calculations in the next section will use ECCO2 data [Menemenlis et al., 2005] on 22–24 December 2002. The seawater Boussinesq approximation [Young, 2010], employed by ECCO2, changes the above theory in two important, yet easy-to-handle ways. First, Boussinesq ocean models swap the in-situ density $\rho$ for a constant reference density $\rho_{B}$ in the horizontal momentum equations, so (1)–(4) become

[TABLE]

Now, the task of a GSF is to bring $\utilde{\rho}$ inside the gradient with $\utilde{z}$. It is easier to work with $\rho$ rather than $\nu$. Second, the in-situ density (still used in the vertical momentum equation, perhaps reduced to hydrostatic balance) is calculated from a Boussinesq equation of state $B$ that uses depth $z$ rather than the in-situ pressure [Young, 2010]. Specifically,

[TABLE]

Thus the neutral surface relation (5) becomes

[TABLE]

where $\rho_{z}=\partial_{z}B(S,\theta,z)$. In the seawater Boussinesq approximation, the multivalued function of pressure, $\hat{\nu}$, gets replaced by a multivalued function of depth, $\hat{\rho}$. The essential ideas, though, are unchanged. See the appendices for Boussinesq formulas and their numerical discretisation.

7 Numeric comparison of geostrophic velocity errors

The exact GSF for neutral surfaces, $\utilde{\Psi}^{(n)}$, is inexact on approximately neutral surfaces. Also, $\utilde{\Psi}^{(tb)}$ differs slightly from $\utilde{\Psi}^{(n)}$ because of the requirement to be well-defined, and because the form (16) of the empirically fit functions is limiting. Nonetheless, $\utilde{\Psi}^{(tb)}$ may estimate the geostrophic velocity more accurately than other GSFs. The closer to neutral the surface, the more likely $\utilde{\Psi}^{(tb)}$ is to outperform other GSFs. Numerical tests are required to assess the value of $\utilde{\Psi}^{(tb)}$, as well as $\utilde{\Psi}^{(oM)}$.

7.1 Setup of numerical tests

Three types of approximately neutral surfaces are computed:

$\sigma_{1}$- or $\sigma_{2}$-surfaces, i.e. isosurfaces of potential density [Wüst, 1935]; 2. 2.

$\tau^{\prime}$-surfaces, i.e. modified topobaric surfaces [ S19, and this manuscript]; and 3. 3.

$\omega$-surfaces [Klocker et al., 2009].

Two of each of these surfaces is computed, intersecting ($180^{\circ}\mathrm{E}$, $0^{\circ}\mathrm{N}$) at

$-997.99\text{\,}\mathrm{m}$, and 2. 2.

$-1988.60\text{\,}\mathrm{m}$.

These are the depths that the $\omega$-surface heaved to from its initial depth at ($180^{\circ}\mathrm{E}$, $0^{\circ}\mathrm{N}$) of $-1000\text{\,}\mathrm{m}$ and $-2000\text{\,}\mathrm{m}$. All six surfaces are masked to exclude where any of them rise into the mixed layer101010 The mixed layer depth is taken as the depth, found by linear interpolation, at which potential density referenced to $100\text{\,}\mathrm{d}\mathrm{b}\mathrm{a}\mathrm{r}$ equals that at $-15\text{\,}\mathrm{m}$ (the second shallowest grid cell) plus $0.03\text{\,}\mathrm{k}\mathrm{g}\,\mathrm{m}^{-3}$ [Dong et al., 2008]. , where non-neutral dynamics take over. The modified topobaric surfaces differ slightly from regular topobaric surfaces; the area-weighted root-mean-square of the fictitious diapycnal diffusivity (McDougall and Jackett 2005; Klocker et al. 2009; S19) over the regular topobaric surfaces are (i) $1.16\text{\times}{10}^{-5}\text{\,}\mathrm{m}^{2}\,\mathrm{s}^{-1}$ and (ii) $7.34\text{\times}{10}^{-7}\text{\,}\mathrm{m}^{2}\,\mathrm{s}^{-1}$, whereas for modified topobaric surfaces these numbers are (i) $1.15\text{\times}{10}^{-5}\text{\,}\mathrm{m}^{2}\,\mathrm{s}^{-1}$ and (ii) $8.69\text{\times}{10}^{-7}\text{\,}\mathrm{m}^{2}\,\mathrm{s}^{-1}$.

On each of these six surfaces, the full velocity $\bm{u}$ is vertically interpolated onto each surface:

$\utilde{\bm{u}}$.

Then, the geostrophic velocity is estimated from five GSFs (all modified for the Boussinesq ocean), as well as from the $z$-level gradient of the pressure:

$\utilde{\bm{u}_{g}}^{(ZH)}=f^{-1}\bm{\hat{k}}\times\nabla\utilde{\Psi}^{(ZH)}$, using (4); 2. 3.

$\utilde{\bm{u}_{g}}^{(Cu)}=f^{-1}\bm{\hat{k}}\times\nabla\utilde{\Psi}^{(Cu)}$, using (22); 3. 4.

$\utilde{\bm{u}_{g}}^{(MK)}=f^{-1}\bm{\hat{k}}\times\nabla\utilde{\Psi}^{(MK)}$, using (4); 4. 5.

$\utilde{\bm{u}_{g}}^{(oM)}=f^{-1}\bm{\hat{k}}\times\nabla\utilde{\Psi}^{(oM)}$, using (5); 5. 6.

$\utilde{\bm{u}_{g}}^{(tb)}=f^{-1}\bm{\hat{k}}\times\nabla\utilde{\Psi}^{(tb)}$, using (2.6) while satisfying (14); 6. 7.

$\utilde{\bm{u}_{g}}^{(\nabla_{z}p)}=f^{-1}\bm{\hat{k}}\times\rho_{B}^{-1}\utilde{\nabla_{z}p}$, from (26).

For each velocity estimate $\bm{u}^{(i)}$ given by (a)–(g), the error

[TABLE]

is calculated, where, following (26),

[TABLE]

is taken as the “true” geostrophic velocity. Because the transformation from (2) to (3), i.e. transforming $\nabla_{z}p$ to $\nabla_{r}p$, is the first step for all GSFs, $\utilde{\bm{u}_{g}}$ is a more useful “truth” than $\utilde{\bm{u}_{g}}^{(\nabla_{z}p)}$. But $\utilde{\bm{u}_{g}}^{(\nabla_{z}p)}$ also has a fair claim to be the true geostrophic velocity, so the difference between $\utilde{\bm{u}_{g}}^{(\nabla_{z}p)}$ and $\utilde{\bm{u}_{g}}$ represents the precision with which we may know the actual geostrophic velocity. The actual geostrophic velocity is ambiguous below this precision, and there is little point reducing $\bm{\epsilon}$ below that for (g).

Two metrics of $\bm{\epsilon}$ will be calculated: the area-weighted $\ell_{1}$ and $\ell_{2}$ norms, respectively given by

[TABLE]

and

[TABLE]

where $\epsilon=[\epsilon^{(x)},\epsilon^{(y)}]$ is a 1D array concatenating the zonal and meridional errors within a given mask, and $A_{i}$ is the area of the grid cell at the location of $\epsilon_{i}$. The mask, which is common for all (a)–(g) on a given surface, excludes regions not connected to the main ocean, and excludes 1*∘* on either side of the equator where geostrophy is invalid (also recall the mixed layer has already been excluded). These metrics serve for comparing GSFs in Section 7.3, but first Section 7.2 maps just the zonal geostrophic velocity errors, $\epsilon^{(x)}$. (The meridional errors are qualitatively similar.) The metrics $\lVert\epsilon^{(x)}\rVert_{1}$ and $\lVert\epsilon^{(x)}\rVert_{2}$ are defined as in (31) and (32), but taking $\epsilon=[\epsilon^{(x)}]$ only.

7.2 Mapping geostrophic velocity errors

Figures 2–4 map $\epsilon^{(x)}$ on different surfaces. First consider Fig. 2(b)–(f), on the $\sigma_{1}$-surface. $\utilde{\Psi}^{(ZH)}$ and $\utilde{\Psi}^{(MK)}$ both perform well in the Indo-Pacific and the subtropical Atlantic, estimating $\utilde{\bm{u}_{g}}$ with accuracy generally below $\mathcal{O}($1\text{,}\mathrm{m}\mathrm{m},\mathrm{s}^{-1}$)$. Here, $\utilde{S}$, $\utilde{\theta}$, and $\utilde{z}$ are near their reference values, so errors are small — recall (4). In the North Atlantic/Arctic and Southern Ocean, $\utilde{\Psi}^{(ZH)}$ and (to a slightly lesser extent) $\utilde{\Psi}^{(MK)}$ perform poorly, estimating $\utilde{\bm{u}_{g}}$ with accuracy of $\mathcal{O}($10\text{,}\mathrm{m}\mathrm{m},\mathrm{s}^{-1}$)$, because $\utilde{\delta}$ varies considerably in these regions as the $\sigma_{1}$-surface shoals up to about $-100\text{\,}\mathrm{m}$ and to the sea surface, respectively. For $\utilde{\Psi}^{(Cu)}$, errors are small again in the Indo-Pacific and subtropical Atlantic where the surface remains within about $50\text{\,}\mathrm{m}$ of its reference depth of $-997.99\text{\,}\mathrm{m}$, but errors are large again in the North Atlantic/Arctic and Southern Ocean, both because of large gradients of $\utilde{S}$ and $\utilde{\theta}$ here and because the surface is far from the reference depth — recall (23). For $\utilde{\Psi}^{(oM)}$, the reference depth is a function of $\delta$, but the $\sigma_{1}$-surface is far from neutral, so this functional relationship is not very tight, evident by the significant scatter between $\utilde{z}$ and $\utilde{\delta}$ in Fig. 2$(*)$. Again, $\utilde{\Psi}^{(oM)}$ exhibits its largest errors in the North Atlantic and Southern Ocean, and to some degree in the Indian ocean, though it performs well in the Arctic; overall $\utilde{\Psi}^{(oM)}$ estimates $\utilde{\bm{u}_{g}}$ more accurately than do $\utilde{\Psi}^{(ZH)}$, $\utilde{\Psi}^{(Cu)}$, or $\utilde{\Psi}^{(MK)}$. Being geographically dependent, $\utilde{\Psi}^{(tb)}$ keeps errors low, globally, estimating $\utilde{\bm{u}_{g}}$ with accuracy generally below $\mathcal{O}($0.1\text{,}\mathrm{m}\mathrm{m},\mathrm{s}^{-1}$)$. Its highest errors are also found in the North Atlantic/Arctic, and Southern Ocean, where the $\sigma_{1}$-surface is furthest from neutral.

What level of geostrophic velocity error should be deemed acceptable? The ageostrophic zonal speed (Fig. 2a) is typically in the range $1\text{\,}\mathrm{m}\mathrm{m}\,\mathrm{s}^{-1}100\text{\,}\mathrm{m}\mathrm{m}\,\mathrm{s}^{-1}$; estimates of the geostrophic velocity should be at least as accurate as this. Most GSFs tested here pass this test in most places, but $\utilde{\Psi}^{(ZH)}$, $\utilde{\Psi}^{(Cu)}$, and $\utilde{\Psi}^{(MK)}$ fail this test in the North Atlantic/Arctic and Southern Ocean. On the other hand, it is unreasonable to wish the geostrophic velocity error be less than the difference between calculating the geostrophic velocity on $z$-levels or in the $\sigma_{1}$-surface (Fig. 2h). Only $\utilde{\Psi}^{(tb)}$ pushes close to this threshold of precision essentially everywhere.

Now, consider Fig. 3, on the upper $\tau^{\prime}$-surface. The errors in estimating the geostrophic velocity from all GSFs (b)–(f) are smaller on this $\tau^{\prime}$-surface than on the $\sigma_{1}$-surface, because the $\tau^{\prime}$-surface is closer to neutral. This reduces errors for $\utilde{\Psi}^{(oM)}$ and $\utilde{\Psi}^{(tb)}$ because their empirically fit functions become more accurate: notice the reduced scatter of $\utilde{z}$ vs. $\utilde{\delta}$ in Fig. 3$(*)$. In fact, $\tau^{\prime}$-surfaces make this relationship exact: there is actually zero scatter between $\utilde{z}$ vs. $\utilde{\delta}$, and the apparent scatter in Fig. 3$(*)$ is actually an inability to distinguish between thousands of regions (arcs of the Reeb graph) at the scale shown. As for $\utilde{\Psi}^{(ZH)}$ and $\utilde{\Psi}^{(MK)}$, the better neutrality of the $\tau^{\prime}$-surface tends to reduce $|\nabla\utilde{\delta}|$. To see this, expand \nabla\utilde{\delta}=\utilde{\nu_{S}}\nabla\utilde{S}+\utilde{\nu_{\theta}}\nabla\utilde{\theta}+\big{(}\utilde{\nu_{p}}-\partial_{p}\mathcal{V}(S_{0},\theta_{0},\utilde{p})\big{)}\nabla\utilde{p} (reverting to the non-Boussinesq form for familiarity), then note the first two terms are close to zero for nearly neutral surfaces, by (6).

Finally, consider Fig. 4, on the lower $\omega$-surface. Again, $\utilde{\Psi}^{(ZH)}$, $\utilde{\Psi}^{(Cu)}$, and $\utilde{\Psi}^{(MK)}$ perform well in the Pacific where the surface properties remain near the reference properties. In the subtropical Atlantic though, $\utilde{z}$ is about $250\text{\,}\mathrm{m}$ shallower than in the Pacific, and this depth difference increases errors for $\utilde{\Psi}^{(Cu)}$ over most of the subtropical Atlantic. Also, this surface samples North Atlantic Deep Water, so $\utilde{S}$ and $\utilde{\theta}$, and hence $\utilde{\delta}$, in the Atlantic differ substantially from their Pacific reference values, which causes considerable errors for $\utilde{\Psi}^{(ZH)}$ across the entire Atlantic; $\utilde{\Psi}^{(MK)}$ is similarly affected but to a lesser degree. Once again, $\utilde{\Psi}^{(ZH)}$, $\utilde{\Psi}^{(Cu)}$, and $\utilde{\Psi}^{(MK)}$ exhibit large errors in the Southern Ocean, particularly near the outcrop where the surface is farthest from its reference depth, and the surface properties furthest from their reference values. As $\omega$-surfaces are extremely neutral, $\utilde{z}$ is tightly related to $\utilde{\delta}$, so $\utilde{\Psi}^{(oM)}$ keeps its reference depth close to $\utilde{z}$, and $\utilde{\Psi}^{(oM)}$ performs well, globally. Again, $\utilde{\Psi}^{(tb)}$ produces the best estimate of $\utilde{\bm{u}_{g}}$, keeping errors low, globally.

For most GSFs, $\bm{\epsilon}$ is high where the geostrophic velocity itself is high, particularly in the Southern Ocean and North Atlantic, and to a lesser degree in the Kuroshio current and near the equator. We could instead study the relative error (divide $\bm{\epsilon}$ by $|\utilde{\bm{u}_{g}}|$), but this creates the opposite problem: errors tend to be high where $\utilde{\bm{u}_{g}}$ is near zero. Nonetheless, the results are qualitatively the same for relative and absolute errors (not shown), so we proceed with the latter.

7.3 Quantitative error measurements

We now turn to more quantitative comparisons of the geostrophic velocity errors. Figure 5 shows $\lVert\epsilon\rVert_{1}$ and $\lVert\epsilon\rVert_{2}$, as well as boxplots providing the area-weighted 0.05, 0.25, 0.5, 0.75, and 0.95 quantiles of $\epsilon$ for the seven velocity estimates on the three surface classes at both depths. For numeric clarity, $\lVert\epsilon\rVert_{1}$ and $\lVert\epsilon\rVert_{2}$ are also listed in Table 1. The following discussion focuses mostly on $\lVert\epsilon\rVert_{2}$.

For $\utilde{\Psi}^{(ZH)}$, $\lVert\epsilon\rVert_{2}$ is about $1\text{\,}\mathrm{m}\mathrm{m}\,\mathrm{s}^{-1}1.5\text{\,}\mathrm{m}\mathrm{m}\,\mathrm{s}^{-1}$ on the $\tau^{\prime}$- and $\omega$-surfaces, but about $7\text{\,}\mathrm{m}\mathrm{m}\,\mathrm{s}^{-1}35\text{\,}\mathrm{m}\mathrm{m}\,\mathrm{s}^{-1}$ on the $\sigma_{1}$- and $\sigma_{2}$-surfaces. These errors are reduced for $\utilde{\Psi}^{(MK)}$ by a factor of about 2.5 on the $\sigma_{1}$- and $\sigma_{2}$-surfaces, whereas on the $\tau^{\prime}$- and $\omega$-surfaces this factor is about $34$ on the upper surfaces and about 1 (no reduction) on the lower surfaces. For the orthobaric Montgomery potential $\utilde{\Psi}^{(oM)}$ this factor is roughly $3.55.5$ on the upper surfaces, and $729$ on the lower surfaces. For the topobaric geostrophic streamfunction, $\utilde{\Psi}^{(tb)}$, this reduction factor is in the range $1745$ on the non-$\tau^{\prime}$-surfaces (it is hardly fair to compare against $\utilde{\Psi}^{(tb)}$ on the $\tau^{\prime}$-surface where it is exact). Indeed, $\utilde{\Psi}^{(tb)}$ improves upon previous results across all measures shown in Fig. 5.

The worst performance by $\utilde{\Psi}^{(tb)}$ is on the $\sigma_{1}$- and $\sigma_{2}$-surfaces, which are furthest from neutral. On $\omega$-surfaces, $\utilde{\Psi}^{(tb)}$ estimates $\utilde{\bm{u}_{g}}$ about as well as does $\utilde{\bm{u}_{g}}^{(\nabla_{z}p)}$. On the $\tau^{\prime}$-surfaces, $\utilde{\bm{u}_{g}}^{(tb)}$ exhibits $\lVert\epsilon\rVert_{1}$ and $\lVert\epsilon\rVert_{2}$ about an order of magnitude smaller than those for $\utilde{\bm{u}_{g}}^{(\nabla_{z}p)}$. In fact, $\utilde{\Psi}^{(tb)}$ is exact on $\tau^{\prime}$-surfaces, and $\utilde{\bm{u}_{g}}^{(tb)}-\utilde{\bm{u}_{g}}$ reveals errors not in $\utilde{\bm{u}_{g}}^{(tb)}$ but in $\utilde{\bm{u}_{g}}$. Specifically, $\utilde{\bm{u}_{g}}$ involves a finite difference of $\utilde{z}$ multiplied by a simple average of $\utilde{\rho}$, whereas $\utilde{\bm{u}_{g}}^{(tb)}$ involves a finite difference of two integrals of $\hat{\rho}$ with respect to depth. These would be equivalent if $\hat{\rho}$ were an affine linear function of $\utilde{z}$, but this is false for two reasons. First, $\hat{\rho}$ is multivalued, and this matters when the finite difference selects points in different domains of the branches of $\hat{\rho}$. Second, each branch of $\hat{\rho}$ is, by our choice (16), a quadratic function of depth (from integrating an empirically fit affine linear function) plus a rational function of depth (from the equation of state).

Table 2 lists the percentage of surface area within which the estimated geostrophic velocity exceeds (in magnitude) the ageostrophic velocity. For $\utilde{\Psi}^{(ZH)}$, this area occupies almost 7% of the $\sigma_{1}$-surface, while $\utilde{\Psi}^{(Cu)}$ and $\utilde{\Psi}^{(MK)}$ reduce this area to about 3%, and $\utilde{\Psi}^{(oM)}$ reduces it to about 2%. For $\utilde{\Psi}^{(tb)}$, this area is a mere 0.3%. Results are similar on the other five surfaces. It is therefore important to use an accurate GSF, particularly when working on global surfaces, lest the ageostrophic velocity be swamped by errors in estimating the geostrophic velocity.

From a numerical perspective, the trouble caused by the ill-defined nature of the exact GSF for neutral surfaces tends to be rather small. For example, on the $\omega$-surface intersecting ($180^{\circ}\mathrm{E}$, $0^{\circ}\mathrm{N}$, $-1988.60\text{\,}\mathrm{m}$), two topobaric GSFs are estimated using the methods of Section 3b but excluding the cycle constraints (14). They differ only in that they integrate $\hat{\nu}_{a}$ differently on each arc $a$ that defines a cycle in the cycle basis: one integrates up from $p_{l_{a}}$ and the other down from $p_{h_{a}}$. The geostrophic velocity estimated by these two GSFs differs at 41,140 grid points, out of a total of 868,954 grid points on this surface (4.7%). The (unweighted) $\ell_{2}$ norm of the velocity difference over the grid points where they differ at all is a mere $3.6\text{\times}{10}^{-4}\text{\,}\mathrm{m}\,\mathrm{s}^{-1}$.

Finally, the orthobaric Montgomery potential is indeed more accurate than the orthobaric GSF. On the $\omega$-surface through ($180^{\circ}\mathrm{E}$, $0^{\circ}\mathrm{N}$, $-1988.60\text{\,}\mathrm{m}$), $\lVert\epsilon\rVert_{2}=$1.4\text{\times}{10}^{-4}\text{,}\mathrm{m},\mathrm{s}^{-1} for $\utilde{\Psi}^{(oM)}$, whereas for the orthobaric GSF (fitting $\utilde{\delta}$ to $\utilde{p}$ as a cubic spline with 12 pieces), $\lVert\epsilon\rVert_{2}=$2.7\text{\times}{10}^{-4}\text{\,}\mathrm{m}\,\mathrm{s}^{-1}. The utility of the orthobaric Montgomery potential becomes more apparent for surfaces that include both the Southern Ocean and the Arctic: on the $\omega$-surface through ($180^{\circ}\mathrm{E}$, $0^{\circ}\mathrm{N}$, $-997.99\text{\,}\mathrm{m}$), $\lVert\epsilon\rVert_{2}=$3.8\text{\times}{10}^{-4}\text{,}\mathrm{m},\mathrm{s}^{-1} for $\utilde{\Psi}^{(oM)}$, whereas for the orthobaric GSF, $\lVert\epsilon\rVert_{2}=$9.2\text{\times}{10}^{-4}\text{\,}\mathrm{m}\,\mathrm{s}^{-1}.

8 Summary

The exact GSF on a neutral surface, whose existence has long been known [McDougall, 1989], has been derived. It is defined using path integrals, along neutral trajectories, of the specific volume as a function of pressure. On a (hypothetical) well-defined neutral surface, the specific volume is a multivalued function of the pressure, and its geographic structure is described by the Reeb [1946] graph of the pressure on the surface [S19]. The path integrals defining the GSF are equivalently described by a sum of simple integrals determined by a walk through the Reeb graph. Islands, and other holes in the neutral surface, can create cycles in the Reeb graph that may be walked in either direction. That is, they cause the path integral to be path-dependent, and the exact GSF for a neutral surface is actually a multivalued function of geographic position. When mapping the GSF, its multivalued nature appears as discontinuities emanating from islands and other such holes. Though the gradient of the GSF, and hence the geostrophic velocity, is unaffected by this multivalued nature, it does hamper the use of streamlines for flow visualization or other analyses.

The problem of these discontinuities is overcome by the topobaric GSF, which everywhere approximates the exact GSF for neutral surfaces, but is well-defined. Moreover, the topobaric GSF can be used on any approximately neutral surface. On a general approximately neutral surface, the specific volume is empirically fit as a multivalued function of pressure, subject to a set of constraints determined by cycles in the Reeb graph, that ensure the topobaric GSF is well-defined. Numerical tests reveal that the topobaric GSF estimates the geostrophic velocity on potential density surfaces and $\omega$-surfaces far more accurately than any other known GSF. On topobaric surfaces, the specific volume is given a priori as a multivalued function of pressure. Numerical tests confirm that the topobaric GSF is exact on modified topobaric surfaces.

The Montgomery [1937] potential has also been generalized, furthering the progress made by Zhang and Hogg [1992]. This leads to an alternative formulation of the exact GSF for neutral surfaces, which is amenable to an approximation wherein the pressure on a surface is empirically fit as a single-valued function of the specific volume anomaly on the surface. This approximation, called the orthobaric Montgomery potential, is easy to compute and improves upon all previously known GSFs for approximately neutral surfaces.

Acknowledgements

The author thanks David Marshall for helpful discussions, Trevor McDougall for early encouragement to pursue this topic, two reviewers for their constructive feedback, and authors of quality software, including Andreas Klocker ($\omega$-surfaces), David Gleich (GAIMC), Jonas Lundgren (splinefit), and Harish Doraiswamy and Vijay Natarajan (ReCon). GJS was supported by the Clarendon Scholarship, and the Canadian Alumni Scholarship at Linacre College, University of Oxford. MATLAB software to compute the topobaric geostrophic streamfunction, the orthobaric Montgomery potential, and other results herein is available from the author’s website.

Appendix A Boussinesq geostrophic streamfunctions

As discussed in Section 6, the Boussinesq approximation swaps the roles of $p$ and $z$, and of $\nu$ and $\rho$. It is helpful to re-define $\delta$ as the in-situ density anomaly,

[TABLE]

where $\rho_{0}=B(S_{0},\theta_{0},z)$ is the in-situ density at the local depth but at a reference salinity $S_{0}$ and reference potential temperature $\theta_{0}$. The derivations of Sections 2–5 proceed, in essence, unchanged. The results are as follows.

Hydrostatic balance is integrated to obtain the pressure,

[TABLE]

The exact GSF for a neutral surface (8) becomes,

[TABLE]

and the graph integral form (2.6) transforms similarly.

For the Montgomery potential, (18) becomes

[TABLE]

for some constant depth $Z$.

For the Zhang and Hogg [1992] GSF, (4) becomes

[TABLE]

for another constant $z_{0}$.

For the orthobaric Montgomery potential, (5) becomes

[TABLE]

where $\hat{z}$ is a function of $\delta$.

For the Cunningham [2000] GSF, (22) becomes

[TABLE]

For the McDougall and Klocker [2010] GSF, the Boussinesq form is most easily found by noting how terms transform in the Montgomery potential: $g\utilde{z}$ transforms to $\rho_{B}^{-1}\utilde{p}$, and $\utilde{p}(\utilde{\nu}-\nu_{0})$ to $g\rho_{B}^{-1}\utilde{z}(\utilde{\rho}-\rho_{0})$. One can use these rules to determine how the $T_{b}/(12\,\utilde{\rho})$ term in $\utilde{\Psi}^{(MK)}$ transforms to the Boussinesq form, pivoting on their Eq. 55. The Boussinesq version of (4) is

[TABLE]

where again $T_{b}/\utilde{\rho}=$2.7\text{\times}{10}^{-15}\text{,}\mathrm{K}^{-1}\mathrm{P}\mathrm{a}^{-2}\mathrm{m}^{2}\mathrm{s}^{-2}$$ is treated as constant.

Appendix B Numerical methods

The numerical methods to calculate the various GSFs from discretised model data require some care. They depend, first, on how the geostrophic velocity is calculated. The model has discrete depth levels $z_{k}<0$ for the centre of each tracer cell, where the salinity $S_{i,j,k}$ and potential temperature $\theta_{i,j,k}$, and thus the in-situ density $\rho_{i,j,k}$, are known in each water column $(i,j)$. For each water column, the model calculates the pressure $p_{i,j,k}$ at each $z_{k}$ by integrating hydrostatic balance, (34), using $\rho_{B}$ as the density between $z=0$ and the free surface, $\rho_{i,j,1}$ as the density between $z=0$ and $z=z_{1}$, and otherwise using trapezoidal integration:

[TABLE]

This must be extended to give the pressure $p_{i,j}(z^{\prime})$ at an arbitrary depth $z^{\prime}$. Linearly interpolating the pressure implies, via hydrostatic balance, that the in-situ density is uniformly $(\rho_{i,j,k}+\rho_{i,j,k+1})/2$ between $z_{k}$ and $z_{k+1}$. Then $\utilde{\rho}$ would discontinuously jump where the surface crosses grid points — most unconscionable — and so too would the geostrophic velocity given by (30). Instead, the in-situ density is continuously extended to $\rho_{i,j}(z^{\prime})$, then hydrostatic balance gives the pressure as

[TABLE]

where $k$ is such that $z_{k}\leq z^{\prime}<z_{k+1}$. This form ensures the pressure is continuous, i.e. $p_{i,j}(z_{k})=p_{i,j,k}$.

Now suppose we define $\rho_{i,j}(z^{\prime})$ by linear interpolation:

[TABLE]

The discretised version of the middle equation in (26) — rather, its meridional component, up to a constant scaling factor, and suppressing zonal grid indices for brevity — is

[TABLE]

where $\bar{z}_{j}=(\utilde{z}_{j-1}+\utilde{z}_{j})/2$. In (44), $\utilde{z}_{j-1}$ and $\utilde{z}_{j}$ are assumed to lie in the same range, $[z_{k},z_{k+1}]$. (If they are not, the error term differs, but no discontinuities appear.) Dividing (44) by the meridional distance between water columns, the first line is the discretised $z$-level pressure gradient, the second line is the discretised in-surface pressure gradient plus a contribution that undoes the sloping gradient, and the third line is a third-order error. Zonal discretisations are analogous.

Now, neutral errors on a surface are traditionally calculated using $\nabla\utilde{S}$ and $\nabla\utilde{\theta}$, which has an equivalent form using $\nabla\utilde{\rho}$ and $\nabla\utilde{z}$ only if $\utilde{\rho}=B(\utilde{S},\utilde{\theta},\utilde{z})$ — see the discussion around Eq. 28 of S19. For consistency with this, instead of (43), we define $\rho_{i,j}(z^{\prime})=B(S_{i,j}(z^{\prime}),\theta_{i,j}(z^{\prime}),z^{\prime})$, then define $S_{i,j}(z^{\prime})$ and $\theta_{i,j}(z^{\prime})$ by linear interpolation analogous to (43). This introduces extra errors to (44) due to non-linearity in $B$, but they are considerably smaller than the discretization error, the third line in (44).

Whereas the integral of in-situ density is designed to match the model’s pressure, the integral of the reference in-situ density has the sole purpose of satisfying

[TABLE]

To best satisfy this numerically, $\rho_{0}$ is trapezoidally integrated using a mesh with $0.1\text{\,}\mathrm{m}$ spacing, much finer than the model’s $z_{k}$ levels. (Where the reference $S$ and $\theta$ are not constants, as in $\utilde{\Psi}^{(Cu)}$, the model’s $z_{k}$ levels are used.)

The remaining integrals, namely the affine linear part of $\utilde{\Psi}^{(tb)}$, and the $\hat{p}$ part of $\utilde{\Psi}^{(oM)}$, are handled analytically.

The TEOS-10 routines to calculate $\utilde{\Psi}^{(ZH)}$ and $\utilde{\Psi}^{(MK)}$ are not used, for they assume a too-recent equation of state and a non-Boussinesq ocean, whereas ECCO2 uses the Boussinesq form of the Jackett and Mcdougall [1995] equation of state. Instead, $\utilde{\Psi}^{(ZH)}$ and $\utilde{\Psi}^{(MK)}$ are calculated from scratch, using the aforementioned vertical interpolation and integration methods.

ECCO2 uses a C-grid, so the zonal velocity $u$ lives on the west face of a tracer cell, whereas the zonal geostrophic velocity lives on the south face of a tracer cell. To evaluate $\utilde{u}$ — the zonal component of (a) in Section 7.1 — the 3D $u$ field is brought to the south face by the same weighted average that the model uses when calculating the Coriolis acceleration, then linearly interpolated to $\bar{z}_{i,j}$. The meridional velocity is handled analogously.