This study provides observation-based estimates, determined by inverse methods, of horizontal and isopycnal eddy diffusion coefficients KH and KI, respectively, the small-scale mixing coefficient D, and the diathermohaline streamfunction Ψ. The inverse solution of Ψ represents the ocean circulation in Absolute Salinity SA and Conservative Temperature Θ coordinates. The authors suggest that the observation-based estimate of Ψ will be useful for comparison with equivalent diagnostics from numerical climate models. The estimates of KH and KI represent horizontal eddy mixing in the mixed layer and isopycnal eddy mixing in the ocean interior, respectively. This study finds that the solution for D and KH are comparable to existing estimates. The solution for KI is one of the first observation-based global and full-depth constrained estimates of isopycnal mixing and indicates that KI is an order of magnitude smaller than KH. This suggests that there is a large vertical variation in the eddy mixing coefficient, which is generally not included in ocean models. With ocean models being very sensitive to the choice of isopycnal mixing, this result suggests that further investigation of the spatial structure of isopycnal eddy mixing from observations is required.
A key role of the ocean in the climate system is to store and redistribute heat, freshwater, and biogeochemical tracers such as carbon dioxide. The manner in which this occurs depends strongly on the strength of mixing processes. Although ocean models are sensitive to the strength and the spatial distribution of the parameterized mixing, our knowledge and understanding of mixing processes from observations is limited.
In the ocean interior, mesoscale eddies stir properties along neutral tangent planes (Iselin 1939). Here, a neutral tangent plane is the direction along which a water parcel can move a small distance without experiencing a vertical buoyancy force (McDougall 1987a). As a result, ocean mixing is parameterized as 1) isopycnal downgradient tracer diffusion due to stirring by mesoscale eddies by means of an isopycnal eddy diffusion coefficient KI and 2) small-scale isotropic downgradient turbulent diffusion by means of a turbulent diffusion coefficient D (Redi 1982; Griffies 2004; McDougall et al. 2014). Here, mesoscale refers to length scales of the order of 50–300 km, while small scale refers to processes such as breaking internal waves and shear-induced turbulence and is of the order of millimeters to tens of meters. Mesoscale mixing in the mixed layer is alternatively defined as surface diabatic mixing (Tandon and Garrett 1997) and is parameterized as horizontal downgradient tracer diffusion due to stirring by mesoscale eddies by means of a horizontal eddy diffusion coefficient KH. The isotropic nature of D is discussed in McDougall et al. (2014) but has previously been regarded to be dianeutral (Redi 1982; Griffies 2004) or vertical in the small-slope approximation.
a. Small-scale mixing
Small-scale mixing leads to a downward flux of heat that is required to warm and subsequently upwell abyssal dense water (Ganachaud and Wunsch 2000; Sloyan and Rintoul 2000) and plays a crucial role in closing the ocean circulation (Munk 1966; Munk and Wunsch 1998; Wunsch and Ferrari 2004). Sources of small-scale mixing are (for example) wind-generated mixing by breaking of near-inertial waves near the surface and in the ocean interior (D’Asaro 1985; Alford 2001; Alford et al. 2012), dissipation of the internal tides (Nycander 2005), and generation and dissipation of lee waves (Nikurashin and Ferrari 2011). Estimates of small-scale mixing have been obtained from measurements using microstructure instruments (Polzin et al. 1997; Garabato et al. 2004), tracer release experiments (Ledwell et al. 1993), inferred from lowered ADCP and CTD profiles (Kunze et al. 2006), or calculated using moored profilers (Alford et al. 2006) or floats (Whalen et al. 2012; Meyer et al. 2015).
Microstructure observations have indicated that, in the ocean interior, D increases toward the ocean floor with maximum values exceeding D = 1 × 10−3 m2 s−1 in regions of rough topography. The combination of all small-scale mixing processes results in a highly spatially variable structure for D (Waterhouse et al. 2014), with values ranging between O(10−7) and O(10−3) m2 s−1. Strong sensitivity of ocean and climate models to differences in the magnitude of the parameterized vertical mixing argues that more studies are required to further constrain these estimates (Simmons 2004; Melet et al. 2013).
b. Mesoscale mixing
Eddy transports dominate dispersion of particles and mixing of tracers on large spatial and temporal time scales. In numerical models the eddy diffusion coefficients KH and KI parameterize tracer diffusion due to eddies that are not resolved. As a rule of thumb, the higher the resolution of a model, the more eddy transport is resolved. Therefore, higher-resolution models tend to need lower diffusion coefficient (Okubo 1971). In this study, we are considering the eddy diffusivity to parameterize tracer diffusion (Redi 1982), which is not the same as the quasi-Stokes streamfunction parameterization for eddy-induced advection [or Gent–McWilliams (GM) parameterization; Gent and McWilliams 1990; Gent et al. 1995; McDougall and McIntosh 2001].
Estimates of mesoscale diffusivities have been based on tracer release experiments (Ledwell et al. 1993, 1998), float dispersion (Zhurbas and Oh 2004; LaCasce et al. 2014), satellite data (Holloway 1986; Abernathey and Marshall 2013; Klocker and Abernathey 2014), or, most recently, a combination of Argo and the ECCO2 state estimate (Cole et al. 2015). The observations clearly illustrate a highly spatially variable pattern with enhanced mixing in western boundary currents and a general decay of the mesoscale diffusivities with depth, with the exception of some subsurface maximums and values ranging between 0 and O(104) m2 s−1.
Allowing for spatial variation of eddies in climate models can reduce systematic drift (Ferreira et al. 2005; Danabasoglu and Marshall 2007). Differences in the KI value used in numerical climate models can lead to differences in the simulated global surface temperature change by the end of the twenty-first century of up to 1°C (Pradal and Gnanadesikan 2014), influence the stability of the overturning circulation (Sijp et al. 2006), and affect processes due to the nonlinear equation of state, such as cabbeling and thermobaricity (McDougall 1987b; Iudicone et al. 2008; Groeskamp et al. 2016). It is therefore essential that coarse-resolution climate models use the most accurate parameterization of the spatial varying effect of eddies.
c. Estimating mixing from an ocean climatology
In this study, we will provide a global, observation-based quantification of circulation and mixing. This is obtained by applying an inverse method to an ocean climatology, which is based on spatial and temporal averaging of observations of the ocean’s hydrography. This is possible because the spatial and temporal structure of the ocean’s hydrography is set by the processes acting upon it (e.g., tides, winds, surface freshwater and heat fluxes, and mixing). As a result the ocean’s hydrography holds information about these processes and the circulation resulting from it. The inverse model is used to extract and estimate the ocean dynamics that are embedded in a climatology.
Historically box inverse models have been used to estimate the absolute velocity from ocean hydrography (Wunsch 1978). With increasing spatial and temporal resolution of available ocean climatologies, it is now possible to obtain estimates of more complicated (higher order) processes, such as ocean mixing. However, this requires inverse methods that are specifically designed for that purpose. Some box inverse methods have estimated D but often contain unknowns at the boundaries that require both dynamical constraints and conservation statements, increasing the complexity of the system and sensitivity to prior estimates, leading to uncertainties in the obtained transport rates and especially D (Zika et al. 2010a).
2. The thermohaline inverse method
This section provides a conceptual explanation of the diathermohaline streamfunction (Groeskamp et al. 2014a) and subsequently the derivation of the thermohaline inverse method (THIM; Groeskamp et al. 2014b). For detailed mathematical derivations and assessment of both the diathermohaline streamfunction and the THIM, we refer the reader to Groeskamp et al. (2014a,b) and Hieronymus et al. (2014). The THIM estimates mixing and circulation in Absolute Salinity SA and Conservative Temperature Θ coordinates. Here, Conservative Temperature is proportional to potential enthalpy (by the constant heat capacity factor , in J kg−1 K−1), representing the heat content per unit mass of seawater (McDougall 2003; Graham and McDougall 2013). Absolute Salinity is measured on the Reference Composition Salinity Scale (Millero et al. 2008) and is an approximation to the mass fraction of dissolved material in seawater (g kg−1; McDougall et al. 2012) and is the salinity variable of the IOC et al. (2010) thermodynamic description of seawater.
a. The diathermohaline streamfunction
The Walin framework (Walin 1982) originally considered water mass transformation from one temperature to another. The Walin framework has also been applied to calculate water mass transformation as the change in a fluid parcel’s density (Marshall 1997; Marshall et al. 1999; Nurser et al. 1999; Marsh 2000; Iudicone et al. 2008; Badin and Williams 2010; Badin et al. 2013). Here, we define a water mass transformation as a change in the (SA, Θ) values of a water parcel (Speer and Tziperman 1992; Speer 1993, 1997). Water mass transformation is then directly related to heat and salt fluxes by mixing and air–sea interactions (i.e., thermohaline forcing), as these are the only mechanisms that can change a water parcel’s SA and Θ.
Groeskamp et al. (2014a) showed that nondivergent transformations (those that do not cause convergences or divergences of the volumes of water masses) can be represented by a streamfunction in (SA, Θ) coordinates known as the diathermohaline streamfunction. We will use the symbol Ψ for this streamfunction. In Groeskamp et al. (2014a) and Zika et al. (2012), a subscript SAΘ was used to distinguish it from streamfunctions in other coordinates. In Groeskamp et al. (2014a), a superscript “dia” was also used to distinguish the diathermohaline streamfunction that describes the total transformation across SA and Θ surfaces from those that describe the flow normal to SA and Θ surfaces (advective thermohaline streamfunction Ψadv; Döös et al. 2012; Zika et al. 2012) and the movement of those surfaces separately (local thermohaline streamfunction Ψloc). Using numerical models, Hieronymus et al. (2014) and Groeskamp et al. (2014a) provided a thorough analyses of the circulation cells of the diathermohaline streamfunction Ψ.
b. The thermohaline inverse method
Groeskamp et al. (2014b) developed the THIM to be applied to an ocean climatology and allow for observation-based estimates of diffusion coefficients and Ψ based on estimates of air–sea fluxes and a priori estimates of the spatial distribution of mixing. The THIM exploits the fact that circulation in (SA, Θ) coordinates requires a freshwater and/or heat flux convergence that can only be provided by boundary fluxes and/or ocean mixing.
Consider an incompressible ocean where the total volume of water over a particular temperature and salinity range is constant when integrated over some sufficiently long time period such that
Above we have used a conditional volume integral over the whole ocean, where SA is between and and Θ is between and . Equation (1) allows us to define a two-dimensional streamfunction Ψ, describing the flow in (SA, Θ) coordinates. Trends in water mass volumes can be accounted for in this method (Groeskamp et al. 2014b) but will not be considered here as we use an annual climatology of ocean temperature and salinity and surface fluxes (see section 3).
The transformation across an isohaline (surface of constant SA) between two isotherms (lines of constant Θ) can be represented by the change in Ψ with respect to Θ at constant SA such that (Fig. 1)
Here, is a mixing tensor (m2 s−1), and fs = SA(E − P − R), if the equivalent salt flux due to evaporation E, precipitation P, and runoff R, such that ∇fs is the net salt flux divergence (g kg−1 s−1) into the volume and LS is a correction term (m3 s−1), which accounts for the discretization in the SA coordinate [see Eq. (4)]. The transport across an isotherm between two isohalines is likewise (Fig. 1)
where , with Q as the boundary heat flux (W m−2) and ρ0 as a constant reference density (kg m−3), such that ∇fΘ is the rate of cooling of a water parcel due to the divergence of the heat flux (K s−1). The quantity LΘ is a correction term (m3 s−1), which accounts for the discretization in the Θ coordinate.
The LS and LΘ terms are referred to as the local response term in Groeskamp et al. (2014a), defined as
where and are volume-weighted averages within and . These terms account for the change in salt and heat content that causes a change in SA and Θ but does not cause a transformation across the limits of the integration. They are a correction for the discretization that vanishes in the continuous limit where ΔSA and ΔΘ tend to zero. For a continuous form of Eqs. (2) and (3), see Hieronymus et al. (2014). The LS and LΘ terms turn out to be small, as we use a relatively high resolution in (SA, Θ) coordinates in this study.
We do not account for the explicit transport of water through the atmosphere, as this is very small relative to the transformation due to air–sea fluxes (evaporation, precipitation, and runoff are effectively represented as the addition and removal of salt).
The mixing tensor is in reality a function of space and time and represents the anisotropy, for example, between mixing along and across isopycnal surfaces. To make the problem of estimating mixing coefficients tractable, we use prior estimates of the spatial distribution of three varieties of mixing: horizontal KH, isopycnal KI, and isotropic D. We represent each of these by tensors that vary in space (, , and , respectively) and solve for three constants that multiply those tensors (, , and , respectively), such that
For brevity we will use the symbols , , and to represent the terms, which multiply , , and Dinv, respectively, and FS to represent the contribution of surface fluxes, such that we can rewrite Eq. (6) as
For the purposes of our inverse method Ψ, , , and Dinv are unknown, while the A, F, and L terms can be estimated from observations. In section 4, we will describe the inversion that solves for the unknowns Ψ, , , and Dinv. In the following section 3 we describe the source of our estimates of the surface flux and SA and Θ distributions and mixing tensor structures that allow us to calculate the A, F, and L terms in Eqs. (7) and (8).
3. Data and estimated terms
a. Ocean climatology
The observation-based climatology for all latitudes north of 75°S is given by the CSIRO Atlas of Regional Seas 2009 (CARS; Ridgway et al. 2002; Ridgway and Dunn 2003). CARS provides in situ temperature T, practical salinity SP at a resolution of 0.5° × 0.5° grid spacing, and 79 vertical levels for each month of a standard year. Here, we use monthly averaged values and apply the International Thermodynamic Equation Of Seawater—2010 (TEOS-10) software (IOC et al. 2010; McDougall and Barker 2011) to convert the CARS data to SA and Θ. The resulting data are stabilized (such that density increases with depth) using a minimal adjustment of SA and Θ, within the measurement error (Jackett and McDougall 1995; Barker and McDougall 2016, manuscript submitted to J. Atmos. Oceanic Technol.).
b. Boundary forcing
The freshwater (E − P − R) and heat (Q) fluxes are taken from the product of Yeager and Large (2008) commonly known as CORE2 (although CORE2 specifically relates to the atmospheric state from which the fluxes are derived; Large and Yeager 2009). The Yeager and Large (2008) product has a spatial resolution of 1° × 1° grid spacing and is converted to a standard year by averaging surface fluxes for each calendar month and interpolating onto the CARS grid.
c. Horizontal and isopycnal mixing
We separate the effect of the spatial structure of mixing by mesoscale eddies into a tensor that operates on horizontal gradients within the mixed layer and one operating on isopycnal gradients below the mixed layer.
We base our spatial variation of the magnitude of eddy diffusion on the recent estimate of Cole et al. (2015), providing the scalar field KCole. Cole et al. (2015), from Argo profiles between the base of the winter mixed layer to approximately 2000 m, used salinity anomalies and gradients to resolve the three-dimensional structure of isopycnal mixing on a 3° × 3° grid between 60°S and 60°N.
Near the surface, the estimate from Cole et al. (2015) is comparable to estimates obtained from satellite data (Abernathey and Marshall 2013; Klocker and Abernathey 2014), so we also use Cole’s estimates to set the structure of horizontal mixing at the sea surface.
To remove spikes in these data, we impose a maximum diffusivity of 2.5 × 104 m2 s−1. The absolute values from Cole et al. (2015) reflect the eddy mixing coefficient required to parameterize unresolved eddy fluxes of a 3° × 3° resolution climatology. We aim to diagnose our own mixing value based on our 0.5° × 0.5° resolution climatology. So for the structure tensor we divide the Cole estimates by 2.5 × 104 m2 s−1 to yield a scalar function between 0 and 1 (Fig. 2b).
For a smooth transition between horizontal mixing near the surface and isopycnal mixing in the interior we apply a linear transition from the base of the mixed layer (at z = hm; z decreases downward from 0 at the sea surface) to 500 m beneath it. This is achieved using the tapering function β such that
Here, hm is calculated using TEOS-10 software based on de Boyer Montégut et al. (2004) (Fig. 2a). To include the entire pycnocline, the thickness of the transition layer is chosen to be 500 m. This represents a simplified version of the approaches of both Ferrari et al. (2008, 2010). Our structure tensors for horizontal and isopycnal mixing are then
where Sy and Sx are the slopes of locally referenced potential density surfaces in the zonal and meridional directions, respectively, and . Our formulation of uses the small-slope approximation as is typical in ocean modeling (Griffies 2012). Further details of the way isopycnal gradients are calculated are given in appendix A.
To extend Cole et al.’s (2015) diffusivity estimates to the seafloor we extrapolate vertically using
where hmax(x, y) is the deepest level at which a diffusivity estimate exists at a particular x and y location, and we have chosen a decay of 75% over 1500 m, based on the horizontal average of the original data. Mixing values KCole in the mixed layer and at high latitudes are calculated using a nearest neighbor extrapolation. These results are then interpolated onto the CARS grid (Fig. 2d).
d. Small-scale diffusion
Small-scale diffusion is parameterized by multiplying the small-scale diffusion coefficient D with isotropic tracer gradients. To include the effect of increased small-scale isotropic mixing near topography, we will use a structure tensor based on Bryan and Lewis (1979). Hence, is given by
Here, || = 0.3 at the surface, increasing to 1.3 for depths greater than 2500 m, mimicking the effect of increased, small-scale mixing near the bottom (Fig. 2c). Note that the distribution of small-scale diffusion used in numerical climate models by Bryan and Lewis (1979) was determined by multiplying by 10−4 m2 s−1. Based on an inverse calculation we will find our own estimate for the magnitude of this isotropic mixing.
e. Estimation of terms in Eqs. (7) and (8)
We use CARS and the mixing tensors to estimate the terms in square brackets in Eqs. (7) and (8). In CARS, SA and Θ for each month are given at each (x, y, z) point around which six interfaces are defined to enclose a volume. Based on the CORE2 forcing, the local sources of heat fΘ and freshwater fS are determined. These are then multiplied by the volume of each grid, and then all points within the defined SA and Θ interval are summed. Likewise for the L terms local changes are determined between months, and these are then averaged in time.
To calculate the diffusive flux divergences due to the mixing tensors, the gradient of SA and Θ is calculated at each of the six interfaces of each box. These are then multiplied by the area of the interface and the component of the diffusive flux normal to the interface based on the mixing tensors. The divergence is then the sum of these fluxes over all six interfaces. The volume integral is then the conditional sum of these divergences.
4. The inverse technique
Equations (7) and (8) describe the relationship between the transformation of water across isohalines and isotherms as a result of a convergence of salt and heat, respectively, defined at . For the inverse model, Eq. (1) dictates that the sum of the transports through each of the four interfaces of the grid in (SA, Θ) coordinates has to add up to zero. The inverse method therefore solves for the diathermohaline streamfunction Ψ defined at the corners of a grid in (SA, Θ) coordinates (Fig. 1, black colored Ψ and ΔΨ). Hence, the transports related to each ΔΨ are calculated using Eqs. (1)–(8) for volumes in (SA, Θ) coordinates simply shifted by either ΔSA or ΔΘ.
Now divide the Θ and SA range into and grids, respectively. This gives unknowns, where the 3 represents the unknown diffusion coefficients (, , and ) and and streamfunctions for the Θ and SA range, respectively. For each grid two unique equations can be constructed, leading to equations. From this we can combine Eqs. (7) and (8) into a single matrix expression , where x is a 1 × M vector that contains the unknowns; is a N × M matrix that contains 0, ±1, and the A terms; and b is a 1 × N vector that contains the F and L terms. For the analyses, we have chosen grid sizes in coordinates to be g kg−3 and °C for [31.8, 37.5] g kg−1 and Θ ∈ [−2, 31]°C, respectively. The resolution is high enough to distinguish different water masses in the ocean’s interior and provide approximately equal number of equations in both SA and Θ directions, such that both have an equal influence on the solution.
Appendix B details exactly how the inverse model is applied, including the different column and row weighting of the matrix and the a priori initial estimates applied, and leaves N = 34 582 and M = 20 600. In section 5, we discuss the results obtained from the inverse estimate based on an inversion with one particular choice of weighting parameters, as described in appendix B. The reported standard deviation associated with the results is calculated using Cp [Eq. (B2)]. It is shown that the inverse model provides a better solution then a priori estimates. Results of a sensitivity study to changes in the prior estimates and weighting are provided and discussed in appendix C. The solution described here is within the physically realistic limits and is well constrained by the data and is only moderately sensitive to different weighting of the matrix.
The surface freshwater and heat fluxes lead to water mass transformation rates in cubic meters per second per Θ bin, giving units of Sverdrups (1 Sv ≡ 106 m3 s−1) per kelvin, and Sverdrups per SA bin, giving Sverdrups per gram per kilogram, respectively. Water mass transformation due to salt diffusion by small-scale mixing processes in (SA, Θ) coordinates is given by . A similar calculation is done for the horizontal and isopycnal eddy diffusion and for the heat terms. Water mass transformation for surface and diffusive fluxes are blue (red) when directed toward lower (higher) salinities or temperatures (Fig. 3).
a. Surface forcing
The distribution of the surface freshwater fluxes show three areas of freshening and one area of salinification (Fig. 3a), generally acting to broaden the salinity distribution (Zika et al. 2015). For C, an area that we associate with the (sub) tropical mixed layer, less saline water is freshened even more and salty water increases in salinity even more. The freshening at lower temperatures (C) is also acting upon already relatively less saline water, associated with higher latitudes in the Southern Hemisphere ( g kg−1) and in the North Pacific ( g kg−1; Fig. 3a).
The distribution of the surface heat flux shows a general tendency for warming of already warm waters (C; Fig. 3e). There are sloping bands of heating and cooling aligned next to each other for g kg−1 and °C. This banding is a result of seasonal cycles at midlatitudes. When considering a constant , the heating is at a higher temperature than the cooling. As with freshwater forcing and salinity, water mass transformation due to surface heat fluxes acts to broaden the temperature distribution. This is also valid for the strong warming/cooling dipole at °C and g kg−1, related to heat fluxes at high latitudes.
b. Small-scale diffusion
As mixing acts to destroy tracer gradients, the small-scale mixing results in an opposite structure to that of surface forcing (Figs. 3b,f), especially for C. This places strong constraints on the inverse estimate of × 10−5 m2 s−1 (Table 1). Note that our inferred diffusivity is roughly half the value of 10−4 m2 s−1 used by Bryan and Lewis (1979) to multiply the vertical profile [Eq. (12)].
As a result of the Bryan–Lewis structure, we find that for the upper 2000 m × m2 s−1. For the depths greater than 2000 m, we find that D increases to × 10−5 m2 s−1. We used 2000 m to separate interior from upper ocean, as this is the depth to which small-scale mixing is required to upwell bottom waters before isopycnal upwelling becomes relevant, reducing the required mixing to close the circulation (Toggweiler and Samuels 1998; Sloyan and Rintoul 2001). The upper-ocean estimate corresponds very well with the observed values by Waterhouse et al. (2014) and the inverse estimates of Lumpkin and Speer (2007). The deep value corresponds to that of Lumpkin and Speer (2007); it is about half of canonical 1.3 × 10−4 of Munk (1966) and about an order of magnitude lower than that of Waterhouse et al. (2014) and Ganachaud and Wunsch (2000).
c. Horizontal eddy diffusion
The water mass transformation rates due to horizontal eddy diffusion in the mixed layer also generally oppose the surface forcing (Figs. 3c,g). We obtain m2 s−1, which multiplies the structure function that corresponds to a spatial resolution of 0.5° × 0.5°. Recall that Cole et al. (2015) found profiles with a maximum 3 times this value of 25 000 m2 s−1. Their estimates represented a lower resolution and are therefore expected to have larger values than those presented here.
d. Isopycnal eddy diffusion
By definition the isopycnal eddy diffusion only affects the interior circulation of the ocean, reducing its spread in (SA, Θ) coordinates significantly (Figs. 3d,h). We obtain m2 s−1, which multiplies the structure function that corresponds to a resolution of 0.5° × 0.5°, that is, 20 times smaller than that in the mixed layer. This will lead to small values in the ocean interior, as a result of the decay with depth, as observed by Cole et al. (2015). A local inverse estimate of Zika et al. (2010b) obtained small values for KI in the interior. Even though the Zika et al. (2010b) estimate is a local estimate, it indicates that KI may significantly decay with depth, suggesting that our result is plausible.
e. The streamfunctions
An estimate of the nondivergent circulation in (SA, Θ) coordinates is given by our solution for (Fig. 4a). The term represents net transformation; it is the sum of the Eulerian flow normal to surfaces of constant and minus the rate at which those surfaces move (if a property surface moved with the flow there would be no transformation). Groeskamp et al. (2014a) showed that these two effects could be represented by distinct streamfunctions and , respectively. The interpretation of is unchanged from that thoroughly explained in Groeskamp et al. (2014a), who used the same CARS data. Hence, we focus our attention on and .
1) Tropical cells
For C, an area in (SA, Θ) coordinates related to the upper 200 m of the ocean (Fig. 5) in the (sub) tropics, the surface freshwater and heat fluxes are balanced by small-scale and mesoscale mixing processes in the mixed layer. Although strong small-scale mixing near the surface is also supported by model-based studies (Ferrari and Ferreira 2011; Hieronymus et al. 2014), we find stronger near-surface mesoscale mixing than those studies. The residual water mass transformation results in a strong residual anticlockwise circulation (Fig. 4b), which may be considered to be the observation-based equivalent of the tropical cell, identified as a shallow wind-driven (near) equatorial cell (Döös et al. 2012).
Hieronymus et al. (2014) found that the tropical cell (estimated from a numerical ocean model) is composed of an Atlantic and an Indo-Pacific circulation. We also find a slight extension into the Atlantic; however, the main difference is that the extension of this cell is at similar temperatures but is centered at g kg−1. This cell exists only in the upper 100 m of the freshest area in the Indian Ocean (Fig. 5) and is most likely related to the circulation in the Bay of Bengal, where freshwater discharge due to the summer monsoon lead to the observed low salinity (Schott and McCreary 2001; Schott et al. 2002).
2) Global cell
The large clockwise-rotating (blue) circulation cell in (Fig. 4b) is a result of a balance between surface forcing, global upwelling through small-scale diffusion (Munk 1966; Gordon 1986; Broecker 1991), and both horizontal and isopycnal eddy diffusion leading to isopycnal upwelling (Toggweiler and Samuels 1998; Sloyan and Rintoul 2001). This cell is very similar in both and and can be considered to be mainly related to advection. This cell connects the different ocean basins with the Southern Ocean (Figs. 4b, 5), making it the observational-based equivalent of the global cell for an ocean model (Zika et al. 2012; Döös et al. 2012; Hieronymus et al. 2014).
3) Cold cells
Both and show two anticlockwise-rotating cells at low temperatures (Figs. 4a,b). For g kg−1 and °C, the related circulation is in the high-latitude Southern Ocean (hereafter referred to as the Antarctic cell). For g kg−1 and °C, the related circulation is in the Arctic and North Pacific (Fig. 5).
The Antarctic cell covers an area in (SA, Θ) coordinates that connects surface dynamics with the interior and is related with strong surface cooling and relatively small freshwater fluxes (Figs. 3a,e). Small-scale and mixed layer mesoscale mixing is observed on the kg m−3 contour at g kg−1 (Figs. 3b,f). Hence, this circulation cell represents the Antarctic Bottom Water (AABW) and Lower Circumpolar Deep Water (LCDW) circulation and can be identified as the observational-based equivalent of the AABW cell as found by Zika et al. (2012) and Döös et al. (2012) for an ocean model. However, the AABW we observe occupies a wider and range than the model-based versions.
The Arctic cell, related to the Arctic and the North Pacific (subpolar) Gyre (NPG) for depths shallower than 200 m (Figs. 4b, 5), is mainly induced by surface freshwater and heat forcing (Fig. 3) and perhaps affected by the Bering Strait flow, which is a northward flow of about 1 Sv (Roach et al. 1995). The North Pacific Gyre rotates clockwise, transporting cold and freshwater equatorward on the eastern side of the basin. Water moving to the equator and subsequently to the western side of the basin will warm and the salinity will increase (evaporation), as a result of surface heat fluxes. At the western boundary, the Kuroshio transports water northward, leading to cooling and freshening by net freshwater flux from precipitation. This NPG cell has not been identified in any of the model-based versions of (Zika et al. 2012; Döös et al. 2012) or (Groeskamp et al. 2014a; Hieronymus et al. 2014).
The cell observed at g kg−1 and Θ = [0, 10]°C is likely a manifestation of the Arctic circulation that may be related to ice dynamics in combination with strong seasonal warming and cooling, leading to a very shallow buoyancy-driven circulation. However, because of sparseness of data in these regions this remains speculative.
Regardless of the fact that we have not included any prior knowledge of the streamfunction structure ( appendix B) the THIM provides an estimate of comparable to that calculated in numerical models (Zika et al. 2012; Döös et al. 2012; Hieronymus et al. 2014), without the use of a priori advective transports.
In this section, we use the inverse estimate of and the diffusion coefficients to calculate energy consumption by small-scale mixing processes. We will then discuss the accuracy of the inverse estimate and the implications of the results.
a. Dissipation of turbulent kinetic energy
Using the Osborn (1980) relationship between the isotropic diffusion of small-scale mixing processes and energy dissipation in combination with our inverse estimate of , we calculate the power used for small-scale mixing using , where is the mixing efficiency, is the Buoyancy frequency squared (s−2), ρ is in situ density [kg m−3, calculated using IOC et al. (2010) software], and V is the volume (m3). These bulk estimates contain errors related to the climatology, choice of , and estimate of ; however, they still provide a worthwhile indication of the power needed to drive the circulation we have inferred.
The global, depth-integrated distribution of P is dominated by the high values of near the surface, resulting in maximum dissipation near coasts and high dissipation at low latitudes [Fig. 6 (left)]. The total globally integrated energy used by small-scale mixing is 1.34 TW, smaller than the 3.5 TW of energy available from tides (Munk and Wunsch 1998). Only 0.17 (13%), 0.31 (23%), and 0.45 TW (33%) are dissipated via small-scale mixing processes below 2000, 1000, and 500 m, respectively [Fig. 6 (right)]. Hence, 0.89 TW (67%) is dissipated in the upper 500 m, which is consistent with the idea that most of the mixing occurs near the surface (Wunsch and Ferrari 2004). The total energy conversion of 0.45 TW for depths greater than 500 m is of the same order as the available energy for mixing through the sum of the conversion of energy from the barotropic tide to internal tide [ TW; Nycander 2005] and that from the geostrophic flow into lee waves [ TW (Nikurashin and Ferrari 2011)].
The value of Dinv is primarily set as a result of compensating surface fluxes, which, through the use of the Bryan and Lewis (1979) structure function, leads to strong constraints on the interior estimate of D. We anticipate that the interior estimate of D can be improved using more complex parameterizations that separate the different physical mechanisms, such as surface mixing from internal wave breaking in the deep ocean (Melet et al. 2013; Nikurashin and Ferrari 2013). This could lead to more detailed estimates of the required energy for the different small-scale mixing processes and shed light on how small-scale mixing influences the global overturning circulation by transporting bottom water to the level where isopycnal upwelling plays a significant role (Toggweiler and Samuels 1998; Sloyan and Rintoul 2001; Ferrari 2014).
b. Uncertainty and improvements to the solution
Volume divergences in (SA, Θ) coordinates are due to the imbalance between surface and diffusive fluxes of (equivalent) salt and heat and would be zero in a steady-state ocean. The THIM simultaneously minimizes volume divergences for each grid over the global (SA, Θ) domain by adjusting the diffusion coefficients, while accounting for an exactly nondivergent . The THIM greatly reduces volume divergence as a result of the a priori estimates of the diffusion coefficients but does not exactly balance the volume divergences (as detailed in appendix B).
To connect with the literature while quantifying the imbalance, we will quantify water mass transformation using the surface-referenced potential density anomaly coordinates. Using the estimates of , , and , we calculate the convergence of the diffusive heat flux:
By replacing by , a similar expression can be found for the salt flux. This allows us to quantify the total diapycnal transformation due to mixing and surface forcing:
where we use a conditional volume integral over the whole ocean, and is between and . We have calculated the total transformation T and the transformation due to surface fluxes TF, small-scale mixing TD, horizontal mesoscale mixing TH, and isopycnal mixing TI (Fig. 7).
We want to emphasize that regardless of the models for mixing that we use, the THIM gives an exactly balanced solution for the total volume transport. This is implicit in solving for the streamfunction . However, there is a mismatch between the implied transformation rates due to each processes when the diffusive terms are quantified based on the estimated coefficients , , and . Inferred total transformation rates in a density framework (Fig. 7) show the integrated effect of this mismatch.
For a steady-state ocean, T would be zero (green line, Fig. 7), which mainly requires TF to be compensated by mixing (TD, TH, and TI). We suggest that the imbalance of T is mainly due to the lack of variability that we allow in D and K, resulting in imbalances within some isopycnal ranges. For example, the imbalance at kg m−3 is most likely because the implemented does not allow for enhanced near-surface mixing, while is obtained by minimizing volume divergences over the whole (SA, Θ) domain. This limits the THIM to balance TF correctly everywhere. That is, our mixing estimates are more tailored for diagnosing the transformations at the median ocean than high .
Between kg m−3, the total water mass transformation is dominated by surface fluxes TF, indicating transformation to lighter density water. This result agrees with that obtained by Iudicone et al. (2008), who used a numerical ocean model with a different surface forcing. The term TI is significant for kg m−3 and is positive and of the same order of magnitude as that of Klocker and McDougall (2010), who calculated TI using the WOCE climatology (Gouretski and Koltermann 2004) with a constant m2 s−1.
Groeskamp et al. (2016) estimate water mass transformation by directly applying the estimates of Cole et al. (2015) to the World Ocean Atlas 2013 (WOA13; Boyer et al. 2013), without imposing a balance on the water mass transformation budget. They show similar patterns for both TH and TI but estimate about twice the magnitude for TH and 5 times the magnitude of TI as found in this study. We note, however, that WOA13 has a 1° × 1° grid spacing, such that a direct application of Cole et al. (2015) (approximate 3° × 3° grid spacing) may therefore result in an overestimation of the transformation rates by mesoscale mixing.
Improvements in the representation of mixing in the inverse model would improve the estimates of water mass transformation in coordinates. A first step would be to loosen the a priori constraints on the spatial variation of the diffusion coefficients, which results from the structure functions. The spatial variation in mixing can be added by additional unknown diffusion coefficients in the inverse model. This can be achieved by 1) combining the THIM with (traditional) density inverse models (as strongly suggested by Fig. 7), 2) increasing the available information content by adding more tracers, and 3) improving the mixing structure functions. Additional improvements can be made by including the effects of short-wave penetration (Iudicone et al. 2008), geothermal heating (de Lavergne et al. 2016), brine rejection (Nguyen et al. 2009), and the northward extent of a significant freshwater flux from melting sea ice on water mass transformation (Abernathey et al. 2016). We believe that using different ocean climatologies or surface flux products will provide a solution within the expected error of the inverse model.
The sensitivity of the estimates due to the aforementioned mixing choices and structure functions is important to study but is beyond the scope of this paper; here, we focus on the first application of the THIM to observations.
c. Can isopycnal mixing be so small?
Section 5e showed that the estimate of compares well with previous studies using model output, providing a first representation of ocean circulation in (SA, Θ) coordinates from observations. Although there is an imbalance between surface and diffusive fluxes, over the whole (SA, Θ) domain, the inferred diffusivities are well constrained. In addition, and are well within the expected range and leads to water mass transformation rates that are comparable to previous studies.
Previous observationally based studies have not provided a global, full-depth estimate of KI. However, estimates from Cole et al. (2015), between the winter mixed layer and 2000 m, are smaller than their estimate for the deep ocean. Generally, we find that is smaller than expected, as most ocean models tend to use KH = KI, while we found that KH/KI = O(20). This means that the maximum interior mesoscale mixing is about 20 times smaller than that in the mixed layer. Is this a realistic estimate?
To test if KI is small because of the choice for , we have also applied the THIM to estimate one single global small-scale diffusivity D and one global isopycnal diffusivity K, without the use of structure functions (i.e., no spatial variation for K or D was included). This resulted in × m2 s−1 and m2 s−1. The smallness of K motivated a split of K into a mixed layer KH and interior KI component, resulting in a similar estimate for D, but KH = 423 ± 4 m2 s−1 and KI = 12 ± 1 m2 s−1. This leads to a reasonable estimate of D and KH (considering CARS’s spatial resolution), but KI still remains small. This suggest that, with or without using , KI must consistently be much smaller than KH in order to obtain realistic water mass transformation in the deep ocean.
Based on the discussion above, increasing the number of scaling coefficients that are solved for in the vertical will most likely result in higher values of just below the mixed layer, where mode waters are formed (Hanawa and Talley 2001; Iudicone et al. 2011; Groeskamp et al. 2016). For the deep ocean, may reduce even further than currently estimated. This strongly suggests that (with the exception of possible subsurface maxima) decays rapidly with depth, leading to very small values in the deep ocean (below 1500–2000-m depth). This is in strong contrast to numerical ocean models that often retain high values of KI up to 1000 m2 s−1 in the deep ocean.
We have applied the THIM to an ocean climatology in combination with surface flux products, resulting in observationally based estimates of horizontal and isopycnal eddy diffusion coefficients KH and KI, respectively, the small-scale mixing coefficient D, and the diathermohaline streamfunction Ψ, which is representative of the global ocean circulation. The THIM provides an inverse solution that minimizes the divergence of processes based on transformations in individual – classes and gives a balanced global circulation in (SA, Θ) coordinates.
Here, Ψ provides an observational-based quantification of the globally interconnected ocean circulation induced by water mass transformation (Broecker 1991) and 1) shows a robust tropical cell with an extension related to circulation in the Bay of Bengal; 2) suggests the existence of a cell related to the North Pacific (subpolar) Gyre and the Bering Strait outflow, not previously identified by model-based calculations of Ψ; 3) shows a wider spread of the AABW cell compared to model studies; and 4) has a general structure comparable to models. Given that Ψ can now be estimated from observations, we suggest that Ψ has emerged as a unique new diagnostic providing the ability to compare and contrast fundamental properties such as circulation, mixing, and heat and freshwater fluxes between models and observations. The Ψ would be valuable as a standard diagnostic for model intercomparison exercises (Danabasoglu et al. 2014, 2016).
Standard relationships between D, the stratification, and energy dissipation show that 1.34 TW of energy is used for small-scale mixing, of which 0.89 TW is used in the upper 500 m of the ocean and 0.17 TW is used below 2000 m. Hence, most mixing occurs near the surface, but small-scale mixing is important below 2000 m to close the global ocean circulation.
Changes in the magnitude or spatial variation of KI in numerical ocean models have a significant influence on reducing systematic drift (Ferreira et al. 2005; Danabasoglu and Marshall 2007), prediction of the global temperature (Pradal and Gnanadesikan 2014), and the stability of the overturning circulation (Sijp et al. 2006). The results presented here suggest that KH/KI = O(20), meaning that the maximum interior mesoscale mixing is about 20 times smaller than that in the mixed layer. This effect is not included in most ocean models but is potentially important for the aforementioned processes. We conclude that it is very likely that isopycnal eddy mixing in the ocean interior is much smaller than currently assumed, although its vertical structure requires further study. The diffusion coefficient for the eddy-induced, quasi-Stokes (GM) advection and the Redi parameter for tracer diffusion are related, suggesting that the results obtained here may also impact the GM parameterization (Abernathey et al. 2013).
Our results provide the first observationally based, globally constrained, full-depth, isopycnal mixing estimate and reveal unexpectedly small values of the isopycnal diffusivity, compared to that commonly quoted in the literature. Little is known about the vertical structure of isopycnal mixing, while it has a large impact on the sensitivity of the climate and future climate change. Further development of THIM to more accurately estimate this structure from observations is therefore essential.
SG gratefully acknowledges support by the NASA Grant NNX14AI46G. BMS was supported by the Australian Climate Change Science Program, jointly funded by the Department of the Environment and CSIRO. JDZ has been supported by the Natural Environment Research Council Project NE/K012932/1. We thank Daniele Iudicone, Cimarron Wortham, Ryan Abernathey, Andreas Klocker, Stephen Griffies, Jonas Nycander, and Ken Ridgway for useful discussions that have improved the manuscript. We thank Sylvia Cole for providing the mesoscale diffusivity data. We thank Paul Barker for providing preliminary hydrographic stabilization software and Terry O’Kane for providing computational facilities in early stages of this project.
Isopycnal Tracer Gradients
The slopes used in Eq. (10), used to calculate isopycnal SA and Θ gradients, are given by for the x direction and for the y direction. Here, ρl is the locally referenced potential density, that is, it is a potential density referred to the pressure at which one calculated the gradient. This method requires dividing by , which can be zero for a gridded climatology. This leads to spikes in the slope and and a general overestimation of the isopycnal gradients. To avoid singularity, we calculate the slope as and , where and are the differences in the height of the isopycnal between the start and end point in the x and y direction, respectively. The values for , , , and are obtained through an interpolation method, described below.
Tracers are defined at the a T grid at grid point . Here, , where I is the total number of longitude location in the ocean. In a similar way j and J represent latitude, k and K represent depth, and τ and T represent time. The tracer is given in .
Consider only longitude and depth, such that is given at the T grid as . Use averaging to obtain the values at , which is at the edges of the T grid. Start with , where κ is a fixed depth. Using linear interpolation, is obtained for the depth , on profile , for which . A surface has now been constructed, from profile to profile , such that and . In a similar way we can calculate and . The resulting and are given at the middepth , at the midpoint between the two profiles . After repeating this for each i, j, and k, linear interpolation is used to provide the values of the slopes back on the original T-grid depths . If is interpolated through the ocean surface or bottom, the slopes are not calculated.
The Exact Setup of the Inverse Model
Here, we follow up from section 4 to explain in more detail how the inverse method is constructed and applied. An estimate for can be obtained using an inverse technique that minimizes the sum of both the total solution error and the total equation error (Menke 1984; Wunsch 1996; McIntosh and Rintoul 1997):
Above, T is the transpose, is the equation error per equation, is a prior estimate of , is the row (equation) weighting matrix, and is the column (variable) weighting matrix. The matrices and are both diagonal with elements and , respectively; is our best estimate of the error between and ; and is our best estimate of the equation error .
Both elements are chosen to obtain an approximately equal influence of each variable and equation on the solution, as discussed in more detail below. An estimate of the random error on the inverse estimate is given by the square root of the diagonal of the posterior covariance matrix:
Below we describe how to obtain the a priori estimates of the solution and weighting, upon which the results in this manuscript are based. However, the variations in the solution due to the prior choice of the weighting coefficients may be larger than the estimate of the random error associated with that particular choice. A sensitivity study to our a priori choices is therefore provided in appendix C and does not change the general conclusions of this paper.
a. The a priori estimates and constraints
The term is the a priori expected magnitude for , which is unknown as it is very difficult to obtain an accurate a priori estimate of that takes into account both sign and magnitude. We therefore use . We can, however, provide a priori estimates of the diffusion coefficients. In steady state, for a particular water mass, divergence of the transformation by surface forcing [the F terms in Eqs. (7) and (8)] is balanced by convergence due to mixing. It follows that their maximum transformation rates should be of comparable magnitude. As such we use an a priori estimate of the diffusion coefficient such that, when multiplied by the A terms in Eqs. (7) and (8), they do not exceed the water mass transformation given by the boundary forcing (the F terms). For the small-scale diffusion coefficient, we use an a prior estimate of
where is the average over the largest 50 (…) values. We have repeated this for the a priori estimate of the horizontal and isopycnal eddy diffusion coefficients and , respectively (Table 1).
Using the prior estimates of the diffusion coefficients, we have omitted equations that are associated with transports less than 1% of the maximum transport, as these equations do not have a signal to noise ratio that adds information to the solution. We have also omitted equations that are isolated in (SA, Θ) coordinates (i.e., do not bound another bin with significant water mass transformation), as such equations are unconstrained.
b. Row weighting
To construct row weighting , we need to define the a priori estimate of the equation error . We will use the fact that we can calculate the a priori estimates of the volume divergence for each (SA, Θ) grid by combining A, F, and L terms in Eqs. (7) and (8) with the a priori mixing coefficients , , and . The a priori volume divergence shows for which grid cells the a priori solution is further away from nondivergence than others. We will assume that the reason is a higher uncertainty in the water mass transformation estimates, and as such we will we assume that the larger is, the larger the related equation error should be, thereby reducing the influence of equations that contribute to larger a priori volume divergence.
The a priori streamfunction difference is given by
such that for each grid in (SA, Θ) coordinates the a priori volume divergence is given by
Here, reflects the magnitude of the imbalance between water mass formation and destruction (Fig. B1b). The larger this imbalance is, the larger the uncertainty will be in the four associated transport equations (one for each side). Equally distributing the imbalance to each of the four contributing equations, and using that each such transport equation is related to two estimates of (two neighboring grids), we obtain for the SA. An equivalent expression can be obtained for the Θ direction.
We note that, after the inverse solution is obtained, one can use
to calculate (Fig. B1). Here, is calculated as in Eq. (B5) but using instead of . As a result, reflects the magnitude of the imbalance between water mass formation and destruction of the inverse estimate, rather than the a priori estimate (Fig. B1b). The integral of over the whole domain is zero (no net gain of volume). For a perfect solution the volume is exactly conserved in each grid and . In general, , meaning that the inverse model provides a better solution than the a priori estimates (Fig. B1c).
c. Column weighting
We will allow for a standard error of twice the prior estimate of the diffusion coefficients, leading to , , and .
The lack of information for can be compensated by a physically appropriate choice of . We can direct the inversion toward a fit within . As (see beginning of appendix B, section a), we want to be about , such that the structure of is taken into account, and the factor 2 allows for a wide enough range not to restrict the inverse model, while small enough not to make it physically unrealistic.
To obtain an approximate structure for without using prior knowledge of , we will assume that it is likely that if the water mass transformation rates are larger for a particular grid in (SA, Θ) coordinates, than the related values are also likely to be larger. Hence, is proportional to the magnitude of the a priori estimate of the diathermohaline transformations [calculated using Eq. (B4)]. Note that a particular value for is on the corners of each grid in (SA, Θ) coordinates and is thus calculated as the average of four such transformation rates. As is only proportional to the water mass transformations, the resulting structure is normalized with a maximum of 35 Sv to obtain (Fig. B1a).
The results of this particular setup of the inversion are presented in section 5. Results of a sensitivity study to changes in the prior estimates and weighting are provided and discussed in appendix C. For the solution discussed, we have obtained and , both within the range , as suggested in Groeskamp et al. (2014b). As , the solution is obtained with more emphasis on minimizing equation error than minimizing toward the prior estimate.
The formal estimate of the random error of the inverse estimate is given by the square root of the diagonal of the posterior covariance matrix [Eq. (B2)]. However, the variations in the solution due to the prior choice of the weighting coefficients may be larger than the estimate of the random error associated with one particular choice. We will test the sensitivity of our solution to changes in our prior estimates of , , and . We are confident that the scaling argument [Eq. (B3)] provides a good prior estimate of , , and within an order of magnitude. We will use four linearly spaced values for , providing a range of exactly an order of magnitude, centered around the a priori estimate. We will use that , , and , thus scaling with the a priori estimate. This guides the inversion to fit the results within the range . There will be more resistance to fit the solution outside this range.
Although we do not change our a priori estimate of (which is 0), we will also change the a priori estimate of . In the example discussed, we use , where is calculated from the a priori estimate of mixing and is normalized to have a maximum of 35 Sv. Here, we use maxima of 20, 35, and 50 Sv. This turns out to only influence the magnitude of the solution for .
This sensitivity study results in 3 × 43 = 192 inverse solutions that all meet the criteria of , as suggested by Groeskamp et al. (2014b). We have calculated the average and standard deviation (Table 1).
We find that m2 s−1, a difference of almost two orders of magnitude between the extrema. The solution for is most sensitive to and and much less sensitive to the choice of and . The solution for KH ∈ [3800, 10 000] m2 s−1 is most sensitive to (in order of significance) and , while and have little influence on the solution. The solution for m2 s−1 is most sensitive to and somewhat sensitive to .
The extremes for KH is less than a factor 3, while that of D is less than a factor 1.5. This shows KI is more sensitive to changes in the prior estimates. However, we note that even when taking the largest estimate of KI and the smallest estimate of KH, the ratio KH/KI = 3800/1400 = O(3) and still implies significantly reduced magnitude of mesoscale mixing in the interior. Hence, the sensitivity study does not change the key result of this study.