VENM: An Algorithm to Accurately Calculate Neutral Slopes and Gradients.

Mesoscale eddies stir along the neutral plane, and the resulting neutral diffusion is a fundamental aspect of subgrid-scale tracer transport in ocean models. Calculating neutral diffusion traditionally involves calculating neutral slopes and three-dimensional tracer gradients. The calculation of the neutral slope traditionally occurs by computing the ratio of the horizontal to vertical locally referenced potential density derivative. However, this approach is problematic in regions of weak vertical stratification, prompting the use of a variety of ad hoc regularization methods that can lead to rather nonphysical dependencies for the resulting neutral tracer gradients. Here we use a VErtical Non-local Method "VENM," a search algorithm that requires no ad hoc regularization and significantly improves the numerical accuracy of calculating neutral slopes, neutral tracer gradients, and associated neutral diffusive fluxes. We compare and contrast VENM against a more traditional method, using an independent objective neutrality condition combined with estimates of spurious diffusion, heat transport, and water mass transformation rates. VENM is more accurate, both physically and numerically, and should form the basis for future efforts involving neutral diffusion calculations from observations and possibly numerical model simulations. ©2019. The Authors.

Introduction

In the ocean, the neutral direction is the direction along which fluid elements can move a small distance without experiencing a buoyant restoring force. The idea that ocean properties are advected and mixed predominantly along neutral directions has its origins in empirical observations dating back at least to Iselin (1939), who was led to this insight from the similarities between geographically separated vertical casts in a salinity‐temperature diagram. Recent presentations such as in section 7.2 of Griffies (2004), section 2 of McDougall and Jackett (2005), and section 1 of McDougall et al. (2014) argue that the very small dissipation of turbulent kinetic energy measured in the ocean interior by microstructure measurements supports the relevance of orienting tracer mixing according to neutral and dianeutral directions. It is assumed that mesoscale turbulence produces little small‐scale dissipation; thus, mesoscale mixing is constrained to occur predominantly in the neutral plane.

When mesoscale processes are not resolved, as in coarse‐resolution ocean climate models, the associated tracer mixing must be parameterized. The calculation of the neutral direction (neutral slopes) forms a fundamental element of the subgrid‐scale parameterizations used in such models, as well as in observation‐based watermass analyses. Neutral slopes arise in the formulation of both the advective (or “skew”) component of parameterized mesoscale transport (Gent & McWilliams, 1990; Griffies, 1998; McDougall & McIntosh, 2001), as well as the neutral diffusive component (Griffies et al., 1998; Redi, 1982; Solomon, 1971). Evidence from simulations has shown that numerical ocean models are rather sensitive to details of both the advective (Gnanadesikan, 1999; Hirst et al., 1996; Marshall et al., 2017) and diffusive components (Gnanadesikan et al., 2015; Pradal & Gnanadesikan, 2014) of mesoscale transport. An accurate representation of the neutral direction is therefore essential for an accurate implementation of the related ocean physics.

The Importance of Neutral Slopes

Neutral directions are the natural generalization of isopycnal surfaces, with the neutral direction tangent to the locally referenced potential density (ρ ) surface (McDougall, 1987a). Due to a nonzero neutral helicity, which is an effect of the nonlinear equation of state, the neutral tangent plane (NTP) is always well defined locally, while the extensive neutral surfaces have a nontrivial topology (McDougall & Jackett, 1988; Stanley, 2019). The resulting unit normal to the NTP, or dianeutral unit vector, is defined according to the gradient of locally referenced potential density ρ (McDougall et al., 2014) where the slopes of the NTP are given by Here S A is the Absolute Salinity McDougall et al. (2009), Θ the Conservative Temperature (Graham & McDougall, 2013; McDougall, 2003), and the locally referenced potential density is computed as in situ density with respect to the local pressure p ; α=α(S A,Θ,p) is the thermal expansion coefficient, and β=β(S A,Θ,p) is the haline contraction coefficient with units of K −1 and , respectively.

Traditional numerical methods use a finite‐difference discretization based on the local grid stencil (neighboring grid values) to calculate the the normal direction n of equation (2) as the basis for calculating neutral directions and related physical parameterizations (Figure 1; Griffies et al., 1998; Madec, 2015; Marshall et al., 1997; Shchepetkin & McWilliams, 2005). Due to the local nature of the traditional method of calculation, we refer to this as the “local method.” The local methods have nontrivial numerical issues related to this calculation. Due to the key role of eddy‐induced neutral stirring and neutral diffusion for physical parameterizations of unresolved tracer transport, we conjecture that the same issues and sensitivities arise for parameterized tracer transport in prognostic ocean models.

Figure 1

Conceptual explanation of the vertically nonlocal method and local‐method for computing neutral slopes and gradients. The red dots on the south and north casts represent the tracer points at which data are provided. The yellow stars are the T‐grid data averaged on the vertical interface between the two casts. The purple square is where we wish to obtain values of the neutral slope and gradient. The black surface is an example of how a neutral surface could be oriented, which would lead to N 2≈0, and vertically nonlocal method obtains the approximate neutral surface as shown in green and avoids the problems with a locally small N 2.

Content of this Paper

Here we use VErtical Non‐local Method (VENM; section 2.2), a vertically nonlocal method for calculating neutral slopes and tracer gradients based on a search algorithm, rather than on the direct calculation of equation (2). VENM is based upon a method (and code) that has been used by Jackett and McDougall (1997) to construct Neutral Density γ n‐surfaces (Appendix A) and later by Klocker et al. (2009) to construct ω surfaces. In this paper, VENM provides estimates of neutral slopes and gradients locally, and not a means to estimate global neutral surfaces. Additionally, we here reveal for the first time the significant impact on ocean physics of adopting a VENM‐based method rather than a more traditional local method. These results directly impact ocean data analyses and, by conjecture, numerical ocean modeling.

The different methods are tested using a gridded climatology (section 4.1) and subsequently compared and contrasted against an independent neutrality condition (section 4.4). Thereafter, we compare the resulting estimates of water mass transformation (section 4.5) and heat transport (section 4.6). We show that VENM is significantly more skilled, both physically and numerically, compared to the local methods.

We also embedded a comparison to slopes and gradients along γ n surfaces and potential density surfaces referenced to 2,000 dbar, given by  kg/m3. The results for σ 2 and γ n are computed using an algorithm based on VENM that allows us to estimate the effects of nonneutrality when choosing σ 2 or γ n as an approximation to a neutral surface (Appendix B). This issue is important for numerical isopycnal modeling, which typically use σ 2 for the vertical coordinate, or data‐based analyses, which often use γ n as an estimate for the neutral direction.

Methods for Calculating Neutral Diffusion

The Local Method

Neutral diffusion is a fundamental part of the parameterization of subgrid‐scale mixing in models that do not resolve transient mesoscale eddies. A central component of neutral diffusion is a rotated diffusion operator that aligns tracer gradients according to neutral directions (Griffies et al., 1998; McDougall & Church, 1986; McDougall et al., 2014; Redi, 1982; Solomon, 1971; Veronis, 1975). The resulting downgradient tracer flux is given by with K the neutral diffusivity (generally a tensor), C the tracer concentration, and ∇N the gradient operator aligned according to the neutral direction. We use the notation to indicate “along the neutral tangent plane (NTP).” How one computes ∇N C has nontrivial implications for the resulting diffusion operator. We can write the neutral tracer gradient in the form of an operator that projects out that portion of the three‐dimensional tracer gradient aligned with the NTP (Griffies et al., 1998; McDougall et al., 2014; Olbers et al., 1986), in which case Expanding equation (4) into a matrix‐vector form leads to the following expression for the downgradient neutral diffusion tracer flux which reveals the role of the Redi diffusion tensor N acting on the three‐dimensional tracer gradient ∇C (Redi, 1982). We made the common assumption that the neutral diffusivity is isotropic in the NTP (see Smith & Gent, 2004  and Fox‐Kemper et al., 2013,  for discussions of anisotropic neutral diffusion). The three‐dimensional convergence of this downgradient diffusive flux then renders the neutral diffusion operator acting on the tracer C

For the numerical implementation of equation (6) it is common practice to use the small‐slope approximation to N of equation (5) (Gent & McWilliams, 1990). In the small‐slope approximation, it is assumed that , meaning that the angle between the neutral direction and the horizontal is small, leaving

It is important to retain the (3,3) component in equation (7), as it combines with the vertical diffusivity arising from dianeutral processes. The small‐slope approximation is motivated both numerically and physically. Numerically, the small‐slope approximation is less expensive to calculate. Physically, it is motivated since mesoscale eddies generally mix horizontally rather than neutrally when encountering regions of near vertical neutral directions, such as in the mixed layer, in which case neutral diffusion reverts to horizontal diffusion (Danabasoglu et al., 2008; Ferrari et al., 2008; 2010; Treguier et al., 1997).

McDougall et al. (2014) showed that the small‐slope approximation to the full tensor is not exactly in the direction of the correct neutral tracer gradient, but the related error is tiny and of negligible physical consequence. The neutral tracer gradient is then given by where

The three‐dimensional convergence of this downgradient small‐slope diffusive flux then renders the neutral diffusion operator acting on the tracer C

In equation (8) we used the two‐dimensional projected nonorthogonal neutral tracer gradient ∇n C, which was first introduced by Starr (1945) and has since been widely used in geophysical fluid theory and modeling. Isopycnal ocean models (Bleck, 1978a, 1978b) compute isopycnal tracer diffusion by taking the two‐dimensional convergence of −K ∇n C along isopycnal layers. This approach contrasts to that of geopotential coordinate models, which compute the 3D‐convergence of as in equation (10). The two approaches are mathematically identical under the small‐slope approximation, and so the different approaches arise from convenience rather than fundamentals. Whether one uses ∇n C or ∇N C, the integrity of a numerical neutral diffusion scheme relies on an accurate representation of the neutral slope, with the slope calculation the focus of this paper.

Fundamental to the Redi tensor is the computation of the slope of the NTP relative to the horizontal plane (S and S ). The local method calculates the slope of the NTP, or “neutral slope,” as the ratio of horizontal to vertical derivatives of the locally referenced potential density (equation (2)). In practice the calculation of the neutral slope according to equation (2) can be problematic, with problems arising from division by ∂ρ l/∂z in regions where the vertical stratification is weak (Figure 1). This division results in an unbounded neutral slope, thus requiring a numerical regularization such as those proposed by Cox (1987), Gerdes et al. (1991), and Danabasoglu and McWilliams (1995). Depending on the regularization chosen, one can realize an overly small or overly large neutral tracer gradient. Numerical simulations are sensitive to the (arbitrary) choice for the slope regularization (e.g., Gnanadesikan et al., 2007; Griffies et al., 2005; Ferrari et al., 2008), thus leading to a rather unfortunate sensitivity to a generally ad hoc numerical scheme. In addition, slope regularizations such as that proposed by Cox (1987) can lead to spurious dianeutral transports (though the methods of Gerdes et al., 1991, and Danabasoglu & McWilliams, 1995, avoid this problem).

VENM

The use of ad hoc slope regularizations is unsatisfying both physically and numerically. We thus seek an alternative approach to calculate neutral slopes and related neutral tracer gradients. For this purpose, rather than starting from equation (2) for the neutral slope, we start from the mathematically identical form

This expression for the slope is based on computing the change in depth of a neutral direction when moving laterally along a NTP; hence, indicates “along the NTP, at constant y.” Correspondingly, neutral tracer gradients can also be obtained by directly computing the change in tracer values when moving (laterally) along a NTP following equation (8), leaving

We make use of a search algorithm to calculate the slope of the neutral direction (equation (11)) and the neutral tracer gradients (equation (12)). For a finite resolution data set, such as from a numerical model or an observational analysis, the goal of the algorithm is to find the intersection of a NTP between adjacent vertical casts (Figure 1; details are given in section 3.2). Knowing the intersection of one NTP on two adjacent vertical casts allows us to estimate the height and tracer values at these locations, which can subsequently be used to calculate both the neutral slopes and neutral tracer gradients. Notably, the search generally spans multiple vertical grid cells, which contrasts to the vertically local methods based on the slope calculation starting from equation (2). Very steep NTPs may intersect the surface or bottom before reaching the adjacent column, with these NTPs effectively vertical.

In the following sections, we compare the three methods defined below to calculate neutral tracer gradients. We quantify and compare the related sensitivity of diagnosed physical processes, such as vertical heat transport and water mass transformation.

The local method. For the local method, equation (8) is the basis for a discrete approximation of the tracer gradients along NTPs. Both the calculation of the slopes (according to equation (2)) and tracer gradients are based on the local grid stencil discretization of the vertical and horizontal derivatives. The neutral slopes and tracer gradient of the local method are represented as

VENM: a vertically nonlocal method. VENM uses an algorithm that finds the intersections of a NTP between two vertical casts. The difference between tracer values and height at both intersections provides a direct estimate of the neutral tracer gradient and slope according to equations (11) and (12). The details of this algorithm are described in section 3.2. VENM can reach vertically to distant grid cells for estimating the NTP, thus earning the nonlocal moniker. The neutral slopes and tracer gradient of VENM are represented as and S nloc, respectively. Although S nloc is not needed to calculate neutral gradients, it can be used for other purposes such as parameterizing eddy‐induced stirring (e.g., Gent et al., 1995).

The hybrid method. The hybrid method employs equation (8) with the small‐slope approximation, to calculate the neutral gradients. However, it combines neutral slopes resulting from VENM (rather than the local method from equation (2)), with Cartesian tracer gradients calculated using the local grid stencil. The neutral slopes and tracer gradient of the hybrid method are represented as , where uses S nloc. The hybrid method serves to illustrate whether the local method suffers only from numerical issues arising from calculating slopes according to equation (2) or also suffers from the introduction of local extrema in Cartesian tracer gradients that can lead to numerical errors. In other words, it assesses if accurate slopes alone are enough to avoid the issues arising in the local method. Note that due to the independent means by which the slopes and Cartesian tracer gradients are obtained, the hybrid method is not self‐consistent. That is, there may be small amounts of unphysical diffusion as a result of this inconsistency as detailed in Griffies et al. (1998). We thus introduce this method only as a halfway point between the traditional fully local method and VENM.

Slope Regularization

Coarse‐resolution prognostic ocean models commonly make use of neutral diffusion schemes with neutral slopes based on equation (2), that is, the local method. For these models, one generally sets the maximum neutral slope according to numerical time stepping stability constraints rather than on numerical accuracy considerations (Cox, 1987; Griffies, 1998; Lemarié et al., 2012). However, numerical stability does not imply numerical accuracy. Indeed, the calculation of neutral slopes using the local approach is inaccurate for slopes steeper than the grid aspect ratio Δz/Δx. The reason is that the method makes use only of information that is vertically local and so neutral slopes steeper than the grid aspect ratio extend beyond a single adjacent cell in the vertical. A local calculation of neutral slopes steeper than the grid aspect ratio makes use of an extrapolation, which is less accurate than interpolation.

VENM interpolates as part of estimating neutral slopes (except for cases when the NTP outcrops into the surface or incrops into the bottom; Appendix A). Slope regularization is intrinsic to VENM as the maximum slope is limited by the ratio of the depth of the ocean (H) to the horizontal resolution of the data set. It replaces the ad hoc slope regularization required by the local method with a regularization based on using only information available from the discrete grid, without extrapolation. By searching over an arbitrary number of vertical grid cells in adjacent columns, VENM has more skill in estimating steep neutral slopes than the vertically local method (Figure 1). That is, where there could be a significant neutral tracer gradient over the grid scale, there may locally be a weak gradient.

Without some form of quantification of the effects of the different methods on ocean physics, it is not obvious that one approach is fundamentally better than the other. Although the local method is exact in principle, its finite‐difference realization suffers from numerical issues (extrapolation and regularization). While VENM could be considered “less fundamental” due to the vertical nonlocality in the calculation of the neutral directions, it suffers far less from numerical inaccuracies and arbitrariness, resulting in more neutral and physically sensible results.

As shown in this paper, the above technical details have important implications for the calculation of diagnostic analyses of observation‐based data. By extension, we conjecture that the same issues and sensitivities arise for parameterized tracer transport in prognostic ocean models. It is beyond our scope to pursue this conjecture in this paper.

Implementation of VENM

The Finite‐Difference Expressions

VENM's algorithm (section 3.2) calculates neutral slopes (S) and neutral gradients of S A, Θ, and pressure p. The algorithm is specifically designed for a z‐grid coordinate discretization of a tracer C and provides neutral gradients and slopes on predetermined locations between two vertical casts. The neutral slopes are calculated using the finite‐difference approximation of equation (11):

Here ensures that Δz is the change in height of the NTP over the horizontal distance Δx. The search algorithm determines the vertical distance Δz, whereas the spatial resolution provides the horizontal distances Δx and Δy. As a result, the maximum slopes that can be obtained are when Δz=H equals the depth of the ocean; that is, a NTP that extends from the bottom of one column to the top of an adjacent column. Smaller effective S max arises for points that are within the ocean interior. For a 1° spatial resolution and a 5,000‐m deep ocean, this provides S max≈0.05 as an example for the natural neutral slope regularization embedded in VENM.

In the algorithm we make use of the small‐slope approximation of the neutral gradients of equation (8):

Here , and are the unit normal vectors in the x, y, and z directions, respectively. Discretizing equation (14) renders

Again, indicates the tracer concentration difference in the x direction at constant y and z. Hence, is the gradient of C along the NTP in the x direction. Using equations (14), (15), (16), we rewrite ∇N C as

In section 3.2 we discuss how the algorithm provides , , , and .

The Search Algorithm

Here we present the algorithm to calculate neutral slopes ( and ) and neutral gradients ( and ), where C will later be replaced by S A, Θ, or p. This section is mostly conceptual, with more technical details provided in Appendix A. We emphasize that the algorithm as presented here is developed for gravitational stabilized vertical casts. We start by considering a north‐south transect (j+1 to j, where j is the latitudinal index), at constant longitude, with sufficient vertical resolution (Figure 1), and

We construct a finite length NTP between these two neighboring vertical casts, such that C and C represent the values of tracer C at the intersection of the NTP with the j (south) and j+1 (north) cast, respectively. Note that C and C are not required to be at the same height. A NTP is defined with respect to a local reference pressure p . We choose to define the NTP with respect to the midpressure , at the midlatitude between the casts. Consequently, the related neutral gradients are also defined at the (p m and y m). We apply two main steps: (1) iteration to find a neutral surface and (2) iteration to find the NTP that crosses through the predetermined grid point.

Step 1: Finding a NTP

In the first step of the algorithm we find the finite‐difference NTP. To do so, the intersection of the NTP with the south cast (j) is fixed at a T‐grid, that is, the locations where tracers are given (Figure 11 of Appendix A). We provide an initial guess of the height of the intersection of the NTP with the north cast (j+1). From both south and north intersections, we apply an adiabatic and isohaline displacement of fluid parcels (i.e., Θ and S A are constant, but p is not) to their midpressure p m at y m. The difference in the specific volume of seawater (υ=ρ −1) between the two displaced fluid parcels, with respect to their midpressure, is then given by

Figure 11

Description of a T‐grid box.

Here Δυ measures whether the parcels are on the same NTP with respect to p m. A larger value of |Δυ| indicates that both parcels are not on the same NTP (different ρ ). Hence, by finding an intersection for which Δυ=0, we can consider both intersections to be on the same ρ surface; that is, they are on the same NTP. Finding this intersection (where Δυ=0) is the key concept that defines VENM.

Note we deliberately chose to use υ instead of ρ. This choice is motivated by two reasons. First, we argue that it is more appropriate to define a quantity relative to a conserved property rather than to a nonconserved property. Mass is conserved by fluid parcels whereas volume is not. Thus, the more suitable property to consider is υ rather than its reciprocal, ρ. (For the full range of oceanic S A, Θ, and p combinations, Graham & McDougall, 2013, showed that turbulent mixing always destroys υ but does not always produce ρ.) Second, the McDougall and Barker (2011) and Roquet et al. (2015) software calculates υ (S A, Θ, and p) using a simple polynomial, making it is computationally more efficient to find the partial derivatives of υ than to find the corresponding partial derivatives of ρ.

The algorithm finds an approximate NTP over a finite distance Δy, between two vertical casts. As a NTP is only locally defined (equation (2)), we do not expect to find precisely Δυ=0. Rather, the NTP is generally determined within a predefined accuracy. Given an accuracy υ crit, an algorithm searches iteratively for the intersection on the north cast (with the south‐cast intersection fixed) that satisfies

Based on the following arguments, we judge this accuracy to be more than sufficient for our purposes. For a hydrostatic fluid, the buoyancy frequency is given by . Combining with Δρ≈−ρ 2 Δυ and Δυ=υ crit leads to an accuracy |Δp|≈g 2 ρ 2 υ crit/N 2≈0.1 dbar for a weakly stratified region with N 2=1×10−7 s−2. This result can be interpreted as a mismatch between the two NTPs at the midpoint pressure p m of at most 0.10 m in regions of weak stratification. Increased accuracy requires more computational iterations and thus more cost, with an accuracy of 0.10 m sufficient for our purposes. Once a NTP is constructed, we can use the values of the depth z of both end points to calculate the slope, while we use the values of the tracers (S A, Θ, and p) to calculate the gradients.

Step 2: Finding the NTP at a Target Pressure

In the previous step, we found an approximate NTP connecting a T‐grid point on the south cast to a particular location on the north cast. In other words, p on the south cast was fixed, while the pressure p on the north cast varied. As a result, p , where the neutral slopes and gradients are defined at, is most likely not on the same pressure as the T‐grid itself.

For the calculations of budgets, and tracer transports into and out of the T‐grid volumes, the neutral slopes and gradients are required to be defined at the middle of the interfaces surrounding the T‐grid (Griffies, 2012; Madec, 2015). This midpoint on a vertical interface is at the same pressure as the T‐grid but at a different longitude or latitude. The midpoint on a horizontal interface is at the same longitude and latitude but at the average pressure between two vertically spaced T‐grid points (Figure 11 of Appendix A). The goal of step 2 is to obtain neutral slopes and gradients that are defined at the specified pressure.

Consider the case of a neutral slope and tracer gradient in the y direction. We describe two options for how to proceed. The first option constructs NTP's for each T‐grid point (numbers will depend on the vertical resolution) and calculates the related neutral slopes and gradients on the interfaces. The neutral slopes and gradients are then interpolated exactly onto the predetermined pressures that correspond to the same pressure as the T‐grid. The second option uses an iterative process in which we also move the starting location vertically along the south cast, in order to search for a NTP that crosses through the predetermined pressures on the interface.

Because NTPs are defined locally, we adopt the second option, as this provides a more local and thus more accurate method. We iterate toward an accuracy of 0.50 m, meaning that NTPs have to cross the predetermined pressures within 0.50‐m distance. The numerical details of this iterative process are provided in Appendix A. After the algorithm has found the intersection of the NTP with the two casts, for the given criteria, the z, S A, Θ and p at the intersections are determined to provide and . The exact same method is applied for the x direction to obtain and . From the combination of the tracer gradients and slopes in the x and y direction and some horizontal averaging, we can calculate the vertical component of the tracer gradient (see equation (18)) on the horizontal cell faces.

Comparison of the Different Methods

To emphasize the importance of properly estimating neutral slopes and gradients, we here present diagnostics to compare the five different methods. We compare neutral slopes and Θ‐gradients obtained from the VENM, hybrid, and local methods as described above, and using σ 2 and γ n surfaces as approximate neutral surfaces (Appendix B). Also, the neutrality condition represented as a fictitious diffusivity, watermass transformation rates, and ocean lateral and vertical eddy heat transports are compared.

Data

We use monthly means of the World Ocean Atlas 2013, which is a set of objectively analyzed (1° grid) climatological fields of in‐situ temperature, practical salinity, and other tracers at standard depth levels for the World Ocean (Boyer et al., 2013). TEOS‐10 software (IOC et al., 2010; McDougall & Barker, 2011) is applied to convert the data to S A and Θ, achieve static stability (Barker & McDougall, 2017), and obtain Neutral Density (not available north of 64° N, Jackett & McDougall, 1997). For this data set we calculate slopes and gradients using the three methods described above. We include neutral slopes and gradients along σ 2 surfaces, using a method based on VENM (see Appendix A). We have excluded the Arctic region from this analysis, as the data in that area does not allow for the calculation of physically realistic gradients with any method. For sections (4.4)–(4.6), we make use of mesoscale eddy diffusivity, K, that parameterizes stirring by mesoscale eddies along the NTP. We use K = 1,000 m2/s for illustrative purposes, recognizing that there are a number of theories and observations suggesting that the diffusivity is flow dependent (e.g., Abernathey & Marshall, 2013; Cole et al., 2015; Groeskamp et al., 2017; Klocker & Abernathey, 2013; Roach et al., 2018). It is beyond our scope to investigate sensitivity to flow dependent diffusivities, with our study concerned instead with sensitivity to the calculation of the neutral slope and tracer gradients.

Comparing the Neutral Slopes

The most fundamental variables to compare between the local method and VENM are the neutral slopes in the x and y directions. A slope limit of S max=0.01 is commonly used in climate models and will thus be used for the local method (Cox, 1987; Danabasoglu & McWilliams, 1995, Griffies, 1998, 2004; Lemarié et al., 2012). We note that the maximum slopes found with the local method reach up to values of 200, that is, vertical. Applying the slope capping of to the local method affects roughly 2.8% and 1.8% of the slopes in x and y directions, respectively. For σ 2 and γ n, the slopes in the x direction were capped to 0.02. This limit was based on the maximum slopes found in the y direction, which were not capped. For VENM and the hybrid method, no regularization is applied to either the slopes or the tracer gradients.

Steep neutral slope regions often correspond to strong vertical transport and ventilation. These regions are thus important for the calculation of tracer gradients when using equation (5). Physical processes in such regions, for example, diffusive fluxes, will be particularly sensitive to the neutral slope.

A map and transect of the neutral slopes in January (Figures 2 and 3) reveal that the slopes resulting from the local method, VENM, and σ 2 have very similar spatial patterns. Results for γ n are very similar to those obtained by VENM (Appendix B). The slopes vary in magnitude more irregularly and on smaller spatial scales in the x direction than in the y direction. However, the slopes in the local method are notably spikier compared to the smooth patterns from VENM and from σ 2. Clear examples of such behavior can be found by comparing the slopes in the Southern Ocean both horizontally and at different depths (South of 40° S; Figures 2 and 3).

Figure 2

A map (depth = 300 m) and transect (longitude = 215.5) of S for January. The dashed‐blue line indicates the location of the transect (in the map) and the depth for which the map is plotted (in the transect). Positive slopes indicate neutral surfaces sloping upward to the north (for S ) or east (for S ). Black contours indicate Neutral Density γ n contours. VENM = vertically nonlocal method.

Figure 3

As Figure 2, but for S . VENM = vertically nonlocal method.

Comparing the Neutral Gradients

VENM, the hybrid method, local method, σ 2 and γ n provide five estimates of neutral tracer gradients. We discuss the neutral Θ gradients (∇NΘ) in x, y, and z direction (equation (18)), for the different methods (Figures 4, 5, and 6). Results for γ n are very similar to those obtained by VENM (Figure 12) The gradients for S A have very similar patterns and characteristics, while the gradients for P are very similar to the slopes, but with the opposite sign (not shown). The general patterns for the gradients for all methods are very similar. Gradients in the x direction (Figure 4) are more irregular than those in the y direction (Figure 5). The component of the neutral gradients in both x and y directions are about 3 orders of magnitude larger than the component of the gradient in the z direction (Figure 6).

Figure 4

A map (depth = 300 m, top) and transect (longitude = 215.5, bottom) of the neutral Θ gradients in the x direction , for January. The transect is limited to 2‐km depth, in order to show the more interesting parts. The dashed‐blue line indicates the location of the transect (in the map) and the depth for which the map is plotted (in the transect). Positive Θ gradients increase from west to east. Black contours indicate Neutral Density γ n contours. VENM = vertically nonlocal method.

Figure 5

As Figure 4, but for . When Θ gradients are positive, this indicates and increase of Θ when moving along the surface from along the surface south to north. VENM = vertically nonlocal method.

Figure 6

As Figure 4, but for . When Θ gradients are positive, this indicates and increase of Θ when moving along the surface from the shallow to the deep part of the neutral tangent plane. VENM = vertically nonlocal method.

Figure 12

A map (depth = 300 m, left column) and transect (longitude = 215.5, right column) of the Θ gradients along neutral density surfaces, in the x direction (top row), in the y direction (middle row), and in the z direction (bottom row), for January. The transect is limited to 2‐km depth, in order to show the more interesting parts. The dashed‐blue line indicates the location of the transect (in the map) and the depth for which the map is plotted (in the transect). Positive Θ gradients increase from west to east. Black contours indicate neutral density γ n contours. This figure is comparable to Figures 4, 5, 6.

It is important to consider the structure of the Θ gradients in x and y directions, south of 50° S (transect), between the surface and about 600‐m depth. These are some of the largest neutral gradients found in the ocean, due to the steep neutral slopes and strong temperature changes across the fronts of the Antarctic Circumpolar Current (Figures 4 and 5). VENM leads to smoother and arguably more realistic and accurate gradients than any of the other methods. As the neutral Θ gradient in the z direction is a linear combination of the gradients in the x and y directions, inaccuracy in the x and y directions translates into inaccuracies in the z direction. These inaccuracies are problematic for the parameterized vertical eddy heat fluxes (section 4.6).

Fictitious Diffusion

We make use of the neutrality condition as a measure of the skill for a method in estimating neutral slopes and neutral gradients. An exact NTP satisfies the balance between temperature and salinity gradients

Here L has units of inverse length, and ρ L can be interpreted as the change in locally referenced potential density, per meter, arising from nonneutrality. Any discrete approximation of the NTPs is inexact, thus leading to L≠0.

To help gauge the severity of nonneutrality, we convert L≠0 into a fictitious diffusivity D f as defined by Klocker et al. (2009)

The term ∇N z−∇a z measures the difference between the slope of the exact NTP (∇N z) and the approximate NTP (∇a z). When the approximate NTP is not exactly neutral, then diffusion along the approximate NTP leads to a dianeutral transport through the exact NTP, with that spurious diffusion measured by the fictitious diffusivity D f. To calculate D f, we approximate the slope difference ∇N z−∇a z, as the height difference along a NTP (Δz N) minus the height difference of the approximate NTP (Δz a). Using a finite horizontal distance , this height difference between the NTP and the approximate NTP is given by

Now we link this difference to L by noting that implying that there is a change in ρ along the approximate NTP as given by

Here Δρ can be interpreted as the density difference over the finite‐difference NTP. Noting that N 2≈−(g/ρ)(Δρ/Δz), this Δρ can then be related to the “vertical distance” Δz that one needs to move in order to find such a ρ change. Because Δρ is the density difference between the approximate and exact NTP, the related Δz can be interpreted as an approximate height difference between the approximate and exact NTP, which is the same as Δz . We are thus led to where we used equations (23), (24), (25), (26) to obtain the second part. By replacing L nloc with one obtained using gradients from the other methods, we can calculate , , , and .

In the ocean, a small‐scale diffusivity of (10−5 m2/s) is considered to be “background” mixing (Waterhouse et al., 2014). A D f that exceeds the background mixing adds potentially significant unphysical ocean mixing to calculations. The resulting distribution for D f (Figure 7) shows that for VENM 3% of the grid points have a value for which D f>10−5 m2/s. In contrast, the hybrid method has 22%, the local method 22%, and the σ 2 method 23 %, each revealing that about one in every five T‐grid points exceed D f=10−5m2/s. For Neutral Density about 13% of the grid point exceeds D f=10−5 m2/s, which is larger than VENM but less than the other methods. For VENM these few locations are in hot spots in the Southern Ocean, while the large D f values are more widespread for the other methods (Figures 8 and 13 for γ n).

Figure 7

The distribution of the values of the fictitious diffusivity ) for VENM (blue), the hybrid‐method (orange), local‐method (yellow), σ 2 (purple) and γ n (green). The thick vertical black line indicates values for . The percentage of grid points that have a D f>10−6 (rather than 10−5 as presented in the text) is (5%, 32%, 36%, 39% and 24%) for the VENM, Hybrid, Local, σ 2 and γ n methods, respectively.

Figure 8

As figure 4, but for D (m2s−1) for the month January.

Figure 13

The slopes of neutral density surfaces in the x direction (top row, comparable to Figure 2) and in the y direction (middle row, comparable to Figure 3), for January. The bottom row shows D m2/s) for the month January, for neutral density (comparable to Figure 8)

VENM satisfies the convergence constraint, Δυ ≤ υ crit (equation (21)), and so only approximates the neutral condition of α∇NΘ=β∇N S A (equation (22)). Except for the limits of numerical precision, and the limits set by step 2 of the search algorithm (section 3  3.2  3.2.2), one could argue that should be zero. However, as shown in Appendix B (see also Jackett & McDougall, 1997), satisfying Δυ ≤ υ crit is only approximately the same as satisfying α∇NΘ=β∇N S A. In this study, we assessed the neutrality of by means of α∇NΘ=β∇N S A, while the algorithm determines this neutrality by means of satisfying Δυ ≤ υ crit. As this is not exactly the same, this leads to nonzero values of that are unrelated to the accuracy of the algorithm but related to the method of assessment. may therefore be more accurate than shown in Figure 7.

Cabbeling and Thermobaricity

A globally integrated quantitative measure of the effect of the differences in neutral gradients on ocean physics can be obtained by calculating the water mass transformation WMT (Groeskamp et al., 2019; Walin, 1982) due to cabbeling and thermobaricity. We choose the use of cabbeling and thermobaricity, because it has important consequences for the formation of certain key ocean water masses (Evans et al., 2018; Groeskamp et al., 2016, 2017, Iudicone et al., 2008; Marsh, 2000; Nycander et al., 2015; Talley & Yun, 2001; Urakawa & Hasumi, 2012). However, its quantification depends strongly on an accurate calculation of neutral tracer gradients. Cabbeling and thermobaricity arise when Θ and S A mix along neutral directions in the presence of a nonlinear equation of state. Density of the resulting mixture is a nonlinear function of Θ, S A, and P and is generally not the same as the original density (Foster, 1972; McDougall, 1984, 1987b; Witte, 1902). Neutral diffusion along curved density isolines leads to an associated dianeutral transport that can be quantified in terms of WMT (Klocker & McDougall, 2010). Here the WMT due to cabbeling (T cab) and thermobaricity (T thb) in neutral density coordinates is given by (Groeskamp et al., 2016; Iudicone et al., 2011):

Here is the integral over the volume for which , such that is the WMT of γ n into a larger density. Positive (negative) values indicate transformation into denser (lighter) water (Figure 9). The dimensionless ratio b=|∇γ|/|∇ρ | results from the definition of γ n (Jackett & McDougall, 1997). Here C b and T b are the cabbeling and thermobaricity coefficients as defined by (McDougall, 1984, 1987a). The WMT rates (in Sv, where 1 Sv = 106 m3/s) are calculated and binned using monthly means. By replacing ∇NΘ and ∇N S A in equations (28) and (29) with those obtained by the different methods, their respective WMT rates are obtained ( , , , , and , and the same for T Thb).

Figure 9

Water Mass Transformation due to Cabbeling (upper panel) and Thermobaricity (lower panel), using neutral gradients from VENM (blue), the hybrid‐method (orange), local‐method (yellow), σ 2 (purple) and γ n (green). Positive (negative) values indicate transformation into denser (lighter) water. We note that WMT for σ 2 is calculated in σ 2 coordinates. The σ 2 axis is subsequently scaled to γ n coordinates using a polynomial fit between σ 2 and γ n.

The WMT due to cabbeling using neutral gradients from the local method or σ 2 leads to transport of 131 and 139 Sv, respectively. This is of the same magnitude as estimates of the Antarctic Circumpolar Current of 134–174 Sv (Donohue et al., 2016; Whitworth & Peterson, 1985). Such amounts of cabbeling would produce large amounts of water of Neutral Density γ n≈28.25 kg/m3, which is not observed. There is also no known mixing (Groeskamp et al., 2017), air‐sea flux (Groeskamp & Iudicone, 2018), or geothermal heating effect that compensates this production (de Lavergne et al., 2015).

Using an isopycnal model, Marsh (2000) finds annual mean WMT with maximum values of about 10 Sv for the Southern Ocean. Klocker and McDougall (2010) estimated about 5 Sv of cabbeling and 1 Sv of Thermobaricity for an annual mean climatology (WOCE; Gouretski & Koltermann, 2004). VENM, the hybrid method and using γ n have a cabbeling maxima of 22, 35, and 38 Sv, respectively. These results are likely higher because Marsh (2000) is not global; the study of Klocker and McDougall (2010) is performed on annual means, while all studies use different values for K. From these estimates and Figure 9, we conclude that the local method produces a physically unrealistic amount of cabbeling and thermobaricity, while VENM, the hybrid method, and using γ n produce results that are physically realistic. Note that both the cabbeling and thermobaricity estimates from VENM and the hybrid method show significant differences, with even a sign difference for thermobaricity for some densities.

Comparing Vertical Heat Transport

The last measure we employ to illustrate and quantify the differences between the methods is the heat flux caused by neutral diffusion. We can project the neutral flux into the lateral and vertical directions and integrate spatially to obtain the corresponding meridional and vertical eddy heat transports: with the heat capacity of seawater (J·kg−1·K−1). Here h(y,z) is the heat transport per meter latitude (W/m), such that H(z) is the integrated global vertical heat transport (W). The respective heat transport estimates are obtained by substituting ∇NΘ in equation (30) with neutral gradients obtained by the different methods.

The peak meridional heat transports change 25% using the local method compared to VENM (at 50° S and 40° N; Figure 10a), while the local method sometimes has the opposite sign to VENM (at 40° S). In addition, we find strong differences in the magnitude and sign of the estimates of h(y,z=1,000) (Figure 10b), which holds when globally integrated and for many different depths (Figure 10c). Figure 7 of Gregory (2000) shows the vertical heat transport of the HADCM2 atmosphere‐ocean general circulation model and has a comparable structure to VENM, the hybrid method and using γ n (Figure 10c). That is, a subsurface peak with positive values in about the first 1,500 m of the column and small negative values below that. Their maximum value is up to 4 W/m2, whereas that of VENM is about 1.5 W/m2.

Figure 10

Eddy heat transport H(y) (a), h(y,z=1,000) (b), and H(z) (c) using neutral gradients from vertically nonlocal method (VENM; blue), hybrid method (orange), the local method (yellow), σ 2 (purple), and γ n (green).

The local method produces a stronger midthermocline downward heat flux than the other three methods. The magnitude and direction of the inferred vertical heat transport depends on the method used to calculate neutral gradients; thus, it is relevant for large‐scale ocean circulation and climate.

Conclusion

In this paper we present VENM, which is a vertically nonlocal method to calculate neutral slopes and tracer gradients applicable to gridded data sets and numerical ocean simulations. The key concept defining VENM is to find a plane that intersects two vertical casts such that the difference in specific volume anomaly between these intersections is (near) zero (equations (20) and (21)). We compared VENM (based on a search algorithm) with local methods (based on a local grid stencil) that are currently implemented in many modeling and data studies. To help examine differences between VENM and the local method, we also introduced the hybrid method that uses the neutral slopes from VENM, combined with Cartesian tracer gradients based on the local grid stencil to calculate neutral tracer gradients.

We used a gridded climatology to calculate neutral slopes and gradients using the different methods. This approach corresponds to that used in observational analyses, such as water mass transformation diagnostics. It also ensures that the water mass structures are realistic. We furthermore made use of a neutrality condition to measure the accuracy of the neutral tracer gradient calculation. Our analysis reveals that VENM outperforms all other methods. Compared to the local method, VENM reduces fictitious diffusivity and produces far more realistic estimates of vertical heat transports and water mass transformation rates. The hybrid method suffers from numerical issues arising from local extremes that are produced when calculating Cartesian tracer gradients using the local grid stencil. That is, even when accurate slopes are somehow available, the fact that tracer gradients are obtained at the local grid stencil will still lead to inaccuracies. As such, this study quantifies a nontrivial sensitivity to calculating neutral slopes and gradients that is fundamental to numerical analyses of data.

Calculation of neutral slopes using the local approach is inaccurate for slopes steeper than the grid aspect ratio Δz/Δx, with the inaccuracies arising from the use of extrapolation for such slopes. The local method in turn suffers from numerical issues that lead to unphysical amounts of energy at the grid scale (reflected by spikes and noise) that are not observed in the ocean. To alleviate problems with the local method requires the use of ad hoc regularization methods that strongly influence resulting calculations. VENM improves on the local method by means of calculating NTPs using a vertically nonlocal searching algorithm, which can span the full depth grid aspect ratio H/Δx. As a result, VENM effectively applies a more accurate interpolation scheme rather than the extrapolation employed by local methods. With numerical models and gridded products aiming to increase their resolution, we expect that neutral slopes will often exceed the grid aspect ratio, thus further exposing the numerical problems with the local method. In contrast, VENM does not suffer under those conditions and will continue to perform well with refined spatial resolution. We do note that VENM does not deal with unstable vertical profiles (density inversions), in which case a different solution will need to be found.

It is not possible to satisfy the local down gradient form of discrete rotated diffusion using a linear, consistent, and local scheme (Beckers et al., 1998). Instead one needs to use the two‐dimensional projected nonorthogonal neutral tracer gradient ∇n C (Bleck, 1978a; 1978b; Starr, 1945). Although VENM is vertically nonlocal and obtains neutral slopes and gradients directly from the algorithm (equation (18)), the underlying physics does rely on the rotated diffusion tensor. As such VENM cannot guarantee downgradient diffusion over each grid cell and could produce tracer extremes. Although we show that VENM would much improve existing schemes, for numerical modeling purposes it may be more suitable to construct an algorithm that uses the vertically nonlocal nature of the VENM approach combined with a NTP‐based convergence scheme to always obey downgradient diffusion.

Isopycnal layered models typically make use of potential density referenced to 2,000 dbar, σ 2, to estimate neutral surfaces (e.g., Dunne et al., 2012). When using gradients of S A, Θ, and p calculated along σ 2 surfaces based on a searching algorithm, the results for water mass transformation, heat transport, and fictitious diffusivities are comparable in magnitude to that of the local method. This suggests that layered ocean models need a more accurate representation of neutral directions than the σ 2 surfaces typically used. Many data‐based studies use Neutral Density γ n to estimate their neutral surface. We here show that, although γ n does a reasonable job, results remain significantly less accurate than those from VENM.

VENM may influence other variables than those presented in this study, including parameterized eddy‐induced advective transport as per the Gent et al. (1995) scheme that is commonly used in mesoscale inactive climate models. Pradal and Gnanadesikan (2014) showed that changing the tracer fluxes (−K∇N C) by changing the mesoscale diffusivity K has a large impact on the climate. Here, instead of changing K, the flux is improved by having a more accurate representation of the neutral gradients ∇N C. Based on the nontrivial sensitivities found here, we conjecture that improving the calculation for neutral gradients will have a significant impact on numerical climate simulations.

VENM is computationally more expensive than the more widely implemented local method due to the use of iterations to find the NTP. However, efforts are underway to implement a VENM‐like method for neutral diffusion, using a computationally efficient method suitable for ocean modeling and allowing for density inversions. This model implementation will facilitate testing prognostic simulations using a neutral diffusion scheme that follows the vertically nonlocal approach proposed here.