Neutral surface topology

TL;DR

This paper develops a topological framework for understanding and approximating neutral surfaces in the ocean, introducing topobaric surfaces that are neutral, computationally efficient, and topologically consistent.

Contribution

It introduces a topological theory of neutral surfaces using Reeb graphs and proposes topobaric surfaces as a new class of approximately neutral, computationally efficient surfaces.

Findings

Topobaric surfaces closely approximate neutral surfaces in the ocean.

Reeb graph topology effectively describes neutral surface structure.

Topobaric surfaces are faster to compute than traditional methods.

Abstract

Neutral surfaces, along which most of the mixing in the ocean occurs, are notoriously difficult objects: they do not exist as well-defined surfaces, and as such can only be approximated. In a hypothetical ocean where neutral surfaces are well-defined, the in-situ density on the surface is a multivalued function of the pressure on the surface, $\underset{\sim}{p}$. The surface is decomposed into geographic regions where there is one connected pressure contour per pressure value, making this function single-valued in each region. The regions are represented by arcs of the Reeb graph of $\underset{\sim}{p}$. The regions meet at saddles of $\underset{\sim}{p}$ which are represented by internal nodes of the Reeb graph. Leaf nodes represent extrema of $\underset{\sim}{p}$. Cycles in the Reeb graph are created by islands and other holes in the neutral surface. This topological theory of neutral surfaces is used to create a new class of approximately neutral surfaces in the real ocean, called topobaric surfaces, which are very close to neutral and fast to compute. Topobaric surfaces are the topologically correct extension of orthobaric density surfaces to be geographically dependent, which is fundamental to neutral surfaces. Also considered is the possibility that helical neutral trajectories might have a larger pitch around islands than in the open ocean.

Equations

$$\nabla \rho ={\rho }_S\nabla S+{\rho }_{\theta }\nabla \theta +{\rho }_p\nabla p,$$

$$\mathbf{N}={\rho }_S\nabla S+{\rho }_{\theta }\nabla \theta ,$$

$${\int }_{\partial \Omega }\mathbf{N}\cdot \mathrm{d}\mathbf{r}={\int }_{\Omega }\nabla \times \mathbf{N}\cdot \mathrm{d}\mathbf{\Omega }.$$

$$0={\int }_{\Omega }\frac{H}{|\mathbf{N}|}\mathrm{d}\Omega ,$$

$$H=({\rho }_{\theta }{\rho }_{Sp}-{\rho }_S{\rho }_{\theta p})\nabla p\cdot \nabla S\times \nabla \theta ,$$

$$0=\mathrm{d}\mathbf{r}\cdot \mathbf{N}=\mathrm{d}\mathbf{r}\cdot ({\rho }_S{\nabla }_{N}S+{\rho }_{\theta }{\nabla }_N\theta ),$$

$${\rho }_S{\nabla }_{N}S+{\rho }_{\theta }{\nabla }_N\theta =0.$$

$$\underset{˜}{{\rho }_S}\nabla \underset{˜}{S}+\underset{˜}{{\rho }_{\theta }}\nabla \underset{˜}{\theta }=0.$$

$$\nabla \underset{˜}{\rho }=\underset{˜}{{\rho }_p}\nabla \underset{˜}{p},$$

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

$$\nabla \underset{˜}{\rho }=\frac{\mathrm{d}\hat{\rho }}{\mathrm{d}p}(\underset{˜}{p})\nabla \underset{˜}{p}.$$

$$\frac{\mathrm{d}\hat{\rho }}{\mathrm{d}p}(\underset{˜}{p})=\underset{˜}{{\rho }_p}.$$

$$\underset{˜}{{\rho }_p}=\hat{\pi }(\underset{˜}{p}),$$

$$0=\hat{\mathbf{k}}\cdot \nabla \underset{˜}{{\rho }_p}\times \nabla \underset{˜}{p}.$$

$$\hat{\rho }(\mathtt{p})={\rho }_c+{\int }_{p_c}^{\mathtt{p}}\hat{\pi }(p^{\prime })\mathrm{d}{p}^{\prime },$$

$$\utilde{\rho}(\bm{x})=\hat{\rho}_{a}\big{(}\utilde{p}(\bm{x})\big{)}\quad\forall\bm{x}\in\mathcal{R}_{a}.$$

$$0=\left(\underset{˜}{{\rho }_{p\theta }}-\underset{˜}{{\rho }_{pS}}\underset{˜}{{\rho }_{\theta }}{\underset{˜}{{\rho }_S}}^{-1}\right)\hat{\mathbf{k}}\cdot \nabla \underset{˜}{\theta }\times \nabla \underset{˜}{p}.$$

$$\utilde{\rho_{p}}(\bm{x})=\hat{\pi}_{a}\big{(}\utilde{p}(\bm{x})\big{)}\quad\forall\bm{x}\in\mathcal{R}_{a},$$

$$\underset{˜}{\rho }(\mathbf{x})-\underset{˜}{\rho }({\mathbf{x}}_0)={\int }_{\mathcal{P}}\nabla \underset{˜}{\rho }\cdot \mathrm{d}\mathbf{r}={\int }_{\mathcal{P}}\underset{˜}{{\rho }_p}\nabla \underset{˜}{p}\cdot \mathrm{d}\mathbf{r}.$$

$${\hat{\rho }}_{a_{\mathcal{J}}}(\mathtt{p})={\hat{\rho }}_{a_0}(p_0)$$

$$+$$

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

$${\hat{\pi }}_a(p)=K_a+L_a(p-p_{{\ell }_a})+{\partial }_{p}R(S_0,{\theta }_0,p),$$

$${\hat{\rho }}_a(p)=J_a+K_a(p-p_{{\ell }_a})+\frac{L_a}{2}(p-p_{{\ell }_a}{)}^2+R(S_0,{\theta }_0,p).$$

$$R\Big{(}\mathcal{S}(\bm{x},p^{\prime}),\ \vartheta(\bm{x},p^{\prime}),\ p^{\prime}\Big{)}=\hat{\rho}_{a}(p^{\prime}),$$

$$\rho =B(S,\theta ,z)=R(S,\theta ,-g{\rho }_{B}z),$$

$$\nabla \rho ={\rho }_S\nabla S+{\rho }_{\theta }\nabla \theta +{\rho }_z\nabla z,$$

$$\nabla \underset{˜}{\rho }=\underset{˜}{{\rho }_z}\nabla \underset{˜}{z}.$$

$$\mathbf{ϵ}=\nabla \underset{˜}{\rho }-\underset{˜}{{\rho }_z}\nabla \underset{˜}{z},$$

$$D^f=K\mathbf{s}\cdot \mathbf{s}$$

Full text

††journal: Ocean Modelling

1 Introduction

Strong stratification throughout most of the ocean inhibits vertical motion, largely confining the oceanic flow to a two-dimensional surface, the ideal of which is called a neutral surface [McDougall, 1987a]. These surfaces are far from flat, and it is along these sloping surfaces that oceanic flows efficiently mix tracers (epineutral mixing), whereas tracer mixing across them (dianeutral mixing) is enormously weaker—an idea tracing back to Iselin [1939]. This is a great conceptual simplification, but only useful if we can map the depth, or pressure, of such surfaces. Unfortunately, non-linearity in seawater’s equation of state leads to a path-dependence underlying the definition of neutral surfaces, making neutral surfaces ill-defined [McDougall and Jackett, 1988].

Given this difficulty, physical oceanographers craft well-defined surfaces that are everywhere nearly tangent to the neutral tangent plane, called approximately neutral surfaces. These surfaces are usually isosurfaces of a 3D variable, the earliest being potential density [Wüst, 1935] and specific volume anomaly [Montgomery, 1937]. Lynn and Reid [1968] revealed the highly undesirable property that, far away from its reference pressure, potential density surfaces (isopycnals) in a stably stratified ocean can exhibit unphysical overturns. This problem also affects specific volume anomaly surfaces, far from the reference values.

To overcome this, Reid and Lynn [1971] introduced patched potential density. They map the $\sigma_{4}=45.92$ potential density surface (referenced to $4000\text{\,}\mathrm{d}\mathrm{b}\mathrm{a}\mathrm{r}$) in the tropical Atlantic, and where this surface rises above $-3000\text{\,}\mathrm{m}$ in the North Atlantic, it is patched together with the $\sigma_{2}=37.14$ potential density surface (referenced to $2000\text{\,}\mathrm{d}\mathrm{b}\mathrm{a}\mathrm{r}$). In fact, $\sigma_{2}$ varies somewhat along the length of this $-3000\text{\,}\mathrm{m}$ contour, and 37.14 is chosen to minimize this discontinuity. Similarly, where the $\sigma_{4}=45.92$ surface rises above $-3000\text{\,}\mathrm{m}$ in the South Atlantic, they patch it together with the $\sigma_{2}=37.10$ surface. Noting that a single $\sigma_{4}$ surface is patched together with different $\sigma_{2}$ surfaces in the North Atlantic and Southern Ocean, it is clear that neutral surfaces are not just dependent upon salinity, temperature, and pressure, but also upon geography (latitude and longitude).

In this way, the ocean may be cut, stacked, and arranged into boxes covering certain depth ranges and horizontal areas. Where to make these cuts is not entirely arbitrary. Indeed, the equatorial Atlantic is a good place to make a cut, where a “bridge” region connects saltier North Atlantic waters with fresher South Atlantic waters; salinity and potential temperature are tightly related in these three regions, but the functional relationship differs between regions [de Szoeke and Springer, 2005].

Though not entirely arbitrary, these divisions are also not entirely correct. There is nothing particularly important, thermodynamically, about the equator, nor about any other latitude circle, longitude circle, or depth level. There are thermodynamically important regions, but their shape is not so simple. To study their shape is to study topology.

This paper presents a fresh theoretical perspective on this problem, by studying the topology of hypothetical neutral surfaces that are well-defined. On such surfaces, there is a multivalued functional relationship between the in-situ density and the pressure. Different branches of this multivalued function arise because a level set of pressure on a neutral surface can be the disjoint union of multiple connected components, each of which supports a distinct in-situ density. The important topological information about changes in the connectedness of these level sets is captured, as a collection of nodes and arcs, by the Reeb [1946] graph. This also determines the shape of the regions inside which the aforementioned multivalued function actually is just single-valued. This theoretical tool is then used to develop a new class of approximately neutral surfaces in the real ocean, called topobaric surfaces. Topobaric surfaces are very close to neutral and possess an exact geostrophic streamfunction [Stanley, 2019]. Moreover, they are fast to compute: Computational topology is a young field, but efficient algorithms to compute the Reeb graph have recently been developed [Doraiswamy and Natarajan, 2013].

Though defined over 70 years ago, the Reeb graph has not, to the best of the author’s knowledge, been previously used in oceanography, nor as a way of studying multivalued functional relationships between variables. The Reeb graph was most famously used by Arnol’d [1957] in solving Hilbert’s superposition problem [see also Arnold, 2006], but its primary use of late is in computer graphics and visualization [see Biasotti et al., 2008, for a review].

The paper is structured as follows. The theory of neutral surfaces is reviewed in Section 2, then developed in Section 3 from a topological perspective, discussing the Reeb graph, as well as the role of islands in making neutral surfaces ill-defined. A pedagogical illustration of the multivalued functional relationship and the Reeb graph is given in Section 4. Section 5 discusses topobaric surfaces, from their theoretical description to their numerical calculation, and finally to their evaluation as useful approximately neutral surfaces. Conclusions are given in Section 6. A glossary of graph theory definitions is given in A, and B contains a preliminary analysis of the role of islands in helical neutral motions.

2 Background of neutral surfaces

2.1 Definitions

In the ocean, the salinity $S$, potential temperature $\theta$, and pressure $p$ determine the in-situ density $\rho$ according to a function $R$, the equation of state.222 Everything could be described in terms of Absolute Salinity and Conservative Temperature, and this would be completely equivalent for the theory presented here. This paper speaks of practical salinity and potential temperature simply because these are the inputs to the Jackett and Mcdougall [1995] equation of state that is used by the ocean model whose data shall be analysed. Mathematically, $\rho=R(S,\theta,p)$, where $S$, $\theta$, $p$, and $\rho$ are 3D scalar fields. Using the chain rule, the gradient of in-situ density is

[TABLE]

where $\rho_{S}=\partial_{S}R(S,\theta,p)$, $\rho_{\theta}=\partial_{\theta}R(S,\theta,p)$, and $\rho_{p}=\partial_{p}R(S,\theta,p)$ are 3D scalar fields.333 Often (1) is divided by $\rho$, so as to use the thermal expansion coefficient $-\rho^{-1}\rho_{\theta}$, the haline contraction coefficient $\rho^{-1}\rho_{S}$, and the adiabatic compressibility $\rho^{-1}\rho_{p}$. This complicates further differentiation and so is not used here.

Let $\bm{\hat{i}}$, $\bm{\hat{j}}$, and $\bm{\hat{k}}$ be the eastward, northward, and radial (vertical) unit vectors, respectively.

Consider displacing a fluid parcel of in-situ density $\rho_{0}$ infinitesimally by $\mathrm{d}\bm{r}$. Its new surroundings have in-situ density $\rho_{0}+\nabla\rho\cdot\mathrm{d}\bm{r}$. Using (1), this has contributions from salt, potential temperature, and pressure changes. If the fluid parcel is insulated (meaning the salinity and potential temperature are conserved, commonly called adiabatic), its in-situ density after displacement is only $\rho_{0}+\rho_{p}\nabla p\cdot\mathrm{d}\bm{r}$. The difference between the environmental and parcel in-situ density, $(\rho_{S}\nabla S+\rho_{\theta}\nabla\theta)\cdot\mathrm{d}\bm{r}$, creates a buoyant restoring force. Thus, if the displacement $\mathrm{d}\bm{r}$ is perpendicular to the dianeutral vector

[TABLE]

then the buoyant restoring force is zero. The neutral tangent plane is the plane normal to the dianeutral vector [McDougall, 1987a]. A fluid parcel can move infinitesimally in this plane without experiencing a buoyant restoring force.444 Noting that mixing along the neutral tangent plane does require some kinetic energy, Nycander [2011] defined a vector $\bm{P}$ using the dynamic enthalpy rather than the equation of state. The form of $\bm{P}$ is the same as $\bm{N}$, so the essential ideas of this manuscript will apply equally well to surfaces formed from $\bm{N}$ or $\bm{P}$.

A neutral trajectory is a solution of the Pfaffian differential equation $\bm{N}\cdot\,\mathrm{d}\bm{r}=0$. That is, a neutral trajectory is a path that is always orthogonal to $\bm{N}$ and hence always tangent to the local neutral tangent plane.

A neutral surface is a surface that is everywhere tangent to the neutral tangent plane.

2.2 Non-existence of neutral surfaces

Following McDougall and Jackett [1988], suppose we have a neutral surface that is a well-defined, 2D surface. Consider a subset $\Omega$ of the neutral surface, having boundary $\partial\Omega$, which is a neutral trajectory. The path integral of $\bm{N}$ along $\partial\Omega$ is related to an area integral of the curl of $\bm{N}$ via Stokes’ theorem:

[TABLE]

The LHS is zero by the definition of a neutral trajectory. For the RHS, the surface of integration is a neutral surface so $\mathrm{d}\bm{\Omega}=\mathrm{d}\Omega\,\bm{N}/|\bm{N}|$, and hence

[TABLE]

where $H=\bm{N}\cdot\nabla\times\bm{N}$ is the neutral helicity. As $\Omega$ is arbitrary, it must be that $H=0$ everywhere. Expanding $\bm{N}$ and using the chain rule on $\nabla\rho_{S}$ and $\nabla\rho_{\theta}$ akin to (1), we find

[TABLE]

where $\rho_{Sp}=\partial_{Sp}R(S,\theta,p)$ and $\rho_{\theta p}=\partial_{\theta p}R(S,\theta,p)$ are 3D fields. With a non-linear equation of state $R$, the term in parentheses is non-zero; it is $(\rho_{S}\rho)$ times the thermobaricity [defined by McDougall, 1987b]. Thus $H=0$ if and only if the three vectors $\nabla p$, $\nabla S$, and $\nabla\theta$ are coplanar. As this is not generally true [see McDougall and Jackett, 2007, for analysis of oceanic data], a contradiction is reached. The false assumption was to assume the neutral surface was a well-defined, 2D surface. Even if the ocean is stably stratified ($\bm{\hat{k}}\cdot\bm{N}<0$ everywhere), the neutral trajectory $\partial\Omega$ has returned to its initial geographic location but at a different depth from which it began. This phenomenon is called the neutral helix, and this depth change is called the pitch [McDougall and Jackett, 1988].

A neutral helix can be arbitrarily shrunk, in terms of its radially projected area, to produce a new neutral helix with a different pitch.555If the ocean is neutrally stable somewhere, the neutral tangent plane there contains the radial (vertical) direction. A better version of this argument is to measure the pitch in the $\bm{N}$ direction. This new pitch is a continuous function of the factor by which this area is shrunk, so any desired pitch, between the original pitch and zero, can be found: this is guaranteed by the intermediate value theorem. Thus, neutral trajectories from a given point are a set of points that extend laterally, along the neutral tangent planes, as well as radially (vertically), and so occupy a 3D volume rather than a 2D surface. This is how a neutral surface in an $H=0$ ocean fails to be a well-defined, 2D surface.

2.3 Well-defined surfaces vs. unique-depth surfaces

We have just seen that a necessary condition for a neutral surface to be a well-defined surface is that $H=0$ everywhere on that surface. Jackett and McDougall [1997] asserted that $H=0$ globally is also a sufficient condition to ensure neutral surfaces are well-defined surfaces, and this thinking has persisted [McDougall and Jackett, 2007].

This is correct, if one takes a well-defined surface by the standard mathematical and topological definition: a well-defined surface is a 2-manifold. Roughly speaking, this means it locally resembles 2D Euclidean space: every point on a 2-manifold has a neighbourhood that can be continuously deformed into an open subset of $\mathbb{R}^{2}$, and back again. In a stably stratified ocean, however, we would ideally like neutral surfaces to be unique-depth surfaces: well-defined surfaces whose depth is a single-valued function of geographic location (though undefined where the surface has grounded or outcropped).

In fact, $H=0$ globally is not sufficient to ensure neutral surfaces are unique-depth surfaces in a stably stratified ocean. The aforementioned assertion by Jackett and McDougall [1997] derives from the classic result [Sneddon, 1957] that the Pfaffian differential equation $\bm{N}\cdot\mathrm{d}\bm{r}=0$ (solutions of which are neutral trajectories) is integrable if and only if $H=\bm{N}\cdot\nabla\times\bm{N}=0$. However, this classic result holds only in an infinite domain, not in a domain like the real ocean that is bounded and replete with holes such as islands.

Islands and other such holes are important because neutral helices can exist around them even when $H=0$ everywhere in the ocean. A quick way to see this is to imagine the ocean with $H=0$ everywhere except for some region, then build an artificial island over that region (taking great care to not otherwise disturb the ocean state): helical neutral trajectories that existed before the island construction still exist after it, but now $H=0$ everywhere.

In this case, neutral helices exist but only around islands, so cannot be arbitrarily resized (as in Section 2.2). Each helix has a definite, non-infinitesimal, pitch. In such an $H=0$ ocean, neutral surfaces are well-defined surfaces but not unique-depth surfaces. They resemble a multistorey car park with an interior ramp, as schematised in Fig. 1. With many holes in the neutral surface, there may be many interior ramps. Note that holes in a neutral surface are created not just by islands, but also possibly where the surface outcrops or grounds, even in a flat-bottomed ocean, such as in the bottom-left of Fig. 1. However, a hole in the surface does not necessarily produce a neutral helix, a point we shall return to.

2.4 Gradients on the sphere

A neutral surface is supposed to be a surface that is everywhere tangent to the neutral tangent plane. To mathematise this, the 3D gradient of a tracer $C$ is decomposed into components parallel and normal to the neutral tangent plane, as such: $\nabla C=\nabla_{N}C+(\nabla C\cdot\hat{\bm{n}})\hat{\bm{n}}$ where $\hat{\bm{n}}=\bm{N}/|\bm{N}|$ [McDougall et al., 2014]. Thus, displacements $\mathrm{d}\bm{r}$ in the neutral tangent plane satisfy

[TABLE]

having used $\mathrm{d}\bm{r}\cdot\hat{\bm{n}}=0$ and (2). Since $(\rho_{S}\nabla_{N}S+\rho_{\theta}\nabla_{N}\theta)$ is in the neutral tangent plane and $\mathrm{d}\bm{r}$ is arbitrary in that plane, then $\hat{\bm{n}}$ must be such that

[TABLE]

Whereas $\nabla_{N}C$ is a 3D vector field, it is useful to instead work with a 2D vector field determined only from quantities in the surface. To this end, we use the “projected non-orthogonal gradient” introduced by Starr [1945]. This is used throughout studies of neutral surfaces and commonly denoted $\nabla_{n}\equiv\frac{\partial}{\partial{x}}|_{n}\bm{\hat{i}}+\frac{\partial}{\partial{y}}|_{n}\bm{\hat{j}}$. Applying this gradient to a 3D scalar field $C$ produces a 2D vector field, i.e. $\nabla_{n}C\cdot\bm{\hat{k}}=0$. The partial derivatives in $\nabla_{n}$ sample $C$ from the surface in question, but distances are measured only by their horizontal contribution (i.e. they are projected onto a constant height).

This paper uses $\nabla\utilde{C}$ as an alternative notation for $\nabla_{n}C$. The under-tilde is an operator that restricts a 3D field to the surface in question, then projects it onto a sphere (sharing Earth’s centre and radius). The result of this projection, denoted $\mathbb{M}$, is a subset of the sphere, with holes where the Earth has islands and continents, and where the original surface had grounded or outcropped. For example, $\utilde{p}:\mathbb{M}\rightarrow\mathbb{R}$ is a scalar field (with physical units implied) on $\mathbb{M}$; using a geographic coordinate system, $\utilde{p}(\lambda,\phi)$ is the pressure at longitude $\lambda$ and latitude $\phi$ on the original surface in question. Hence, $\utilde{p}$ will be loosely referred to as the pressure on the surface in question, even though $\utilde{p}$ lives on $\mathbb{M}$. The gradient $\nabla\utilde{p}$ is a standard gradient, calculable in spherical coordinates: distance is measured without regard to the radial (vertical) variations of the original surface in question. When the surface in question is a neutral surface, $\nabla\utilde{C}$ and $\nabla_{n}C$ differ only in that the former specifies a single surface at a time, whereas the latter is a 2D vector field living in 3D space.

McDougall et al. [2014] showed that (7) is equivalently expressed using the projected non-orthogonal gradient (simply change “N” to “n”). With the under-tilde notation, this gives

[TABLE]

A truly neutral surface must satisfy (8) exactly (though this is impossible in the real ocean with $H\neq 0$). Indeed, (8) is identical to the first of two equivalent definitions of neutral surfaces given by McDougall [1987a], the other being identical to

[TABLE]

which is derived like (8) but using (1) to express $\bm{N}=\nabla\rho-\rho_{p}\nabla p$. Defining neutral surfaces by (9), rather than (8), is preferable in this work because pressure is monotonic with depth.

3 Neutral surface topology

In this section, we study the topology of properties on a well-defined neutral surface in a hypothetical ocean in which the neutral helicity is everywhere zero.666 If the surface is not unique-depth, $\mathbb{M}$ can have multiple “layers”, but the following theory works equally well. Alternatively, one could work on the neutral surface itself rather than projecting it to the sphere; then $\mathbb{M}$ is a more general Riemannian manifold, and gradients reformulated in terms of the exterior derivative. This is not pursued here, for pedagogical and practical reasons. In this ideal setting, the exactness of (9) has global topological implications. This insight will be used to form topobaric surfaces (Section 5), which are approximately neutral surfaces in the real ocean with non-zero neutral helicity.

First, we must distinguish between contours and level sets. A level set is the disjoint union of any number of contours. They are defined mathematically as follows.

A path from $\bm{x}\in\mathbb{M}$ to $\bm{y}\in\mathbb{M}$ is a continuous function $\mathcal{P}:[0,1]\rightarrow\mathbb{M}$ having $\mathcal{P}(0)=\bm{x}$, and $\mathcal{P}(1)=\bm{y}$.

The contour of $\utilde{p}$ through a point $\bm{x}\in\mathbb{M}$, denoted $\utilde{p}^{-1}(\bm{x})$, is the set of all $\bm{y}\in\mathbb{M}$ for which there exists a path $\mathcal{P}$ from $\bm{x}$ to $\bm{y}$ having $\utilde{p}(\mathcal{P}(t))=\utilde{p}(\bm{x})$ for all $t\in[0,1]$.

The level set of $\utilde{p}$ at $p^{\prime}\in\mathbb{R}$ is the set of all $\bm{x}\in\mathbb{M}$ such that $\utilde{p}(\bm{x})=p^{\prime}$. Mathematically, it is $\utilde{p}^{-1}(p^{\prime})=\{\bm{x}\in\mathbb{M}:\utilde{p}(\bm{x})=p^{\prime}\}$. The level set of a 3D field is often called an isosurface.

3.1 Single-valued functional relations

To get started, first consider a region in $\mathbb{M}$ where $\nabla\utilde{p}\neq\bm{0}$, and where there is precisely one contour of $\utilde{p}$ for any given pressure value. For example, consider just the light blue region in Fig. 2, surrounding point a and with $\utilde{p}<3$. Since $\nabla\utilde{p}$ is orthogonal to a contour of constant $\utilde{p}$, and similarly for $\nabla\utilde{\rho}$ and $\utilde{\rho}$, (9) implies that the contour of $\utilde{p}$ through any point is parallel to the contour of $\utilde{\rho}$ through the same point. Since this is true at all points, a contour of constant $\utilde{p}$ must be exactly a contour of constant $\utilde{\rho}$: the two contours are the same set of points. Specifying a value of $\utilde{p}$ specifies a unique (by assumption) contour of $\utilde{p}$, which is identical to a contour of $\utilde{\rho}$, upon which $\utilde{\rho}$ is constant. Thus, there is a single-valued functional relationship between $\utilde{\rho}$ and $\utilde{p}$:

[TABLE]

for some function $\hat{\rho}$.

Now, the gradient of (10) yields

[TABLE]

Together, (9) and (11) require

[TABLE]

Not only does this give us $\mathrm{d}\hat{\rho}/\mathrm{d}p$, it says that $\utilde{\rho_{p}}$ is also a function of $\utilde{p}$. Specifically,

[TABLE]

where $\hat{\pi}(\mathtt{p})=(\mathrm{d}\hat{\rho}/\mathrm{d}p)(\mathtt{p})$. (The type-face $\mathtt{p}$ distinguishes the function argument from the 3D scalar field $p$; this will be more necessary later.) This is also evident by cross-differentiating (9), to get

[TABLE]

Thus $\nabla\utilde{\rho_{p}}$ and $\nabla\utilde{p}$ are parallel, and the preceding logic applies.777For the notational convenience of (14), $\nabla\utilde{\rho_{p}}$ and $\nabla\utilde{p}$ are temporarily embedded in 3D space, with zero component in the $\bm{\hat{k}}$ direction. A more rigourous notation is $0=J(\utilde{\rho_{p}}\,,\,\utilde{p})$, where $J$ is the Jacobian.

The relation between $\hat{\rho}$ and $\hat{\pi}$ can equivalently be expressed by integrating $\hat{\pi}$ to obtain

[TABLE]

for some constant pressure $p_{c}$ and constant density $\rho_{c}=\hat{\rho}(p_{c})$.

3.2 Multivalued functional relations

In general $\nabla\utilde{p}=\bm{0}$ in some places, and in general there are multiple disjoint contours for a given value of $\utilde{p}$. Recognizing this, the logic that led from (9) to (10) now reveals that $\utilde{\rho}$ can take different values on each of these disjoint contours, so $\hat{\rho}$ in (10) is actually a multivalued function of pressure.

For example, now consider the region in Fig. 2 with $\utilde{p}<7$. The value of $\utilde{\rho}$ must be constant on the $\utilde{p}=2$ contour surrounding point a, but can be different than the constant value of $\utilde{\rho}$ on the $\utilde{p}=2$ contour surrounding point b. Similarly, there are three disjoint contours of the level set $\utilde{p}^{-1}(6)$, upon each of which $\utilde{\rho}$ may take a different value.

At this point, the reader who wishes to see the oceanographic relevance of this multivalued relationship may jump to Section 4.1, taking in Fig. 3a,b.

Our task is to determine geographic regions such that any level set of $\utilde{p}$ has no more than one connected component (contour) in each region, and thus only one value of $\utilde{\rho}$. Then, there are single-valued functions that satisfy (10) within each region. These regions, their meeting points, and the ways they nest into a global structure are encoded by the Reeb graph.

3.3 The Reeb graph

The Reeb [1946] graph captures the essential topological information about connectedness of level sets of a real-valued function on a topological space. For our purposes, the function is $\utilde{p}$ and the topological space is $\mathbb{M}$. Each contour of $\utilde{p}$ is contracted to a single point in the Reeb graph of $\utilde{p}$.888 More formally, the Reeb graph of $\utilde{p}:\mathbb{M}\rightarrow\mathbb{R}$ is the quotient space $\mathbb{M}/\!\!\sim$, where the equivalence relation $\sim$ is such that $\bm{x}\sim\bm{y}$ if $\bm{y}\in\mathbb{M}$ is an element of the contour of $\utilde{p}$ through $\bm{x}\in\mathbb{M}$.

The Reeb graph of $\utilde{p}$ in the preceding example is shown on the right of Fig. 2. In the preceding example, $\utilde{p}^{-1}(2)$ has two disjoint components (contours), so the Reeb graph of $\utilde{p}$ has two corresponding points at $\utilde{p}=2$. Now consider the $\utilde{p}=2$ contour in the light blue region surrounding point a. Moving along a path in physical space from a point on this contour to lower pressure traces out the light blue curve in the Reeb graph, until this ends when the path reaches the pressure minima in physical space located at a. Or consider the $\utilde{p}$ contours as $\utilde{p}$ is swept upwards towards $\utilde{p}=3$: the two $\utilde{p}$ contours approach each other, finally merging into a single contour at $\utilde{p}=3$. The $\utilde{p}=3$ contour contains the saddle point A, which joins three different curves in the Reeb graph.

These points in the Reeb graph are most usefully expressed by a collection of $N$ nodes and $A$ arcs—a graph. Each node $n$ represents a geographic location $\bm{x}_{n}\in\mathbb{M}$ that is a critical point of $\utilde{p}$. Leaf nodes (nodes of degree one) represent local maxima or minima of $\utilde{p}$, and internal nodes (nodes with degree two or more) represent saddles of $\utilde{p}$. Denote the critical value of node $n$ by $p_{n}=\utilde{p}(\bm{x}_{n})$. An arc $a$ is incident upon two nodes, denoted $\ell_{a}$ and $h_{a}$ (think “low” and “high”), having $p_{\ell_{a}}<p_{h_{a}}$, if there is a path $\mathcal{P}_{a}$ from $\bm{x}_{\ell_{a}}$ to $\bm{x}_{h_{a}}$ that is strictly increasing in $\utilde{p}$ [i.e. $s<t$ implies $\utilde{p}(\mathcal{P}_{a}(s))<\utilde{p}(\mathcal{P}_{a}(t))$] and $\mathcal{P}_{a}(t)$ is not on the contour of any critical point of $\utilde{p}$ for all $t\in(0,1)$. Nothing topologically important happens on these paths between pairs of critical points. The region swept out by the contours intersecting such a path is the (associated) region $\mathcal{R}_{a}\subseteq\mathbb{M}$ to arc $a$. Mathematically, $\mathcal{R}_{a}=\cup_{t\in[0,1]}\ \utilde{p}^{-1}(\mathcal{P}_{a}(t))$. In Fig. 2, these are the coloured regions in physical space corresponding to the coloured arcs in the Reeb graph. Now consider points C and D in Fig. 2. There are infinitely many paths through physical space between these points, but all paths that go left of the island have the same associated region. Similarly all paths that go right have the same associated region. Thus the Reeb graph has two arcs between nodes C and D, which together form a cycle. In general, the Reeb graph has as many arcs between two given nodes as there are paths ($\mathcal{P}_{a}$ as above) with distinct associated regions. Such paths are non-homotopic, meaning they cannot be continuously deformed into one another while remaining in $\mathbb{M}$.

As an aside, do not look too closely at the complex boundaries around the islands in Fig. 2. There ought to be many extrema of $\utilde{p}$ lurking around these boundaries, but for illustrative purposes these are ignored. Complex boundaries do add considerable complexity to the Reeb graph for real oceanic data.

The multivalued functions $\hat{\rho}$ and $\hat{\pi}$ become single-valued when the domain is restricted to an associated region. For each arc $a$ in the Reeb graph, there is a single-valued function $\hat{\rho}_{a}:[p_{\ell_{a}},p_{h_{a}}]\rightarrow\mathbb{R}$ such that

[TABLE]

The functions $\hat{\rho}_{a}$ are called branches of the multivalued function $\hat{\rho}$. In fact, for a perfectly neutral surface, $\hat{\rho}_{a}$ is defined by the data $\{(\utilde{p}(\bm{x}),\utilde{\rho}(\bm{x})):\bm{x}\in\mathcal{R}_{a}\}$. The branches of $\hat{\pi}$ are defined similarly, but extra care is needed at saddle points, discussed next.

To see the Reeb graph and its associated regions on oceanographic data, the reader may jump to Section 4.2, taking in Fig. 3c,d.

3.4 Pressure saddles

How do the branches of the multivalued functions $\hat{\rho}$ relate at the saddle pressures, and similarly for $\hat{\pi}$?

For $\hat{\rho}$, the logic leading to (10) applies perfectly well at saddle points. Contours of $\utilde{p}$ are contours of $\utilde{\rho}$, so $\utilde{\rho}$ is constant along contours of $\utilde{p}$, including those through saddle points. Thus, the branches of $\hat{\rho}$ must match continuously at saddle points. That is, for every internal node $s$, \utilde{\rho}(\bm{x}_{s})=\hat{\rho}_{a}\big{(}\utilde{p}(\bm{x}_{s})\big{)} for all arcs $a$ incident upon node $s$. Actually, critical points are contained in the sets $\mathcal{R}_{a}$, so this condition is already covered by (16).

For $\hat{\pi}$, however, the logic leading to (12) does not apply at saddle points. Combining (9) and (10) actually gives $(\mathrm{d}\hat{\rho}/\mathrm{d}p)(\utilde{p})\,\nabla\utilde{p}=\utilde{\rho_{p}}\nabla\utilde{p}$, trivially satisfied where $\nabla\utilde{p}=\bm{0}$. Is it possible that the branches of $\hat{\pi}$ do not match continuously at the saddle pressures? The logic leading to (12) does apply everywhere on the $\utilde{p}$ contour through the saddle, except at the saddle itself. If we remove the saddle point from this contour, then it has multiple connected components, upon each of which $\utilde{\rho_{p}}$ is constant. To rephrase the previous question, could $\utilde{\rho_{p}}$ take different values on these different components? If the answer were yes, then $\utilde{\rho_{p}}$ would jump discontinuously at this point along a trajectory of constant pressure, which implies a discontinuity of salinity or temperature at this point. Assuming the ocean hydrography is continuous, the answer is no: the $\hat{\pi}$ branches must match continuously at the saddle points.

However, rather than a saddle point, now imagine a saddle region within which $\utilde{p}=p_{s}$, a constant. As $\nabla\utilde{p}=\bm{0}$ here, (9) requires $\nabla\utilde{\rho}=\bm{0}$ too, so $\utilde{\rho}$ is constant within this region, and the functional relationship $\hat{\rho}$ holds. However, (14) is trivially satisfied here, even for non-zero $\nabla\utilde{\rho_{p}}$, so $\utilde{\rho_{p}}$ may be non-constant in this region. This entire “flat” region is a node in the Reeb graph. The branches of $\hat{\rho}$ for the arcs incident upon this node must agree at $p_{s}$, but those branches of $\hat{\pi}$ may be undefined at $p_{s}$, and will disagree in the limit as the pressure approaches $p_{s}$. So, the multivalued function $\hat{\rho}$ must be continuous, but $\hat{\pi}$ may have jump discontinuities at saddle regions.

Ultimately, $\utilde{\rho}$ and $\utilde{\rho_{p}}$ are determined by $\utilde{S}$, $\utilde{\theta}$, and $\utilde{p}$, so it is helpful to discuss the structure of salinity and potential temperature on neutral surfaces. Expanding $\nabla\utilde{\rho_{p}}$ in (14) by the chain rule, akin to (1), and using (8) to re-write $\nabla\utilde{S}=-\utilde{\rho_{\theta}}\utilde{\rho_{S}}^{-1}\nabla\utilde{\theta}$, we find

[TABLE]

The term in parentheses is non-zero, and (17) is of the form (14), not the form (9). So, there is a multivalued function $\hat{\theta}$ for which $\utilde{\theta}=\hat{\theta}(\utilde{p})$ with the same behaviour as $\hat{\pi}$, including the possible jump discontinuities. Of course, one can equally well consider salinity, and find a multivalued function $\hat{S}$ for which $\utilde{S}=\hat{S}(\utilde{p})$, again with the same behaviour as $\hat{\pi}$. This is how $\utilde{\rho_{p}}$ can vary inside a flat pressure region: in such a region, $\utilde{S}$ and $\utilde{\theta}$ can vary in compensatory ways to maintain $\utilde{\rho}$ constant, but cannot simultaneously maintain $\utilde{\rho_{p}}$ constant. Rest assured that, while $\hat{\pi}$, $\hat{\theta}$, and $\hat{S}$ may be discontinuous multivalued functions of pressure, the fields $\utilde{\rho_{p}}$, $\utilde{\theta}$, and $\utilde{S}$ are themselves continuous in space.

As a physical field, not a mathematical construction, we might expect there are no extended regions where $\utilde{p}$ is truly constant. If so, the branches of $\hat{\pi}$ meet continuously and are given by $\hat{\pi}_{a}:[p_{\ell_{a}},p_{h_{a}}]\rightarrow\mathbb{R}$ such that

[TABLE]

analogous to (16). Where a saddle $\bm{x}_{s}$ has constant pressure over an extended region, for each arc $a$ incident upon node $s$, simply exclude $p_{s}$ from the domain of $\hat{\pi}_{a}$ and restrict $\bm{x}$ to the interior of $\mathcal{R}_{a}$.

A final, technical point is that the strict definition of the Reeb graph is rooted in Morse theory, which requires all critical points (of $\utilde{p}$) to be non-degenerate. This and other properties of Morse functions ensure all nodes of the Reeb graph are degree one or three. Of course, a region of constant $\utilde{p}$ is full of degenerate critical points. Moreover, piecewise-linear functions constructed from numerical data often fail to be Morse. Nonetheless the Reeb graph still exists and can be computed, though it may have nodes with degree 2 or 4 or more [Cole-McLaughlin et al., 2003, Doraiswamy and Natarajan, 2013].

3.5 Islands and holes and cycles

Islands and other holes in the neutral surface can impose an additional constraint on $\hat{\pi}$. Consider a neutral trajectory, a path $\mathcal{P}$ in the neutral surface. The change of in-situ density from (a reference point) $\bm{x}_{0}\in\mathbb{M}$ to (any point) $\bm{x}\in\mathbb{M}$ is

[TABLE]

using (9) for the second equality. Having assumed the neutral surface is well-defined, these integrals are path-independent.

The path $\mathcal{P}$ in (19) corresponds to a walk in the Reeb graph of $\utilde{p}$, alternately passing along arcs and through nodes in the order $a_{0}$, $n_{1}$, $a_{1}$, …, $n_{\mathcal{J}}$, $a_{\mathcal{J}}$. Specifically, $\bm{x}_{0}\in\mathcal{R}_{a_{0}}$ and $\bm{x}\in\mathcal{R}_{a_{\mathcal{J}}}$. From (16), we also have $\utilde{\rho}(\bm{x}_{0})=\hat{\rho}_{a_{0}}(p_{0})$ and $\utilde{\rho}(\bm{x})=\hat{\rho}_{a_{\mathcal{J}}}(p)$, where $p_{0}=\utilde{p}(\bm{x}_{0})$ and $\mathtt{p}=\utilde{p}(\bm{x})$. Thus the path integral in (19) becomes a “graph integral”,

[TABLE]

having also used (18). This generalizes (15) when the functions $\hat{\rho}$ and $\hat{\pi}$ are multivalued.

The only way path-dependence can affect (3.5) is via cycles in the Reeb graph. To see this, suppose the Reeb graph has no cycles (a tree), and consider any two nodes. There is a unique walk between these nodes having no repeated nodes, called a straight walk. The result of (3.5) for any walk between these nodes is identical to that for the straight walk, so (19) is path-independent.

If there are holes in $\mathbb{M}$ (such as created by islands), there may be cycles in the Reeb graph. In this case, there may be multiple straight walks between two nodes, corresponding to the path $\mathcal{P}$ navigating one way or another around a hole. Consider a cycle whose walk is $n_{1}$, $a_{1}$, $n_{2}$, $...$, $n_{\mathcal{J}-1}$, $a_{\mathcal{J}-1}$, $n_{\mathcal{J}}=n_{1}$. Without loss of generality, suppose $\bm{x}_{0}=\bm{x}=\bm{x}_{n_{1}}$. Then (3.5) becomes

[TABLE]

This constraint must be satisfied by $\hat{\pi}$ for each cycle in the Reeb graph of $\utilde{p}$.

It may come as a surprise that the existence of holes in $\mathbb{M}$ does not immediately guarantee that the Reeb graph of $\utilde{p}$ will have cycles. A cycle can only occur if there is a $\utilde{p}$ contour that intersects the hole at precisely one end; the contour cannot close on itself, so its other end must intersect a different hole. If only one hole is contained within a $\utilde{p}$ contour (such as the top island in Fig. 2) or between two $\utilde{p}$ contours, that hole does not produce a cycle: any contour intersecting the hole cannot cross the bounding contour(s), and so must intersect the hole twice. As a corollary to this, there are no cycles when $\mathbb{M}$ has a single hole. Note that a bounded ocean containing a single island, as in Fig. 1, is topologically the same as an aqua-planet with two islands.

It is conceivable that the pitch of neutral helices in the open ocean is small ($H\approx 0$ globally), yet around islands or other holes this pitch may be large. This seems plausible if interior ocean dynamics naturally tend to destroy (bring to zero) neutral helicity, as McDougall and Jackett [2007] tentatively suggested. A preliminary analysis of this possibility is given in B.

3.6 Summary

When $H=0$ everywhere, neutral surfaces are well-defined surfaces. On a well-defined neutral surface, the in-situ density, $\utilde{\rho}$, is a multivalued function of the pressure, $\utilde{p}$. The neutral surface can be covered by regions within each of which $\utilde{\rho}$ is a single-valued function of $\utilde{p}$, called a branch of the multivalued function. Each of these regions, and so too each of the branches, is associated with an arc of the Reeb graph of $\utilde{p}$. The matching conditions of these branches are determined by the structure of the Reeb graph. Each internal node of the graph is associated with a saddle of $\utilde{p}$. All branches associated with arcs incident to a common node must match continuously at the pressure associated with that node.

Moreover, the partial derivative of in-situ density with respect to pressure on a well-defined neutral surface, $\utilde{\rho_{p}}$, is also a multivalued function of $\utilde{p}$. It is single-valued within the same regions as above and these branches must integrate to zero around every cycle in the Reeb graph.

Following an oceanographic example in Section 4, this theory will be used in Section 5 to develop unique-depth, approximately neutral surfaces in the real ocean with $H\neq 0$.

4 Illustrative example

4.1 Pressure and in-situ density on an approximately neutral surface

The main computations and tests presented will use high resolution data (spatially and temporally), but for illustrative purposes smoother fields are desired: the OCCA 2004—2006 climatology [Forget, 2010] provides these. The potential density referenced to $1000\text{\,}\mathrm{d}\mathrm{b}\mathrm{a}\mathrm{r}$ is (crudely) calculated using the climatological salinity and potential temperature, then the isopycnal surface intersecting ($180^{\circ}\mathrm{E}$, $0^{\circ}\mathrm{N}$, $1000\text{\,}\mathrm{d}\mathrm{b}\mathrm{a}\mathrm{r}$) is found. Then, the pressure on this isopycnal, $\utilde{p}$, is slightly adjusted to globally minimize the error from neutrality, resulting in an $\omega$-surface [Klocker et al., 2009]. This is the “illustrative surface”. This calculation is crude, essentially treating climatological data as if it were instantaneous, but suffices for illustration. It also ignores the fact that OCCA is a Boussinesq model (the implications of which will be discussed in Section 5.3).

The pressure on the illustrative surface, $\utilde{p}$, is mapped in Fig. 3a. The (thick white) $\utilde{p}=$500\text{,}\mathrm{d}\mathrm{b}\mathrm{a}\mathrm{r} level set possesses two clearly disjoint contours in the Southern Ocean and North Atlantic. If this surface were truly neutral, the *in-situ* density would be constant along each of these pressure contours, yet possibly different between them. This phenomenon gives rise to the multivalued nature of $\hat{\rho}$ evident from the scatter plot of $\utilde{\rho}$ vs. $\utilde{p}$ shown in Fig. [3](#S4.F3)b. (A reference profile $R(S_{0},\theta_{0},\utilde{p})$ is subtracted from $\utilde{\rho}$ purely for illustrative purposes, as $\utilde{\rho}$ vs. $\utilde{p}$ looks essentially linear. The values $S_{0}$ and $\theta_{0}$ are taken as $\utilde{S}$ and $\utilde{\theta}$ at ($180^{\circ}\mathrm{E}$, $0^{\circ}\mathrm{N}$).) Indeed, at $\utilde{p}=$500\text{\,}\mathrm{d}\mathrm{b}\mathrm{a}\mathrm{r}, $\utilde{\rho}$ is very nearly one of two values, each corresponding to one of the two disjoint contours of $\utilde{p}=$500\text{,}\mathrm{d}\mathrm{b}\mathrm{a}\mathrm{r}$$. Actually, there is some scatter of $\utilde{\rho}$ around these values. A potential density surface would show more scatter, whereas a truly neutral surface would show no scatter whatsoever—but the essential multivalued nature would remain.

4.2 The Reeb graph and its associated regions

The Reeb graph of $\utilde{p}$ for the illustrative surface contains 1,370 arcs—too many to be particularly informative. Most of these represent very small regions in physical space, often hugging coasts where local extrema of $\utilde{p}$ are common. Even with climatological fields, the Reeb graph requires simplification. This is a complicated task, and used only for this illustration, so only a brief description is given here; see Stanley [2018] for further details. In brief, small holes are filled in with extremely large values (so contours go around them), then the Reeb graph is calculated, then those filled holes removed from the associated regions. The simplification method that Carr et al. [2010] used on trees is then used on this graph with cycles. It iteratively removes the least important leaf node until only a specified number of arcs remain; it never destroys cycles and never produces a node with only arcs leading to lower pressures, or only arcs leading to higher pressures.

The Reeb graph of $\utilde{p}$ on the illustrative surface is computed and simplified down to 43 arcs. The graph itself is drawn in Fig. 3c with nodes positioned according to the pressure (ordinate) and longitude (abscissa) of their associated critical points.999This graph drawing conveys additional useful information at the cost of some arcs overlapping. Deciding the best placement of nodes and arcs of a Reeb graph is non-trivial [Heine et al., 2011]. The regions $\mathcal{R}_{a}$ associated with each arc $a$ of this simplified Reeb graph are mapped in Fig. 3d. Three arcs cross the (thick white) $500\text{\,}\mathrm{d}\mathrm{b}\mathrm{a}\mathrm{r}$ level. Each of these supports a distinct in-situ density, as seen in Fig. 3b. (At OCCA’s coarse-resolution, Mediterranean outflow on this surface jumps from $280\text{\,}\mathrm{d}\mathrm{b}\mathrm{a}\mathrm{r}$ to $850\text{\,}\mathrm{d}\mathrm{b}\mathrm{a}\mathrm{r}$ between neighbouring grid cells, so the $500\text{\,}\mathrm{d}\mathrm{b}\mathrm{a}\mathrm{r}$ contour in Fig. 3a is hardly visible near Gibraltar, though it exists.) The associated regions for these arcs lie in the Southern Ocean (red), the North Atlantic (blue), and the Mediterranean (maroon). The point on the Reeb graph on the blue arc at $500\text{\,}\mathrm{d}\mathrm{b}\mathrm{a}\mathrm{r}$ represents the entire, but single, $500\text{\,}\mathrm{d}\mathrm{b}\mathrm{a}\mathrm{r}$ contour of $\utilde{p}$ in North Atlantic region. Following this arc upwards corresponds to a path, on the surface, in which the pressure decreases monotonically. It can be followed until $90\text{\,}\mathrm{d}\mathrm{b}\mathrm{a}\mathrm{r}$, where node A is reached, corresponding to a saddle of $\utilde{p}$ at the Grand Banks of Newfoundland (indicated by A in Fig. 3d). All such pressure-monotonic paths in the North Atlantic between $90\text{\,}\mathrm{d}\mathrm{b}\mathrm{a}\mathrm{r}$ (node A) and $865\text{\,}\mathrm{d}\mathrm{b}\mathrm{a}\mathrm{r}$ (node C) are equivalent in the Reeb graph. As pressure increases along these paths, the in-situ density increases identically between all such paths (were this surface truly neutral). At the saddle A, there are two options: the path can be continued monotonically to shallower pressure, ultimately reaching the local minima of $\utilde{p}$ in the North Atlantic (node c); or, the path can descend to higher pressure, ultimately reaching the local maxima of $\utilde{p}$ in Baffin Bay (node e).

Similarly, the point on the Reeb graph on the dark blue arc at $1350\text{\,}\mathrm{d}\mathrm{b}\mathrm{a}\mathrm{r}$ represents the entire, but single, $1350\text{\,}\mathrm{d}\mathrm{b}\mathrm{a}\mathrm{r}$ contour in the South Pacific. As one moves along this arc to node x, the contour contracts to a point, the maxima of $\utilde{p}$ in the South Pacific. Or, moving to shallower pressures, the contour grows, until it intersects another contour of the same level set from the South Indian, at the saddle node W. Note, at $1250\text{\,}\mathrm{d}\mathrm{b}\mathrm{a}\mathrm{r}$ in the South Pacific, the in-situ density seems to branch into two distinct values. This is because a small closed $\utilde{p}$ contour, west of New Zealand, has been added to the dark blue region by the simplification process (compare Fig. 3a and d). Much of the scatter in Fig. 3b is actually due to such merging of regions by the simplification process, rather than poor neutrality of the $\omega$-surface.

Now consider an upward-monotonic path from the South Pacific $\utilde{p}$ maxima (node x) to the Amundsen Sea (node b). The most direct route is equivalent, in the Reeb graph, to a journey around the South Pacific visiting nodes W, V, S, G, B, then b. Alternatively, one can take a South Indian route, through nodes W, V, U, T, Q, G, B, then b. These are two ways around New Zealand. As New Zealand (united, at this depth) represents a hole in $\mathbb{M}$, it creates a cycle in the Reeb graph: V, S, G, Q, T, U, V. Madagascar creates a second cycle in the Reeb graph, U, T, R, U. A third, and final, cycle is the loop around both New Zealand and Madagascar. There are many more holes in $\mathbb{M}$, but only these cycles. Antarctica does not produce a cycle in the Reeb graph, because there exists a contour of $\utilde{p}$ encircling Antarctica, as discussed in Section 3.5. (One can, equivalently, think of the non-Antarctic continents as islands in an ocean bounded by Antarctica.) The simplification procedure pre-emptively removed other, smaller holes before calculating the Reeb graph. No such simplification will be used when computing topobaric surfaces, discussed next.

5 Topobaric surfaces

Studying the topology of a neutral surface in an ocean with zero neutral helicity (Section 3) revealed the existence of a multivalued function $\hat{\rho}$ whose branches satisfy (16). The real ocean has neutral helicity that is non-zero yet small, so surfaces can only be approximately neutral. On such surfaces, $\nabla\utilde{\rho}-\utilde{\rho_{p}}\nabla\utilde{p}$ is non-zero, but small enough that the multivalued function $\hat{\rho}$ still usefully describes the approximately neutral surface, albeit a small error must be added to (16); this error is the scatter in Fig. 3b.

Topobaric surfaces turn this around, by forcing (16) to be exact for a given $\hat{\rho}$. To obtain nearly neutral surfaces, $\hat{\rho}$ must be chosen carefully. Clearly, choosing a constant $\rho_{c}=\hat{\rho}(p)$ is a bad choice. The theory for how to choose $\hat{\rho}$ is given next (Section 5.1), followed by an algorithm to construct topobaric surfaces (Section 5.2). This method is applied to data (described in Section 5.3) to numerically compute topobaric surfaces and compare them with other approximately neutral surfaces (Section 5.4).

5.1 Theory

Because the neutral helicity is non-zero, and because we wish to make unique-depth surfaces in the presence of multiple islands and other holes in $\mathbb{M}$, we cannot construct perfectly neutral surfaces. Starting in the simple setting of Section 3.1 where $\hat{\rho}$ and $\hat{\pi}$ are single-valued, this means that a unique-depth surface cannot simultaneously satisfy the $\hat{\rho}$ relation (10) and the $\hat{\pi}$ relation (13) exactly. We could choose $\hat{\rho}$ upfront then force (10) to hold exactly and not worry about (13). But, for neutrality, what really matters is not the values of $\hat{\rho}$ but its derivative, as (12) shows. So, we choose to satisfy (10) exactly and approximate (13) with an empirically fit $\hat{\pi}$, which is integrated to determine $\hat{\rho}$ according to (15). To see that this maximizes neutrality, combine (10) and (15) and take the gradient via the Leibniz integral rule to find $\nabla\utilde{\rho}=\hat{\pi}(\utilde{p})\nabla\utilde{p}$. The approximate version of (13) then yields $\hat{\pi}(\utilde{p})\nabla\utilde{p}\approx\utilde{\rho_{p}}\nabla\utilde{p}$. Combining these gives an approximate version of the neutrality condition (9). The better $\hat{\pi}(\utilde{p})$ approximates $\utilde{\rho_{p}}$, the closer the resulting surface will be to neutral.101010Another possibility would be to make (13) exact and obtain $\hat{\pi}$ by differentiating an empirically fit $\hat{\rho}$ that approximates (10). But, as above when choosing $\hat{\rho}$ upfront, there is no reason to believe the surface will be nearly neutral. The condition (14) would be guaranteed, but this does not imply neutrality (9).

Translating this into the case of multivalued functions, a topobaric surface, denoted a $\tau$-surface, satisfies the $\hat{\rho}$ relation (16) exactly, where $\hat{\rho}$ is obtained by integrating $\hat{\pi}$ according to (3.5) with a given a reference location $\bm{x}_{0}$, and $\hat{\pi}$ satisfies the cycle constraint (21). To make topobaric surfaces as neutral as possible, $\hat{\pi}$ is chosen to approximate (18) as best as possible. The cycle constraint (21) ensures topobaric surfaces are unique-depth surfaces. Also, the exactness of (16) ensures topobaric surfaces possess an exact geostrophic streamfunction [Stanley, 2019].

Topobaric surfaces allow $\hat{\pi}$ to be discontinuous at all pressure saddles. There are three justifications for this. First, constant pressure regions could exist in the continuous $\utilde{p}$, but not be present in a discrete data representation of $\utilde{p}$. Second, with non-zero neutral helicity (14) becomes $\epsilon=\bm{\hat{k}}\cdot\nabla\utilde{\rho_{p}}\times\nabla\utilde{p}$ for some small scalar field $\epsilon$. Near saddle pressures $|\nabla\utilde{p}|$ is small, so $|\nabla\utilde{\rho_{p}}|$ can be large while maintaining a small $\epsilon$, and thus $\utilde{\rho_{p}}$ can change rapidly near saddle pressures. This differs from $\bm{\epsilon}=\nabla\utilde{\rho}-\utilde{\rho_{p}}\nabla\utilde{p}$, which requires $|\nabla\utilde{\rho}|$ to be small near pressure saddles. Third, numerical tests show that requiring $\hat{\pi}$ to match continuously at the pressure saddles produces surfaces that are less neutral.

This definition of topobaric surfaces is circular: $\utilde{p}$ must be known in order to calculate the Reeb graph of $\utilde{p}$, which is used in defining the multivalued functions $\hat{\pi}$ and $\hat{\rho}$, the latter being an implicit definition for $\utilde{p}$ via (16).

5.2 Methods

To overcome the preceding circular definition, surfaces are built by an iterative algorithm that converges to a topobaric surface, as follows.

Begin with an approximately neutral surface, with pressure $\utilde{p}$, and a reference location $\bm{x}_{0}$. 2. 2.

Compute the Reeb graph of $\utilde{p}$. 3. 3.

Empirically fit the branches $\hat{\pi}_{a}$ using the data $\{(\utilde{p}(\bm{x}),\utilde{\rho_{p}}(\bm{x})):\bm{x}\in\mathcal{R}_{a}\}$ for each arc $a$ of the Reeb graph, subject to the cycle constraints (21). 4. 4.

Obtain $\hat{\rho}$ by integrating $\hat{\pi}$ according to (3.5). 5. 5.

Update $\utilde{p}$ with that satisfying (16), which is a root-finding problem for each water column. 6. 6.

Return to Step 2, unless a convergence test is passed.

Surfaces constructed by a finite number of the above iterations are loosely referred to as topobaric surfaces. The remainder of this section describes each of these steps, in turn.

5.2.1 Initial surface and reference location

The initial surface can be any approximately neutral surface, chosen by the user. The resulting topobaric surface is somewhat dependent on this choice; starting from an $\omega$-surface yields slightly better results than starting from a potential density surface (not shown).

Results shown in this paper all use topobaric surfaces initialized from potential density surfaces, to ensure the method succeeds when initialized from a surface that is not particularly neutral, globally. Also, one advantage of topobaric surfaces is their computational speed, which is defeated if one must first construct an $\omega$-surface.

Next, the user selects a reference location $\bm{x}_{0}$. Then, by linear interpolation in the water column $\bm{x}_{0}$, record $p_{0}=\utilde{p}(\bm{x}_{0})$, $S_{0}=\utilde{S}(\bm{x}_{0})$, and $\theta_{0}=\utilde{\theta}(\bm{x}_{0})$; also record $\rho_{0}=R(S_{0},\theta_{0},p_{0})=\utilde{\rho}(\bm{x}_{0})$. The iteratively updated surface will maintain these properties at $\bm{x}_{0}$. (The initial potential density surface can have any reference pressure, but $p_{0}$ is a good choice.)

5.2.2 Computing the Reeb graph

Fast, robust, and general computation of Reeb graphs from discretised data has only been achieved recently. Carr et al. [2003] developed a fast and general method to compute the contour tree (a Reeb graph with no cycles) in any dimension, based on sorting the data and sweeping through it twice. This is used by Doraiswamy and Natarajan [2013] to compute the Reeb graph, by decomposing the space into a collection of loop-free regions. Their software, called ReCon, is written in Java. It is slightly modified to work with 64-bit floating point numbers, and to communicate directly with MATLAB.

ReCon requires its input to be a simplical mesh. This is a collection of vertices and, in 2D, a collection of triangles. With function values specified on the vertices, a piecewise-linear function may be constructed by linear interpolation. This differs from data on a rectilinear grid, as is common for oceanographic data. Rectilinear data may be bilinearly interpolated, but this non-linear interpolation can introduce new critical points, thereby changing the topology of level sets of $\utilde{p}$. How one constructs the simplical mesh from rectilinear data can matter [Carr et al., 2006]. One method is to add a new vertex at the centre of each rectangle by four-way averaging, then split each rectangle into four triangles. This is how contours are typically computed from rectilinear data by the marching cubes algorithm. However, the associated regions for some arcs of the Reeb graph can contain only these extra vertices. To then fit $\utilde{\rho_{p}}$ to $\utilde{p}$ would require calculating $\utilde{\rho_{p}}$ on these extra vertices; this requires averaging $\utilde{S}$ and $\utilde{\theta}$ onto these extra vertices, but this is undesirable because $\utilde{\rho_{p}}$ is a non-linear function of $\utilde{S}$, $\utilde{\theta}$, and $\utilde{p}$. We use a simple method, splitting each rectangle into two triangles; where a rectangle has one data point missing, only one triangle is produced. Using global oceanographic data, both methods produce Reeb graphs that appear similar.

For simplicity, one connected component of the approximately neutral surface is handled at a time. A graph (not the Reeb graph) is constructed from the simplical mesh, having a node at every vertex of the mesh, and an arc between two nodes when their corresponding vertices share a face of the mesh. The connected component of this graph containing the point $\bm{x}_{0}$ is found [Tarjan, 1972].111111Faster, image manipulation methods are not sufficient because the simplical decomposition can remove the odd grid point, such as those with ground for 7 of 8 neighbours. This guarantees that the Reeb graph is a connected graph.

Finally, all critical points must have unique values. Following standard practice [Doraiswamy and Natarajan, 2013], any duplicate values of $\utilde{p}$ are perturbed by a machine-precision amount until all vertices have unique values.

ReCon computes the Reeb graph in $\mathcal{O}(v\log v+sn)$ time, where $v$ is the number of vertices, $n$ is the number of triangles, and $s$ is the number of saddles of the simplical mesh.

5.2.3 Empirically fit $\hat{\pi}$, with cycle constraints

The goal now is to use the data $\{(\utilde{p}(\bm{x}),\utilde{\rho_{p}}(\bm{x})):\bm{x}\in\mathcal{R}_{a}\}$ to empirically fit a function $\hat{\pi}_{a}$ that approximates (18), for each arc $a$ of the Reeb graph of $\utilde{p}$, subject to the cycle constraint (21) for each cycle.

Rather than finding all cycles, we need only find a cycle basis, which is a minimal set of cycles out of which all other cycles can be produced by “addition”: taking the arcs in one, but not both, of two cycles to produce a third cycle.121212If the original two cycles are disjoint, the third is not a cycle but a more general object, an Eulerian subgraph—see A. If (21) holds for the first two cycles, it will hold for the third. A cycle basis is determined in two steps.

First, find the arc $a_{1}$ having $\bm{x}_{0}\in\mathcal{R}_{a_{1}}$, and let $m_{1}$ and $n_{1}$ be the nodes that $a_{1}$ is incident upon. Perform a breadth-first search, starting with $m_{1}$ as the initially discovered node. This iteratively discovers nodes adjacent to previously discovered nodes (the first step is rigged so that $n_{1}$ is discovered next). The result is a sequence of nodes $n_{1},...,n_{N-1}$ and another sequence of nodes $m_{1},...,m_{N-1}$ such that $m_{j}$ was discovered before its adjacent node $n_{j}$. (The set $\{m_{1},n_{1},..,n_{N-1}\}$ is all nodes in the graph.) The set of arcs $\{a_{1},...,a_{N-1}\}$, where $a_{j}$ is incident upon both $m_{j}$ and $n_{j}$, forms a spanning tree.

Second, for each arc $a$ in the Reeb graph but not in the spanning tree, perform another breadth-first search in the spanning tree starting at a node upon which $a$ is incident, and stopping upon discovery of the other node upon which $a$ is incident. This finds the shortest walk in the spanning tree between the adjacent nodes of $a$. This, together with $a$, gives a cycle. All such cycles form the cycle basis.

Now, the branches $\hat{\pi}_{a}$ can be empirically fit. What form should $\hat{\pi}_{a}$ take? The equation of state $R$ is usually expressed by a rational function of $S$, $\theta$, and $p$, so $\partial_{p}R$ is also a rational function of $S$, $\theta$, and $p$. It might make sense to fit $\hat{\pi}_{a}$ as a rational function of $p$, but a functional form with fewer degrees of freedom is preferable to avoid over-fitting the data. A form with two degrees of freedom will never be under-determined, because the chosen simplical mesh has at least two data points for each arc. Thus, we use the simple form

[TABLE]

where the constants $K_{a}$ and $L_{a}$ are to be determined. The addition of $\partial_{p}R(S_{0},\theta_{0},\utilde{p})$ helps capture some of the non-linear behaviour of $\utilde{\rho_{p}}$ with respect to $\utilde{p}$. For all arcs $a$ that are not in a cycle in the cycle basis, $K_{a}$ and $L_{a}$ are determined by ordinary least squares, fitting $K_{a}+L_{a}\,(\utilde{p}-p_{\ell_{a}})$ to $\utilde{\rho_{p}}-\partial_{p}R(S_{0},\theta_{0},\utilde{p})$ within $\mathcal{R}_{a}$. (That these may each be fit independently is a boon of allowing $\hat{\pi}$ to meet discontinuously at the saddle pressures.) The remaining branches are fit similarly, but as a single, coupled problem that also satisfies the cycle constraints (21) for each cycle; this is done using MATLAB’s lsqlin function. As the pressure difference between adjacent nodes is typically small, affine linear functions perform very well: in practice, the $\partial_{p}R(S_{0},\theta_{0},p)$ term in (22) could be excluded with very little detriment.

5.2.4 Obtain $\hat{\rho}$ by integrating $\hat{\pi}$

Having chosen the branches $\hat{\pi}_{a}$ in (23), they are integrated according to (3.5) to obtain the branches $\hat{\rho}_{a}$. In practice, each branch is first analytically integrated as in (15) to get

[TABLE]

Each branch $a$ has a free constant of integration, $J_{a}$. One of these, $J_{a_{1}}$, is set by requiring $\rho_{0}=\hat{\rho}_{a_{1}}(p_{0})$. The remainder are used to ensure the branches of $\hat{\rho}$ match continuously at the saddle pressures. This is done using the discovery order of the nodes from the previous step, as follows. First, record $\rho_{m_{1}}=\hat{\rho}_{a_{1}}(p_{m_{1}})$. Then, for each $j=2,...,N-1$, determine $J_{a_{j}}$ from $\hat{\rho}_{a_{j}}(p_{m_{j}})=\rho_{m_{j}}$ and record $\rho_{n_{j}}=\hat{\rho}_{a_{j}}(p_{n_{j}})$. Finally, for each arc $a$ in the cycle basis, determine $J_{a}$ from $\hat{\rho}_{a}(p_{\ell_{a}})=\rho_{\ell_{a}}$. The cycle constraint (21) ensures this is identical to determining $J_{a}$ from $\hat{\rho}_{a}(p_{h_{a}})=\rho_{h_{a}}$.

5.2.5 Updating $\utilde{p}$

With all branches of $\hat{\rho}$ chosen, we must now update $\utilde{p}$ to satisfy (16). Specifically, for each arc $a$ and for each geographic position $\bm{x}\in\mathcal{R}_{a}$, we set $\utilde{p}(\bm{x})=p^{\prime}$ where $p^{\prime}$ solves

[TABLE]

where $\mathcal{S}$ and $\vartheta$ are versions of $S$ and $\theta$ with pressure as the vertical coordinate. Mathematically, $\mathcal{S}(\bm{x},p^{\prime})=S(\bm{x},z^{\prime})$ where $z^{\prime}$ solves $p(\bm{x},z^{\prime})=p^{\prime}$, and similarly for $\vartheta$.

Bisection is used to solve (24). Since multiple solutions are possible, an initial guess is provided, based on $\utilde{p}$ from the previous iteration. A small interval around the initial guess is tested for a sign change at its limits. If a sign change is found, bisection proceeds inside this interval. Otherwise, the interval is (geometrically) expanded until a sign change is found and bisection can proceed, or the shallowest and deepest grid cells are reached. If the latter, no solution is found, meaning the updated surface has outcropped or incropped.131313There is currently no capacity for “wetting”, whereby subsequent iterations retest water columns that previously outcropped or incropped. This would require a way to define $\hat{\rho}$ for this water column. Perhaps $\hat{\rho}_{a}$ could be used when the water column is entirely surrounded by a single region $\mathcal{R}_{a}$. At the boundary between regions, perhaps an average of these branches could be used. MATLAB’s code generation is used to turn this into fast C executables (MEX).

5.2.6 Iteration

Recall that $\hat{\pi}$ is allowed to be discontinuous only at the saddle pressures. The iterative method is required because solving for new pressures (step 5) causes the saddle points to change, and thus after an iteration, $\hat{\pi}$ will be discontinuous inside some regions. This causes large errors along the pressure contours that were formerly through pressure saddles. As the whole algorithm iterates, $\utilde{p}$ converges, and discontinuities of $\hat{\pi}$ occur only at the pressure saddles.

The stopping criterion may be chosen by the user. The default is to stop when the root-mean-square change of $\utilde{p}$ is less than $10^{-3}$ dbar. Provided this value is sufficiently small, the results are not sensitive to the choice of this stopping value.

5.3 Data & Boussinesq models

With these methods, topobaric surfaces are constructed and tested using ECCO2 [Menemenlis et al., 2005] data, having $0.25^{\circ}\times 0.25^{\circ}$ horizontal resolution, on 22–24 December 2002. A single archived time-step is chosen to make the task as hard as possible: the Reeb graph of a smoother climatology is considerably simpler. In truth, neutral helices possess a temporal as well as spatial dimension [Klocker and McDougall, 2010], which we are not considering here.

Boussinesq models, such as ECCO2, calculate the in-situ density from a Boussinesq equation of state $B$ that uses depth $z$ (negative and decreasing downward) rather than in-situ pressure. That is,

[TABLE]

where $\rho_{B}$ is the Boussinesq reference density and $g$ the gravitational acceleration [Young, 2010]. The preceding theory is modified to the Boussinesq case simply by swapping $p$ for $z$ and $R$ for $B$. For instance, the in-situ density gradient becomes

[TABLE]

where now $\rho_{S}=\partial_{S}B(S,\theta,z)$, $\rho_{\theta}=\partial_{\theta}B(S,\theta,z)$, and $\rho_{z}=\partial_{z}B(S,\theta,z)<0$. Also, $\nabla z=\bm{\hat{k}}$. The neutral surface relation (9) becomes

[TABLE]

Now $\hat{\rho}$ and $\hat{\pi}$ are functions of $z$, not $p$. However, the essence of the theory is unchanged. (In practice, the topobaric code trivially switches to the Boussinesq case by using $B$ instead of $R$, and internally treating $z$ as positive and increasing downwards, like $p$.)

5.4 Results

To assess topobaric surfaces, they are compared against five other approximately neutral surfaces. Computation time is briefly discussed, but mostly the comparison rests on neutrality.

5.4.1 Six classes of approximately neutral surfaces

Six types of approximately neutral surfaces will be constructed:

potential density surfaces [Wüst, 1935], isosurfaces of $\sigma_{1}=B(S,\theta,$-1000\text{,}\mathrm{m}$)$ or $\sigma_{2}=B(S,\theta,$-2000\text{,}\mathrm{m}$)$; 2. 2.

in-situ density anomaly surfaces [Montgomery, 1937], isosurfaces of $\delta=\rho-B(S_{\delta},\theta_{\delta},z)$ where $S_{\delta}$ and $\theta_{\delta}$ are constants; 3. 3.

neutral density surfaces, isosurfaces of $\gamma^{n}$ [Jackett and McDougall, 1997] 4. 4.

$\sigma_{\nu}$-surfaces, “orthobaric”; 5. 5.

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

$\tau$-surfaces, topobaric.

The “orthobaric” surface is not actually an isosurface of the de Szoeke et al. [2000] 3D orthobaric density; quotation marks around “orthobaric” help indicate this distinction. Rather, it is a special case of a topobaric surface but with the whole surface fit together in a single region, and $\hat{\pi}$ fit as a cubic spline with knots at $z_{max}$, $-200\text{\,}\mathrm{m}$, $-1500\text{\,}\mathrm{m}$, $-1800\text{\,}\mathrm{m}$, and $z_{min}$, where $z_{max}$ and $z_{min}$ are the shallowest and deepest depths found on the surface, respectively. This allows us to test the importance of geography in topobaric surfaces. Similarly, the $\sigma_{1}$-, $\sigma_{2}$-, and $\delta$-surfaces are not constructed as isosurfaces of 3D scalar fields using vertical interpolation, but rather as solutions of a non-linear equation in each water column, found by bisection, e.g. at each $\bm{x}$, solving for $z$ in $\rho(\bm{x},z)-B(S_{\delta},\theta_{\delta},z)-\delta=0$ for some constant (isovalue) $\delta$.

The $\omega$-surface is constructed first; it is initialized from a $\sigma_{2}$-surface intersecting ($180^{\circ}\mathrm{E}$, $0^{\circ}\mathrm{N}$, $-2000\text{\,}\mathrm{m}$), but it heaves during its iterative procedure, finishing at ($180^{\circ}\mathrm{E}$, $0^{\circ}\mathrm{N}$, $-1988.60\text{\,}\mathrm{m}$). The other five surfaces are constructed to intersect this latter point. Specifically, this means (a) $\sigma_{2}=1036.9551~{}\mathrm{kg}~{}\mathrm{m}^{-3}$, (b) $\delta=-3.4759\times 10^{-3}~{}\mathrm{kg}~{}\mathrm{m}^{-3}$ with $S_{\delta}=$34.6568\text{,}\mathrm{p}\mathrm{s}\mathrm{u} and $\theta_{\delta}=2.0899^{\circ}$C taken as mean values on the $\omega$-surface between $55^{\circ}$S and $50^{\circ}$S (in the Southern Ocean, the nexus of the other oceans) , and (c) $\gamma^{n}=27.9248$. Moreover, the “orthobaric” surface (d) and the topobaric surface (f) are initialized from the isopycnal (a) with $\bm{x}_{0}$ = ($180^{\circ}\mathrm{E}$, $0^{\circ}\mathrm{N}$) and reference depth $z_{0}=$-1988.60\text{\,}\mathrm{m} (the analogue of the reference pressure $p_{0}$ for the Boussinesq case).

Figure 4 shows the associated regions for the Reeb graph of $\utilde{z}$ for the $\sigma_{2}$-surface that initializes the topobaric surface calculation. This is part of step 2 in the first iteration to create the $\tau$-surface. Some large eddies can be seen as individual regions, and zonally elongated structures can be seen in the Southern Ocean. As there are over 40,000 arcs, the full structure is far beyond comprehension. It is possible to simplify the graph as in Section 4, but the present goal—to produce as neutral a surface as possible—is best met by keeping all the fine-scale structure of the Reeb graph.

5.4.2 Computation time

The topobaric surface (f) is constructed in five iterations, and the convergence is quite rapid: the root-mean-square change of $\utilde{z}$ after the first iteration is $15\text{\,}\mathrm{m}$, then $90\text{\,}\mathrm{c}\mathrm{m}$, $5.0\text{\,}\mathrm{c}\mathrm{m}$, $2.8\text{\,}\mathrm{m}\mathrm{m}$, $0.45\text{\,}\mathrm{m}\mathrm{m}$, after iterations two through five. On a (single) 2.2 GHz processor, the $\tau$-surface was computed in under 50 seconds. Using the original code, the $\omega$-surface requires about half a day (50 iterations). Vectorizing the original MATLAB code reduced this to about 30 minutes. The $\omega$-surface code remains under development, so a precise speed comparison is not the present goal.

Topobaric surfaces have the advantage that the Reeb graph (which is quick to compute) mostly decouples the problem. Each arc $a$ of the Reeb graph that is not on a cycle (which is the vast majority of arcs) requires fitting an affine linear function in $\hat{\pi}_{a}(z)$; this involves a dense $n$ by 2 matrix, where $n$ is the number of water columns in the region $\mathcal{R}_{a}$. Arcs on cycles are fit in a larger, coupled problem to satisfy the cycle constraint. This involves a block matrix with as many submatrices (dense $n$ by 2 matrices as above) as there are arcs in the cycle basis, arranged diagonally, plus one equality constraint for each cycle in the cycle basis. Then, updating $\utilde{z}$ is completely decoupled, being a single problem per water column, though non-linear.

In contrast, the bottleneck for $\omega$-surfaces is finding the least-squares solution of a coupled linear equation for the entire ocean. This involves a sparse $m$ by $n$ matrix, where $n$ is the number of water columns in the rectilinearly gridded ocean and $m\approx 2n$. This matrix is banded with five non-zero entries per row, plus a row of ones at the bottom to conserve density.

5.4.3 Neutrality

The neutrality of the six surfaces is now compared. In fact, neutrality is one of three desirable properties of a quasi-conservative variable (such as potential temperature in a dry atmosphere); the other two are material conservation and the existence of an exact geostrophic streamfunction [de Szoeke and Springer, 2009]. Material conservation cannot be assessed because topobaric surfaces are 2D, not 3D structures. Topobaric surfaces possess an exact geostrophic streamfunction [Stanley, 2019]. Thus, attention rests on neutrality.

The error from neutrality (27) is measured by

[TABLE]

which is zero for a perfectly neutral surface. Numerically, $\nabla\utilde{\rho}$ and $\nabla\utilde{z}$ are calculated by non-centred finite differences, and $\utilde{\rho_{z}}$ is evaluated from the equation of state using $\utilde{S}$, $\utilde{\theta}$, and $\utilde{z}$ averaged between the two grid points involved in the aforementioned finite difference. This is a third order accurate discretisation: expanding both terms in this discretisation of $\nabla\utilde{\rho}$ using a Taylor series about the averaged $\utilde{S}$, $\utilde{\theta}$, and $\utilde{z}$ reveals the quadratic terms cancel identically.

A second measure of error, and one for which we have a more familiar numeric sense, is the diapycnal diffusivity caused by the isopycnal diffusivity when the surface is not perfectly aligned with the neutral tangent plane. This is called the fictitious diapycnal diffusivity [McDougall and Jackett, 2005a, Klocker et al., 2009], expressed as

[TABLE]

where $K$ is an isopycnal eddy diffusivity, taken as a representative constant $K=1000~{}\mathrm{m}^{2}~{}\mathrm{s}^{-1}$, and

[TABLE]

is the slope difference between the neutral tangent plane and the approximately neutral surface. The slope of the neutral tangent plane, $\nabla_{n}z$, is a vector field in 3D space, hence the under-tilde is used to restrict it to the surface in question. Expressing $\bm{N}=(N_{1},N_{2},N_{3})$ in the $(\bm{\hat{i}},\bm{\hat{j}},\bm{\hat{k}})$ basis, an explicit formula is $\nabla_{n}z=-(N_{1}/N_{3})\bm{\hat{i}}-(N_{2}/N_{3})\bm{\hat{j}}$. A more useful expression for $\bm{s}$ is found as follows. Use the standard coordinate transformation $\nabla_{n}\theta=\nabla_{z}\theta+\partial_{z}\theta\,\nabla_{n}z$ and similarly for $\nabla_{n}S$. Here, $\nabla_{z}$ is a horizontal gradient at constant depth. Multiply these by $\rho_{\theta}$ and $\rho_{S}$ respectively, then sum to cancel the neutral $\theta$ and $S$ gradients, by (8). Repeat this for $\nabla\utilde{\theta}$ and $\nabla\utilde{S}$ (for which this cancellation is not complete). Subtracting the two results yields (dropping most under-tildes for visual clarity)

[TABLE]

where $N$ is the Brunt–Väisälä frequency, $N^{2}=-g\rho^{-1}\partial_{z}\sigma_{L}$, where $\sigma_{L}=B(S,\theta,\utilde{z})$ is the locally referenced potential density. Numerically, $D^{f}$ requires the components of $\bm{s}$ on the same grid, so now $\bm{\epsilon}$ in (31) is computed using centred differences for $\nabla\utilde{\rho}$ and $\nabla\utilde{z}$, and $\utilde{\rho_{z}}$ is computed from the equation of state using $\utilde{S}$, $\utilde{\theta}$, and $\utilde{z}$ averaged between the two adjacent water columns (maintaining the third order accuracy), not the central column (tempting as that is). Piecewise Cubic Hermite Interpolating Polynomials (PCHIPs) are used to evaluate $\partial_{z}\sigma_{L}$. As $D^{f}\propto N^{-4}$, global statistics of $D^{f}$ can be overwhelmed by a single grid point with $N\approx 0$, such as in mode water, so $N$ is artificially increased to a minimum value of $2\times 10^{-4}~{}\mathrm{s}^{-1}$.

Figure 5 maps $D^{f}$ for these six surfaces, with a common mask applied so $D^{f}$ is shown only where it is valid on all six surfaces. The $\sigma_{2}$-surface performs well (low $D^{f}$) where it is near the reference depth of $-2000\text{\,}\mathrm{m}$, such as in the Pacific and Indian Oceans, and to a lesser extent in the Atlantic Ocean; however, it performs very poorly in the Southern Ocean where it rises to the sea-surface. In contrast, the $\delta$-surface performs decently in the Southern Ocean where the salinity and potential temperature are close to their reference values, but poorly elsewhere, particularly in the North Atlantic. The reference values for these two surfaces can always be chosen so they perform well in a limited geographic region, but they struggle globally. The $\gamma^{n}$-surface overcomes this problem, yielding reasonably low $D^{f}$ globally. The underlying $4^{\circ}$ by $4^{\circ}$ grid of the WOCE atlas from which neutral density interpolates [Jackett and McDougall, 1997] does produce artefacts, as can be seen in the Southern Ocean where $D^{f}$ changes sharply. The $\sigma_{\nu}$-surface exhibits high $D^{f}$ in both the Southern Ocean and North Atlantic, but not so high as for the $\sigma_{2}$- or $\delta$-surfaces in these regions. Both the $\omega$-surface and $\tau$-surface exhibit very small $D^{f}$ globally, and both exhibit their largest $D^{f}$ in the Southern Ocean and North Atlantic, right where neutral helicity is largest [Klocker et al., 2009, Fig. 4c].

Whereas the $\omega$-surface spreads $D^{f}$ errors quite smoothly over the global ocean, the $\tau$-surface exhibits much more eddy-scale, filamentary structure in $D^{f}$. Some similar structure is visible on the $\sigma_{2}$-surface in the Pacific, where it is extremely neutral, suggesting this structure is partly real. However, it may also be caused by neighbouring grid points that are in regions associated with different arcs. This is a numerical difficulty arising from the finite difference underlying the calculation of $\bm{\epsilon}$: in the limit as the grid spacing goes to zero, $\bm{\epsilon}$ would not compare grid points across regions (except on contours through $\utilde{p}$ saddles). This issue is exacerbated when neighbouring grid points are in regions associated with arcs that are not incident to a common node in the Reeb graph, as can happen for filamentary regions. The topobaric algorithm could be conceivably modified to smooth these numerical errors out, for instance by inflating each associated region by one grid point in every direction when fitting $\hat{\pi}$. This has not been tested, but even without such algorithmic enhancements, the $\tau$-surface performs very well.

For a more quantitative comparison, Fig. 6 shows histograms of $\bm{\epsilon}$ and $D^{f}$ for the above six surfaces, and for a further set of these six surfaces that intersect ($180^{\circ}\mathrm{E}$, $0^{\circ}\mathrm{N}$, $-997.99\text{\,}\mathrm{m}$)141414 Again the $\omega$-surface is calculated first, initialized from a $\sigma_{1}$-surface intersecting ($180^{\circ}\mathrm{E}$, $0^{\circ}\mathrm{N}$, $-1000\text{\,}\mathrm{m}$). For the other surfaces, (a) $\sigma_{1}=1032.0053~{}\mathrm{kg}~{}\mathrm{m}^{-3}$, (b) $\delta=-1.4026\times 10^{-2}~{}\mathrm{kg}~{}\mathrm{m}^{-3}$ with $S_{\delta}=$34.3389\text{,}\mathrm{p}\mathrm{s}\mathrm{u} and $\theta_{\delta}=3.0272^{\circ}$C, (c) $\gamma^{n}=27.5400$, (d) the “orthobaric” surface uses a cubic spline for $\hat{\pi}$ with knots only at 0, -200, and $-6000\text{\,}\mathrm{m}$, and (d,f) the “orthobaric” and topobaric surfaces have reference location $\bm{x}_{0}$ = ($180^{\circ}\mathrm{E}$, $0^{\circ}\mathrm{N}$) and reference depth $z_{0}=$-997.99\text{\,}\mathrm{m}, and are initialized from the isopycnal (a). . The area-weighted $l_{1}$ (mean absolute error) and $l_{2}$ (root mean square error) norms are listed as inset tables in Fig. 6. The data underlying each histogram and norm is best thought of as a 1D array. For $D_{f}$, this is an array of all $D_{f}$ values at the grid cell centres within a common mask where $D_{f}$ is valid on all six surfaces intersecting a common depth at ($180^{\circ}\mathrm{E}$, $0^{\circ}\mathrm{N}$). For $\bm{\epsilon}$, this array is roughly twice as long, containing the zonal component of $\bm{\epsilon}$ within a similarly constructed common mask, followed by the meridional component of $\bm{\epsilon}$ within a similarly constructed common mask.

First consider the surfaces at roughly $-2000\text{\,}\mathrm{m}$ in the Pacific (bottom row). Compared to the $\sigma_{2}$-surface, the $\gamma^{n}$-surface reduces the metrics $||\bm{\epsilon}||_{1}$, $||\bm{\epsilon}||_{2}$, $||D^{f}||_{1}$, and $||D^{f}||_{2}$ by factors of 7.9, 6.9, 107, and 243, respectively. Whereas $||\bm{\epsilon}||_{1}$ involves a sum of zonal and meridional components of $\bm{\epsilon}$, $||D^{f}||_{1}$ involves a sum of $\bm{\epsilon}\cdot\bm{\epsilon}=\epsilon^{2}$, and hence punishes surfaces with more heterogeneous errors. Indeed, the $\sigma_{2}$-surface’s poor performance in the Southern Ocean is particularly crippling to its $D^{f}$ metrics. This punishment becomes even more severe for $||D^{f}||_{2}$, which involves sums of ${\epsilon}^{4}$. The weighting $g/(\rho N^{2})$ actually lowers $D^{f}$ for the $\sigma_{2}$-surface: the stratification $N^{2}$ is larger in the Southern Ocean as these surfaces rise towards the sea-surface (they tend to avoid low $N^{2}$ mode water, which occupies little space in density coordinates). The $\tau$-surface further improves upon the $\gamma^{n}$-surface, reducing $||D^{f}||_{2}$, for instance, by a factor of 29. Also, the $\tau$-surface, with its geographic dependence, vastly outperforms its geographically independent cousin, the $\sigma_{\nu}$-surface: $||D^{f}||_{2}$ differs by a factor of 448. Still, the $\omega$-surface performs best, with $||\epsilon||_{2}$ and $||D^{f}||_{2}$ smaller by factors of 4.3 and 5.9, respectively, than those of the $\tau$-surface.

This analysis does not assess which of neutral density and orthobaric density is superior, a subject of some debate [McDougall and Jackett, 2005a, de Szoeke and Springer, 2009]. Indeed, “orthobaric” surfaces are computed here essentially from the topobaric algorithm, rather than from the algorithm (and dataset) of de Szoeke et al. [2000]. Moreover, material conservation has not been evaluated. Topobaric surfaces, though, take the best of both worlds: they contain geographic dependence which enables neutral density to be more neutral, while retaining the good theoretical properties of orthobaric density.

Ideally, the fictitious diapycnal diffusivity is less than the true diapycnal diffusivity. Taking a representative value of $10^{-5}~{}\mathrm{m}^{2}~{}\mathrm{s}^{-1}$ for the latter, the fraction of ocean over which this is false, for surfaces (a) through (f), is 4.9%, 7.7%, 1.0%, 8.5%, $2.3\times 10^{-5}$, and $8.5\times 10^{-4}$, respectively. Indeed, the histograms of Fig. 6 have a long tail towards high $D^{f}$ for the $\sigma_{2}$- and $\delta$- and $\sigma_{\nu}$-surfaces. The area where $D^{f}$ exceeds this threshold is very small for both $\omega$- and $\tau$-surfaces.

The comparison between the $\omega$- and $\tau$-surfaces is similar on the surfaces around $-1000\text{\,}\mathrm{m}$ in the Pacific, with $||\bm{\epsilon}||_{1}$, $||\bm{\epsilon}||_{2}$, $||D^{f}||_{1}$, and $||D^{f}||_{2}$ smaller for the $\omega$-surface by factors of 2.8, 4.5, 7.8, and 2.5, respectively. In terms of these metrics, the $\gamma^{n}$ surface performs worst at this depth, but this is due to a block of troublesome points south of Java, where the depth of the $\gamma^{n}$-surface changes by over $-1000\text{\,}\mathrm{m}$ between neighbouring grid points. This unrealistic behaviour is caused by the neutral density software’s underlying $4^{\circ}\times 4^{\circ}$ grid, which must have had different bathymetry. The histograms in Fig. 6 reveal the $\gamma^{n}$-surface, aside from these points, again outperforms the $\sigma_{1}$-, $\delta$-, and $\sigma_{\nu}$-surfaces. In general, neutral density performs somewhat worse than it could because the Levitus [1982] climatology underlying neutral density differs from the ECCO2 state, and neutral density uses an older equation of state [Millero et al., 1980].

As $\omega$-surfaces minimize $||\bm{\epsilon}||_{2}$, there can be no hope151515 Actually, $\omega$-surfaces minimize $||\bm{\epsilon}||_{2}$ without changing its two-dimensional curl in the surface. This approach rests on an approximate relationship derived by McDougall and Jackett [1988] that relates the 2D curl of $\bm{\epsilon}$ in an approximately neutral surface to $HN^{-2}$. To the extent that this approximation is good, the 2D curl of $\bm{\epsilon}$ is set by the ocean hydrography, hence $\omega$-surfaces do not attempt to change it: the Klocker et al. [2009] algorithm takes an initial surface with an initial $\bm{\epsilon}$, and finds $\Phi^{\prime}$ to minimize $||\bm{\epsilon}+\nabla\Phi^{\prime}||_{2}$. The quality of this approximation determines the degree to which methods that vary the 2D curl of $\bm{\epsilon}$ might discover surfaces with smaller $||\bm{\epsilon}||_{2}$.

to produce an approximately neutral surface with smaller $||\bm{\epsilon}||_{2}$. Nonetheless, topobaric surfaces are quick to compute and perform admirably, while also also possessing an exact geostrophic streamfunction [Stanley, 2019].

6 Conclusions

The pathways along which the ocean, below the mixed layer, is connected are largely determined by neutral surfaces. The typical oceanic epineutral diffusivity is $\mathcal{O}(10^{3}~{}\mathrm{m}^{2}~{}\mathrm{s}^{-1})$ while dianeutral diffusivity is $\mathcal{O}(10^{-5}~{}\mathrm{m}^{2}~{}\mathrm{s}^{-1})$ in the main thermocline [MacKinnon et al., 2013]. It is therefore crucial to orient the large epineutral diffusion correctly, in the neutral tangent plane, lest a component of it act dianeutrally and swamp the real dianeutral diffusion. This can be quantified as a fictitious dianeutral diffusivity [McDougall and Jackett, 2005b, Klocker et al., 2009]. Relatedly, the lateral velocity that acts in the neutral tangent plane must pierce any approximately neutral surface, creating a flow across that surface. Klocker and McDougall [2010] quantified this flow, finding upwards of $10\text{\,}\mathrm{S}\mathrm{v}$ globally through deep $\sigma_{0}$-surfaces. This diapycnal flow is entirely fictitious, due to $\sigma_{0}$-surfaces’ poor alignment with the neutral tangent plane. They find this is reduced to generally less than $1\text{\,}\mathrm{S}\mathrm{v}$ when using $\omega$-surfaces.

The neutral tangent plane is parallel to a surface of locally referenced potential density [McDougall, 1987a]. Depth-level ocean models long ago rotated their diffusion tensors [Veronis, 1975, Redi, 1982] to align with the neutral tangent plane. This is fairly straightforward because it is a local problem—the neutral tangent plane is well-defined essentially everywhere. However, layered models require a quasi-conservative density variable for their vertical coordinate, which is a global problem, and one that is not well-defined: neutral tangent planes cannot be globally stitched together to form a well-defined neutral surface, because of non-zero neutral helicity [McDougall and Jackett, 1988]. As such, well-defined surfaces can only be approximately neutral. Many oceanic theories and analyses operate on neutral surfaces, or would do so if neutral surfaces were well-defined surfaces with a unique depth. We have shown that this is not guaranteed even in an ocean with zero neutral helicity, as neutral helices can exist around islands and other holes in a neutral surface.

Nowadays, a host of approximately neutral surfaces are available. Many of them are purely thermodynamic variables (functions only of salinity, temperature, and pressure), such as potential density [Wüst, 1935], specific volume anomaly [Montgomery, 1937], orthobaric density [de Szoeke et al., 2000], a rational approximation of neutral density [McDougall and Jackett, 2005b], and thermodynamic neutral density [Tailleux, 2016]. Global isosurfaces of these variables are limited in their neutrality not so much by the non-zero helicity of the real ocean, but overwhelmingly by the fact that they lack any geographic dependence. That is, even in an ocean with zero neutral helicity, a conservative density variable must be a function of latitude and longitude, as well as salinity, temperature, and pressure. Neutral density [Jackett and McDougall, 1997] is such a function, explicitly containing geographic information. Geographic dependence is implicitly built into $\omega$-surfaces [Klocker et al., 2009] by the rectilinear grid on which it operates. As such, $\gamma^{n}$- and $\omega$-surfaces can be close to neutral, globally.

The importance of geography arises because there is a multivalued functional relationship between in-situ density and pressure (equivalently, between practical/Absolute salinity and potential/Conservative temperature) on a neutral surface. This multivalued function has single-valued branches within certain geographic regions, and these branches differ between regions. This was well-known, but the shape (topology) of these geographic regions was unknown.

Much emphasis has been placed on these functional relationships differing between the Northern and Southern hemispheres [McDougall and Jackett, 2005a]. de Szoeke and Springer [2005] advanced the original orthobaric density to use a different virtual compressibility for the North and the South Atlantic. However, the result is discontinuous at the equator, and a fictitious force must be applied to a water parcel crossing this discontinuity to keep it on the same orthobaric density surface [de Szoeke and Springer, 2005]. Following this idea to its limit, de Szoeke and Springer [2009] developed “extended orthobaric density” by segmenting the Atlantic into arbitrarily many latitude bins; this minimizes the discontinuities, but creates more of them. Somewhat similarly, patched potential density [Reid and Lynn, 1971] and generalized patched potential density [Tailleux, 2016] add geographic dependence to potential density by segmenting the ocean into boxes aligned with latitude circles, longitude circles, and depth or pressure levels. However, these are not the correct geographic shapes underpinning the multivalued functional relationship between in-situ density and pressure on neutral surfaces.

The cause of the geographic dependence may also have been unknown. McDougall and Jackett [2005a] correctly stated that, on a neutral surface having pressure $\utilde{p}$, new branches of the functional relationship between in-situ density and pressure open up at points where $\nabla\utilde{p}=\bm{0}$. However, thinking in terms of neutral trajectories in a zonally uniform ocean, they stated that this occurs at extrema of $\utilde{p}$. Extrema of $\utilde{p}$ do have $\nabla\utilde{p}=\bm{0}$, but it is more appropriate to say these points close branches of the multivalued function, and branches open at saddle points of $\utilde{p}$. (Branches often exist that do not enclose any extrema, such as the red $AB$ arc in Fig. 2. A better example is the cycle in Fig. 2, bearing no relation at all to $\utilde{p}$ extrema. If $\utilde{p}$ has one maxima and one minima and $n$ islands, there can be up to $1+3n$ branches.)

The first major advance of this paper is to reveal the entire geographic structure underlying this multivalued function, by use of the Reeb graph of $\utilde{p}$. Each such region is mapped to an arc of the Reeb graph, and nodes represent the critical points (saddles and extrema) of $\utilde{p}$. The structure of the graph—which arcs are incident to which nodes—determines how the geographic regions nest into a global structure. Relationships between branches of the multivalued function are determined by the structure of the graph.

Topobaric surfaces represent the second major advance of this paper. Knowing the geographic regions underlying different branches of the multivalued relation between in-situ density and pressure, these branches can be fit empirically, subject to some matching conditions related to cycles in the Reeb graph. A root-finding problem ensues, solving for the pressure in each water column for which the in-situ density at that pressure exactly matches the single-valued function at that pressure. In this way, an iterative procedure turns an isopycnal (or any approximately neutral surface) into a topobaric surface. Topobaric surfaces have an exact multivalued functional relation between in-situ density and pressure, are very close to neutral, and possess an exact geostrophic streamfunction [Stanley, 2019].

Topobaric surfaces are the topologically correct extension of isosurfaces of orthobaric density to have geographic dependence. Orthobaric density is pycnotropic—a function only of pressure and in-situ density. It employs a pycnotropic virtual compressibility that approximates the real compressibility. On an isosurface of orthobaric density, the in-situ density is an implicit function of pressure, so the virtual compressibility is a function only of pressure, analogous to how topobaric surfaces approximate the real compressibility161616We call $\rho_{p}$ the compressibility, rather $\rho^{-1}\rho_{p}$, a convention also used by de Szoeke et al. [2000]., $\utilde{\rho_{p}}$, by $\hat{\pi}(\utilde{p})$ using a multivalued function $\hat{\pi}$. The discontinuities of extended orthobaric density [de Szoeke and Springer, 2009] are caused by changing the virtual compressibility at latitudinal boundaries. No such discontinuities exist for topobaric surfaces, because the branches of $\hat{\pi}$ change at contours of constant $\utilde{p}$. Future work aims to develop a 3D topobaric density variable.

This topological musing bears on the Lorenz convention, which sets the depth of a surface that has outcropped or incropped to be infinitesimally below the sea-surface or above the sea-floor [Young, 2012]. In this sense, the only holes in the surface are islands—quite a simplification. Is this justifiable? From the topological perspective, two contours of the same level set become one contour that snakes around the hole; points on either side of the hole on the same level set are glued together in topological space. But the fundamental reason for alignment of contours of salinity and potential temperature on a neutral surface is because both tracers are materially advected by the flow, of which there is none through a hole. The Lorenz convention makes sense from the perspective of one fluid column, considering also the density of air and solid Earth, but not from a broader perspective. For example, different water masses exist on different sides of submarine ridges (and so on to smaller features), and it would be wrong to join them.

Using the Reeb graph as a (computational) tool to study multivalued functions arising from equations like (9) is a new endeavour. This is highly translatable to other problems. A multivalued functional relationship between variables is a powerful idea, so fresh insights on other problems may soon be discovered by similar means.

Acknowledgements

The author thanks Chris Hughes and David Marshall for helpful discussions, Trevor McDougall for early encouragement on related work that ultimately led to this paper, two reviewers who helped clarify this paper, Andreas Klocker for sharing his $\omega$-surface software, David Gleich for his graph theory MATLAB toolbox (GAIMC), and Harish Doraiswamy and Vijay Natarajan for their ReCon software. The author was supported by the Clarendon Scholarship, and the Canadian Alumni Scholarship at Linacre College, University of Oxford. MATLAB software to compute topobaric surfaces is available from the author’s website.

Appendix A Graph theory glossary

A graph is a tuple $G=(N,A)$, where $N$ is a set of nodes and $A$ is a set of arcs, and elements of $A$ are subsets of $N$ with precisely two elements. Intuitively, a graph is a collection of nodes, and arcs connecting any two nodes.

When $A$ is a multiset (its elements need not be unique), $G$ is a multigraph. A multigraph can have multiple arcs incident upon the same pair of nodes. The Reeb graph is actually a multigraph, but “graph” is used as shorthand.

Arc $a$ is incident upon node $n$ when $n\in a$.

Two nodes are adjacent if there is an arc incident upon both nodes.

Two arcs are incident if they are both incident upon a common node.

The degree of a node is the number of distinct arcs incident upon it.

A walk is an ordered sequence that alternates between nodes and arcs such that every arc is incident upon both nodes next to it in the sequence. In standard graph theory, a walk must start and end with nodes, but we relax that here.

A (simple) cycle is a walk whose first node is also its last node, and otherwise contains no repeated vertices, and contains no repeated arcs.

A connected graph is a graph such that there is a walk between any two of its nodes.

A tree is a connected graph with no cycles.

A spanning subgraph of a graph $G=(N,A)$ is a graph $H=(N,F)$ with $F\subseteq A$.

A minimum spanning tree is a spanning subgraph that is also a tree.

An Eulerian graph is a graph whose nodes all have even degree.

The cycle space of a graph $G$ is the set of all Eulerian spanning subgraphs of $G$. Two elements of the cycle space may be added by symmetric difference of their arc sets to produce a new element of the cycle space. A cycle basis is a set of simple cycles, a subset of which can be combined by symmetric difference on their arc sets to produce any element of the cycle space.

Appendix B Neutral helix pitch around islands

A neutral trajectory that returns to its starting water column does so, in general, at a different depth from which it began. Is this change—the pitch of a neutral helix—fundamentally different depending on whether the enclosed area is open ocean, or an island/seamount? A preliminary analysis towards this question is shown here, using ECCO2 data [Menemenlis et al., 2005] from 22–24 December 2002.

Given one water column and an initial depth $z_{A}$, a neutral trajectory to a (nearby) water column ends at depth $z_{B}$ if

[TABLE]

where $\overline{z}=(z_{A}+z_{B})/2$, and where $S_{i}(z)$ and $\theta_{i}(z)$ are the salinity and potential temperature as functions of depth in water column $i\in\{A,B\}$ [Jackett and McDougall, 1997]. The functions are taken as piecewise linear interpolants, and (32) is solved by bisection, geometrically expanding outward from an initial guess of $z_{A}$ until a sign change is found.

A neutral trajectory around an island/seamount is discovered by starting just east of an island at a particular depth, then making neutral trajectories one grid cell at a time. The forward direction is initially west. A neutral trajectory in the forward direction is tested; if this fails (no solution, meaning the neutral trajectory grounds or outcrops) the forward direction is rotated $90^{\circ}$ clockwise. This is repeated until a neutral trajectory in the forward direction succeeds. That step is made, and the forward direction is rotated $90^{\circ}$ counter-clockwise. The whole procedure stops when the trajectory returns to the starting water column (or an upper bound of steps is reached, in which case the neutral trajectory has spiralled up or down the topographic slope, and this island must be tested at another initial depth).

Figure 7(a) shows that a neutral helix around Flemish Cap, starting at $-400\text{\,}\mathrm{m}$, descends by $7.8\text{\,}\mathrm{m}$ in one counter-clockwise loop. This is compared to other neutral helices starting at $-400\text{\,}\mathrm{m}$ but at other longitudes and having the same horizontal shape (and thus the same area) as that around Flemish Cap. Most of these open-ocean neutral helices have a pitch very close to $0\text{\,}\mathrm{m}$; the maximum is $1.93\text{\,}\mathrm{m}$. Thus, it seems that Flemish Cap is exceptional amongst duplicates of Flemish Cap’s shape in the open ocean.

Figure 7(b)—(j) show the results of this analysis repeated at other islands/seamounts and starting depths. Moving just $100\text{\,}\mathrm{m}$ deeper, Flemish Cap again appears to be exceptional. Australia and New Guinea at $-700\text{\,}\mathrm{m}$ appears to be exceptional, but not so much so at $-400\text{\,}\mathrm{m}$: a not-insignificant pitch of $-0.52\text{\,}\mathrm{m}$ still exists, but this is smaller than for neutral helices of the same shape at other longitudes, which is not hard given their large area. Fiji at $-650\text{\,}\mathrm{m}$ appears to be exceptional, but not so at $-500\text{\,}\mathrm{m}$. The neutral helix pitch around Madagascar at $-500\text{\,}\mathrm{m}$ and around Kerguelen at $-750\text{\,}\mathrm{m}$ are both tiny. Around Madagascar at $-1000\text{\,}\mathrm{m}$ the neutral helix pitch is larger than equivalent helices’ pitches in the Pacific and Atlantic Oceans, but not in the Indian Ocean. Neutral helices around Kerguelen starting at $-2000\text{\,}\mathrm{m}$ appear to have a moderately exceptional pitch.

Not all islands/seamounts should necessarily be exceptional in this way—only those which produce a cycle in the Reeb graph of $\utilde{z}$. A sense of which islands/seamounts create cycles is gained by looking at depth contours on a potential density surface referenced to the starting depth of the neutral trajectory (not shown). Flemish Cap, Australia and New Guinea, and Fiji should clearly create cycles, because $\utilde{z}$ contours are found emanating from them and intersecting land elsewhere. This is also true of Madagascar at $-1000\text{\,}\mathrm{m}$ and Kerguelen at $-2000\text{\,}\mathrm{m}$, around both of which the neutral helix pitch is large. However, for Madagascar at $-500\text{\,}\mathrm{m}$, a $\utilde{z}=$-450\text{,}\mathrm{m} contour very nearly encircles Madagascar—it barely intersects the northern end of Madagascar, skims the coast of Africa intersecting it only slightly in four places, and just intersects the southern tip of Nazareth Bank to the east of Madagascar. There are a few small seamounts contained within this $\utilde{z}=$-450\text{\,}\mathrm{m} contour, which can create cycles around Madagascar. Still, there is a hint here that Madagascar should not give rise to the path-dependency responsible for neutral helices, which would explain why the pitch of a neutral helix around Madagascar at $-500\text{\,}\mathrm{m}$ is so small, only $-4\text{\,}\mathrm{c}\mathrm{m}$. Similarly, for Kerguelen at $-750\text{\,}\mathrm{m}$, the $\utilde{z}=$-500\text{,}\mathrm{m} and $\utilde{z}=$-700\text{\,}\mathrm{m} contours that flank Kerguelen happen to traverse the entire Southern Ocean without intersecting land (or nearly so—the northern one does skim five small seamounts, which collectively occupy only 35 grid points). Again, this is suggestive that neutral helices around Kerguelen at $-750\text{\,}\mathrm{m}$ should have a pitch near zero, in agreement the numerical calculation of a pitch of only $2\text{\,}\mathrm{c}\mathrm{m}$. However, it is not clear why the neutral helix around Fiji starting at $-500\text{\,}\mathrm{m}$ is so small, as Fiji should create a cycle in the Reeb graph.

The depths, and to some extent the islands, in Fig. 7 were chosen with some care, to illustrate exceptional and unexceptional cases. The actual results seem fairly sensitive to the initial location. For example, the pitch of neutral helices around Madagascar as a function of initial depth is a fairly complicated function; it is generally increasing with depth, but has six $0\text{\,}\mathrm{m}$ crossings shallower than $-1250\text{\,}\mathrm{m}$ and seven $1\text{\,}\mathrm{m}$ crossings between $-1250\text{\,}\mathrm{m}$ and $-2250\text{\,}\mathrm{m}$ (not shown).

Are islands/seamounts—in particular those that intersect depth contours of approximately neutral surfaces so as to create cycles in the Reeb graph—exceptional at producing neutral helices with large pitches? Further work is needed to definitively say, but this preliminary analysis suggests it is a possibility.