Abstract

A regional thermohaline inverse method (RTHIM) is presented that estimates velocities through the section bounding an enclosed domain and transformation rates resulting from interior mixing within the domain, given inputs of surface boundary fluxes of heat and salt and interior distributions of salinity and temperature. The method works by invoking a volumetric balance in thermohaline coordinates between the transformation resulting from mixing, surface fluxes, and advection, and constraining the mixing to be down tracer gradients. The method is validated using a 20-yr mean of outputs from the NEMO model in an Arctic and subpolar North Atlantic domain, bound to the south by a section with a mean latitude of 66°N. RTHIM solutions agree well with the NEMO model “truth” and are robust to a range of parameters; the meridional overturning circulation (MOC), heat, and freshwater transports calculated from an ensemble of RTHIM solutions are within 12%, 10%, and 19%, respectively, of the NEMO values. There is also bulk agreement between RTHIM solution transformation rates resulting from mixing and those diagnosed from NEMO. Localized differences in diagnosed mixing may be used to guide the development of mixing parameterizations in models such as NEMO, whose downgradient diffusive closures with prescribed diffusivity may be considered oversimplified and too restrictive.

1. Introduction

The complex spatial structure of the global ocean circulation often obscures its key physical phenomena and drivers, including mechanisms involving both internal mixing and external forcing. Lagrangian observations have shown that cartoon-like schematics of the thermohaline circulation may oversimplify the true complexity of the flow pathways (Bower et al. 2009; Lozier 2012), yet it remains desirable to seek a view that is both simple and physically accurate. Thermohaline coordinates have long been used to characterize different water masses. Recent studies have shown that expressing circulation and mixing in this nongeographical coordinate system (e.g., Zika et al. 2012; Hieronymus et al. 2014) filters out adiabatic fluctuations and therefore enables insight into diabatic processes that are relevant to the climate system over long time scales.

Globally, it has been demonstrated that a thermohaline-coordinate circulation streamfunction and internal mixing may be diagnosed from observations or model data (Groeskamp et al. 2014a). A companion study further showed that an inverse method, an extension of Walin (1982), may be used to infer the circulation and internal mixing from knowledge of the temperature and salinity combined with surface fluxes—the thermohaline inverse method (THIM) (Groeskamp et al. 2014b, hereafter G14b). THIM elucidates the connectivity of the global transport of heat and freshwater and their drivers, but the trade-off is that detailed regional information may be lost through the transformation from geographical to thermodynamic coordinates, due to the integration of water masses with the same temperature (θ)–salinity (S) properties.

Motivated by the need to understand regional circulation and mixing, this study describes a new regional thermohaline inverse method (RTHIM) and its validation against output from the Nucleus for European Modelling of the Ocean (NEMO) model (Madec 2008), with a particular application to the Arctic and subpolar North Atlantic (Fig. 1). The equivalent transformation in the regional domain once again achieves a simplified view of the complex circulation and mixing, but the integration of water masses retains information specific to the region of interest, rather than this information being a smaller fraction of a global integral. The outputs from traditional “box inverse” models (Wunsch 1978) have been analyzed in thermohaline coordinates to provide insight into the circulation in specific regions such as the Weddell Gyre (Naveira Garabato et al. 2016) and the Arctic Ocean (Tsubouchi et al. 2012), but RTHIM is unique as an inverse method that directly uses a thermohaline-coordinate system applied to a specific region of the ocean, rather than the global ocean, as in THIM. In THIM, water of particular Sθ properties may exist in more than one region. Because of the averaging process, such regional distinction is lost, but in RTHIM the diagnostics remain regionally specific. There is a trade-off between sufficient averaging to obtain a simpler view and dilution of regional information by excessive averaging. Both approaches are arguably valid, depending on the particular question or goal. In addition, previous methods that have solved for interior mixing have done so by imposing a simplified spatial structure on the diagnosed diffusivities. RTHIM does not impose such a structure, instead applying only simple constraints on mixing, which we describe in section 2b.

Fig. 1.

Plan view of the Arctic and subpolar North Atlantic domain to which RTHIM is applied. The background colors show the 1988–2007 time mean of the surface fluxes of (left) freshwater and (right) heat forcing NEMO. Fluxes are positive going into the ocean. The section bounding the domain, defined at a constant pseudo-latitude (red line). The model’s tripolar grid means that the section latitude varies with longitude, with the mean for the section equal to 65.6°N. The white contour is the zero line for the fluxes.

Fig. 1.

Plan view of the Arctic and subpolar North Atlantic domain to which RTHIM is applied. The background colors show the 1988–2007 time mean of the surface fluxes of (left) freshwater and (right) heat forcing NEMO. Fluxes are positive going into the ocean. The section bounding the domain, defined at a constant pseudo-latitude (red line). The model’s tripolar grid means that the section latitude varies with longitude, with the mean for the section equal to 65.6°N. The white contour is the zero line for the fluxes.

Climate models suggest that buoyancy changes in the Arctic and subpolar North Atlantic strongly influence the Atlantic meridional overturning circulation, yet observational evidence is actively being sought to support the model predictions (Lozier et al. 2017). This paper is a first step toward providing relevant robust observational estimates based on RTHIM, which will be used to complement more traditional direct observations provided by new transbasin arrays, such as in the Overturning in the Subpolar North Atlantic Program (OSNAP; www.o-snap.org).

The paper is organized as follows: In section 2 we describe the principles of the new inverse method and detail the background to its development. Section 3 explains how we have applied RTHIM to outputs from the NEMO model for validation. In section 4 we present the results of the validation, and section 5 contains a discussion and conclusions.

2. The regional THIM

a. Distributions in thermohaline coordinates

When working in thermohaline coordinates, it is helpful to consider the distribution of water masses in Sθ space. For example, if we are looking at the transport through a section (e.g., the section indicated by the red line in Fig. 1), the section itself can be divided up using contours of S and θ (Fig. 2a). The section area distribution in Sθ space from a 20-yr time mean of the NEMO fields (Fig. 2b) is calculated by summing the areas of the section that fall between the contours of S and θ, defined as S ∈ (S ± ΔS/2), θ ∈ (θ ± Δθ/2) where ΔS = 0.1 psu and Δθ = 0.2°C. Similarly, the volumetric distribution (Fig. 2c) is calculated by summing the volumes enclosed by contours of S and θ over the whole physical domain. Both distributions are concentrated in comparatively small regions of Sθ space. For the area distribution, 76% of the total section area of 2.6 × 109 m2 is covered by 10% of the occupied (S, θ) bins (i.e., 10% of the bins that have a nonzero associated section area). For the volumetric distribution, 95% of the total volume of 1.8 × 1016 m3 is contained within 10% of the occupied bins.

Fig. 2.

(a) The 1988–2007 time-mean temperature (color shades) and salinity (white contours) from the section defined in Fig. 1. The two basins of the main plot are the (left) Labrador Sea and (right) GIN Seas; the inset shows Bering Strait. (b) A zoomed-in part of the section in the western Norwegian basin with a red polygon indicating a region enclosed by pairs of S and θ contours. (c) The distribution in Sθ space of the section area. The grid spacing is ΔS = 0.2 psu and Δθ = 0.5°C, the same as the contour intervals in (a) and (b). The inset shows a zoomed-in part of the distribution with a red rectangle corresponding to the red polygon in (b). (d) The volumetric distribution in Sθ space for the domain bounded by the section [note that the color scales in (c) and (d) are logarithmic]. A more detailed description of the domain is given in section 3a.

Fig. 2.

(a) The 1988–2007 time-mean temperature (color shades) and salinity (white contours) from the section defined in Fig. 1. The two basins of the main plot are the (left) Labrador Sea and (right) GIN Seas; the inset shows Bering Strait. (b) A zoomed-in part of the section in the western Norwegian basin with a red polygon indicating a region enclosed by pairs of S and θ contours. (c) The distribution in Sθ space of the section area. The grid spacing is ΔS = 0.2 psu and Δθ = 0.5°C, the same as the contour intervals in (a) and (b). The inset shows a zoomed-in part of the distribution with a red rectangle corresponding to the red polygon in (b). (d) The volumetric distribution in Sθ space for the domain bounded by the section [note that the color scales in (c) and (d) are logarithmic]. A more detailed description of the domain is given in section 3a.

Processes acting on and within the ocean affect the distribution of water masses in Sθ space (Groeskamp et al. 2014a; Pemberton et al. 2015). Surface heat fluxes increase the spread of water masses along the θ axis, while freshwater fluxes increase their spread along the S axis. Interior mixing processes on the other hand act to reduce the spread of water masses in Sθ space, since mixing acts down tracer gradients. In addition to the most familiar effects of heating–cooling and freshening–salinification, ice melting causes water masses below to cool. This is particularly relevant along the Atlantic water inflow (Fram Strait and around Svalbard), where surface waters are relatively warm and sea ice melting on top creates large cold freshwater lenses, which get mixed during storms. The competition between boundary and interior processes allows us to diagnose interior mixing in the absence of direct observations, by inferring it from other quantities that are more easily measured (see, e.g., Zika et al. 2015).

b. The volume budget

The aim is to estimate, as outputs from the method, section volume advection and interior diffusive processes (described later), given inputs of surface boundary fluxes and interior S and θ. As in G14b, we begin by considering a volume element bound by pairs of isohaline and isothermal surfaces at time t, V = V(S ± ΔS/2, θ ± Δθ/2, t) (Fig. 3). G14b have Absolute Salinity SA and Conservative Temperature Θ as their coordinates, which is appropriate for application to real ocean observations; here we use practical salinity S and potential temperature θ, which are conserved in NEMO. The conservation equations for volume, salt, and heat mirror those in Eqs. 7–9 of G14b for the global THIM, except that we express the conservation per unit S and θ, and we add terms involving Iadv = Iadv(S, θ), which we define as the volume advection (inward normal velocity component times the area of volume element at the bounding section) into the regional volume at its open boundary per unit S and θ, leading to the following:

 
t(V)+SθUSθdia=1ρ0fm+Iadv,
formula
(1)

where ∂(V)/∂t is the volume tendency per unit (S, θ), ∇USθdia is the diathermohaline volume transport divergence per unit (S, θ), ρ0 = 1024 kg m−3 is a reference density of seawater, fm is the surface boundary mass flux per unit (S, θ), and Iadv is the section volume advection per unit (S, θ);

 
t(VS)+Sθ(SUSθdia)=1ρ0mS+IadvS,
formula
(2)

where ∂(VS)/∂t is the salt tendency per unit (S, θ), ∇(SUSθdia) is the diathermohaline salt transport divergence per unit (S, θ), mS is the salt convergence as a result of interior diffusion per unit (S, θ), and IadvS is the section salt advection per unit (S, θ); and

 
t(Vθ)+Sθ(θUSθdia)=1ρ0cp0(fh+mθ)+Iadvθ,
formula
(3)

where ∂()/∂t is the heat tendency per unit (S, θ), ∇(θUSθdia) is the diathermohaline heat transport divergence per unit (S, θ), cp0 is the specific heat capacity of seawater (4000 J kg−1 K−1), fh is the heat convergence as a result of boundary fluxes per unit (S, θ), mθ is the heat convergence as a result of interior diffusion per unit (S, θ), and Iadvθ is the section heat advection per unit (S, θ). Note that the superscript dia denotes a transport through a surface of constant temperature or salinity.

Fig. 3.

A volume element V illustrating the RTHIM volume budget. The element can be defined between pairs of isotherms (θ ± Δθ/2) and isohalines (S ± ΔS/2), such as is indicated by the polygon in the inset of Fig. 2b. The quantity Iadv(S, θ) is the advective flux through the section into the volume element per unit S, θ; f¯m and f¯h are boundary fluxes of mass and heat, respectively; and m¯S and m¯θ are diffusive fluxes of salt and heat, respectively. The change in volume of the element depends on the movement of its bounding surfaces, Udia, and is governed by Eq. (9). Figure adapted from Groeskamp et al. (2014b).

Fig. 3.

A volume element V illustrating the RTHIM volume budget. The element can be defined between pairs of isotherms (θ ± Δθ/2) and isohalines (S ± ΔS/2), such as is indicated by the polygon in the inset of Fig. 2b. The quantity Iadv(S, θ) is the advective flux through the section into the volume element per unit S, θ; f¯m and f¯h are boundary fluxes of mass and heat, respectively; and m¯S and m¯θ are diffusive fluxes of salt and heat, respectively. The change in volume of the element depends on the movement of its bounding surfaces, Udia, and is governed by Eq. (9). Figure adapted from Groeskamp et al. (2014b).

We begin the separation into inputs and outputs by identifying the component terms of the time-mean diathermohaline volume transport USθdia=(U|Sdia,U|θdia). Here U|Sdia is the transport across an S surface between two θ surfaces per unit θ; and U|θdia is the transport across a θ surface between two S surfaces per unit S. We repeat the algebraic manipulation used in G14b on the analogous equations [Eqs. (1)(3)] but using a different definition of the thermohaline divergence operator (and similar for the gradient operator), ∇ ⋅ [⋅] = (∂/∂S, ∂/∂θ) ⋅ [⋅].The discrete ∂/∂S operator is defined at the central S value of the grid cell corresponding to (S ± ΔS/2, θ ± Δθ/2) and analogous for θ. Note that ∂S/∂S = 1 and ∂θ/∂S = 0; this is similar to when S is swapped with θ. The rightmost terms involving Iadv in Eqs. (1)(3) are found to cancel, leading to the following expressions for each of the components in thermohaline coordinates (an overbar represents a time mean):

 
U|Sdia(S,θ)=U|Ssurf+U|SmixU|Sloc,
formula
(4)
 
U|θdia(S,θ)=U|θsurf+U|θmixU|θloc,
formula
(5)

where

 
U|Ssurf=1ρ0Sfm¯,U|θsurf=1ρ0cp0fhcp0θfm¯
formula
(6)

are due to surface processes;

 
U|Smix=1ρ0m¯S,U|θmix=1ρ0cp0m¯θ
formula
(7)

are due to interior diffusive processes, defined in more detail later; and

 
U|Sloc=1ρ0lS¯,U|θloc=1ρ0cp0lθ¯
formula
(8)

are related to the local response, lS¯=ρ0VS/t¯ and lθ¯=ρ0cp0Vθ/t¯. The local response expresses the rate of change of heat (salt) as a result of the rate of change of temperature (salinity), for a given volume within a discrete (S, θ) grid cell. Equivalently, it is the amount of heat (salt) that can be added or removed from the system, changing the temperature (salinity) of a water mass by an amount whereby it remains within that grid cell.

We may now use Eqs. (4) and (5) to rewrite the time mean of Eq. (1) in terms of desired output estimates (lhs) balanced by specified inputs (rhs):

 
IadvSθUSθmix=Vt1ρ0fm¯+SθUSθsurfSθUSθloc.
formula
(9)

For the sake of simplicity, we follow the approach of Groeskamp et al. (2017) and neglect the local response term, since it is small for the (S, θ) grids used in this study. In addition, the mass flux term fm, is also negligible for the NEMO solution used here so that is omitted too.

At this point we again deviate from G14b by including a more general formulation for the interior mixing, which no longer depends on the assumption of a constant diffusivity with a fixed geographical distribution. We may write the diathermohaline volume transport vector as a result of interior diffusive processes USθmix, as the divergence of a tensor F, whose properties we may constrain to ensure net downgradient mixing:

 
USθmix=SθF=(SFSS+θFSθ,SFθS+θFθθ),
formula
(10)

where F is constrained to be a symmetric diffusive flux tensor in thermohaline coordinates (see the  appendix for an explanation of the symmetric constraint).

Combining this definition of the diffusive flux tensor with Eq. (9), we may define the inverse problem,

 
IadvSθ2F=Vt+SθUSθsurf+ε,
formula
(11)

where ε is the error in the inverse solution, which we aim to minimize. Equation (11) is our volume budget, representing the fact that the time-mean change in volume of a volume element results from a balance between advection, surface fluxes, and interior mixing.

We solve Eq. (11) for Iadv and the elements FC1C2 of F using constrained optimization (MATLAB’s fmincon), with constraints FSS, Fθθ ≥ 0 (net downgradient mixing), and F = FθS (symmetric diffusive flux tensor). Note that because Eq. (11) involves the Sθ2 operator applied to F, it is possible to conceive of an arbitrary, nondivergent gauge N, to be added to F—for example, F′ = F + N where Sθ2N=0 and F′ is the modified F —yet to still obtain the same volume balance: Sθ2F=Sθ2F. Often a Helmholtz decomposition is used to solve similar gauge problems, given appropriate boundary conditions on F. However, there is no appropriate physical metric in thermohaline coordinates for such a decomposition [in Cartesian coordinates (x, y), one metric could be a measure of distance, proportional to (Δx2 + Δy2). Therefore, we emphasize that the interpretation of the RTHIM solution should be focused solely on the role of SθUSθmix, rather than on the individual terms in F.

3. Application to NEMO

a. NEMO model domain

The model data used for the validation of RTHIM come from a 1°-resolution NEMO simulation run at the National Oceanography Centre (NOC) known as ORCA1-N403. It is a z-level Boussinesq global model, with the OPA ocean model (Madec 2008) coupled to the Louvain-la-Neuve Sea Ice Model, version 2 (LIM2) (Timmermann et al. 2005). The horizontal grid is tripolar, giving a resolution of ~50 km over the Arctic Ocean. There are 75 vertical levels ranging in thickness from ~1 m at the surface to ~200 m at the bottom. The model is forced by CORE atmospheric forcing (Large and Yeager 2004), which has 2.5° resolution. We took monthly mean outputs of net surface fluxes of heat and freshwater and of ocean temperature and salinity from 20 yr of the simulation from 1988 to 2007 and averaged them to produce a monthly climatology. The limited resolution of the simulation means that mesoscale features of the circulation are not resolved; however, our focus is to validate the inverse method against an internally consistent system where subgrid mixing is parameterized.

The Arctic and subpolar North Atlantic domain volume to which we apply RTHIM is bounded to the south by a line of constant latitude index on the ORCA1-N403 grid. The model has a tripolar grid, so the constant latitude index line is curved in latitude–longitude coordinates (see red lines in Fig. 1). The flow into and out of the domain is through a full-depth section defined by this line. In order for the balance of Eq. (11) to hold, the water masses that flow through the section must be geographically connected within the model so that they can mix. We therefore exclude ranges from the section corresponding to the Hudson Bay and part of the Canadian Archipelago (68°–94°W; although the latter is connected to the rest of the Arctic in the real ocean, it is not connected in the model), and the Gulf of Bothnia (18°–22°E). We also mask out the surface fluxes and interior water masses associated with these ranges, since these do not contribute to the volume budget for the domain.

b. Sθ grid

We design an irregular S–θ grid for use with RTHIM (see the  appendix for grid node definitions). We saw in Fig. 2 that the majority of the area of the section occupies a small region in S–θ space. The along-section distance Δx that can be resolved in the RTHIM solution is related to the θ grid spacing Δθ by Δx ≈ (∂θ/∂x)−1Δθ, where ∂θ/∂x is the horizontal gradient of θ (the equivalent applies in the S coordinate). An irregular grid allows us to resolve details of the flow through the section while keeping the size of the inverse problem manageable. We want to focus the resolution in the regions of S–θ space where gradients of S and θ are small, with additional focus given to parts of the section where the transports are large. To do this, we first construct cumulative transport functions in S and θ that measure the magnitude of the transport through the section below a contour, based on the absolute model velocities |υ| normal to the section at each model grid cell:

 
|T|S(S)=SminS|υ|ΔA,|T|θ(θ)=θminθ|υ|ΔA.
formula
(12)

Here the sums are over the area of the section for which the salinity (temperature) is below the salinity (temperature) in the bracket. We then plot the transport functions normalized by their totals, which allows the θ and S ranges to be discretized according to an equipartition of |T| (Fig. 4). Each point in θ (or S) identified by the red dashed lines in Fig. 4b (Fig. 4c) as they intersect the x axis then becomes a node in the θ (S) vectors used to construct the S–θ grid.

Fig. 4.

(a) The section θ (colors and black contours) and S (white contours) from the NEMO simulation time mean. As in Fig. 2a, except that the contour intervals have been determined for use in the RTHIM Sθ grid using the transport functions in (b) and (c) as detailed in section 3b. The lines in (b) and (c) represent the normalized transport functions (blue lines), and indicate the selection of the θ and S grid nodes (red dashed lines). In this case, the transport functions have been divided up into 12 approximately equal segments; however, we can choose to construct a grid using more or fewer nodes. The grids constructed from these plots, including additional nodes at either end to ensure binning of the whole range of (S, θ) found within the control volume, are detailed in the  appendix.

Fig. 4.

(a) The section θ (colors and black contours) and S (white contours) from the NEMO simulation time mean. As in Fig. 2a, except that the contour intervals have been determined for use in the RTHIM Sθ grid using the transport functions in (b) and (c) as detailed in section 3b. The lines in (b) and (c) represent the normalized transport functions (blue lines), and indicate the selection of the θ and S grid nodes (red dashed lines). In this case, the transport functions have been divided up into 12 approximately equal segments; however, we can choose to construct a grid using more or fewer nodes. The grids constructed from these plots, including additional nodes at either end to ensure binning of the whole range of (S, θ) found within the control volume, are detailed in the  appendix.

c. Advective flux

During the constrained optimization, we want RTHIM to compute the advection term Iadv from a surface reference-level velocity υref. The Iadv is calculated by assuming a relative velocity υref to be added to υref:

 
υ(x,z)=υref(x)+υrel(x,z)
formula
(13)
 
Iadv(S,θ)=Iref(S,θ)+Irel(S,θ),
formula
(14)

where (x, z) are coordinates in the along-section and vertical directions, respectively. The components of Iadv are found by binning relative velocities at the section:

 
Irel(S,θ)=1ΔSΔθΠ(S)Π(θ)υreldA
formula
(15)

(and similarly for Iref), where υrel(x,z)=υ(x,z)υsurf(x) is the NEMO section velocity relative to the fixed surface velocity. The Π mask only in regions of the section between isohalines S ± ΔS/2 and isotherms θ ± Δθ/2, and the quantities are integrated over the whole section. We construct a (nS × nθ)-by-(nx) matrix A such that Iref = Aυref, where Iref and υref are (nS × nθ)- and (nx)-element column vectors, respectively, and where nS, , and nx are the number of bins used to discretise the S, θ, and x grids. Then Iadv(S, θ) can be calculated using Eq. (14).

In the initial condition, υref is taken as the model values with random noise added between −5 and +5 cm s−1 with five-cell boxcar smoothing applied (equivalent to smoothing over 5° longitude). This initial condition is designed to reflect the uncertainty on the surface velocities that might be expected when applying RTHIM to observations and using surface geostrophic velocity from satellite altimetry data for the initial υref. The 5 cm s−1 is slightly larger than the 3 cm s−1 estimated by Gourcuff et al. (2011) as the uncertainty on altimetrically derived surface geostrophic velocities, but it also represents an upper limit on the uncertainty that might be expected from in situ current meter measurements (see, e.g., Tsubouchi et al. 2012). We also constrain υref in the solution to be within 10 cm s−1 of the NEMO model “truth,” which limits the range of solutions while allowing for a generous amount of uncertainty on our “best guess” by doubling the expected value.

d. Surface flux

Here we describe how the fluxes of heat and freshwater that forced the NEMO run are used to calculate SθUSθsurf for input to RTHIM. The contribution of the surface fluxes to the volume budget, SθUSθsurf, is a prescribed term in the RTHIM solution, assumed perfectly known. This is appropriate for validating with a model that is an internally consistent system; we discuss the implications for applying RTHIM to observations in section 5. The freshwater flux (m s−1) (see Fig. 1, left panel) is divided by the thickness of the model’s top layer, Δzsurf, obtaining fw (s−1), and U|Ssurf is then calculated as follows:

 
U|Ssurf(S*,θ*)=S0ΔSΔθΠ(S*)Π(θ*)fwdV.
formula
(16)

Here, the Π masks in regions on the surface layer between isohalines S* and S* + ΔS and between isotherms θ* ± Δθ/2. The different masking here compared with Eq. (15) for Iadv ensures that ∇USθsurf will coincide with Iadv on the discrete grid, since U|Ssurf is calculated between S nodes and on θ nodes, while U|θsurf (below) is calculated on S nodes and between θ nodes; therefore, SθUSθsurf = (∂/∂S, ∂/∂θ) ⋅ (U|Ssurf,U|θsurf) will be on S and θ nodes, as is Iadv. The heat flux (W m−2) (Fig. 1, right panel) is divided by (ρ0cp0Δzsurf=4.1×106Δzsurf) to obtain fθ (K s−1), and U|θsurf is calculated as follows:

 
U|θsurf(S*,θ*)=1ΔSΔθΠ(S*)Π(θ*)fθdV
formula
(17)

The Π masks regions on the surface layer between isohalines S* ± ΔS/2 and between isotherms θ* and θ* + Δθ.

e. Diffusive flux

We described constructing an initial condition for the reference-level velocities in section 3c. We also require an initial condition for the components of F. Using the NEMO fields, we calculate the x, y, and z components of the along-isopycnal gradients of S and θ and their vertical gradients, assuming small isopycnal slopes (note that this assumption applies only to the initial condition for F and is not required for the solution). We then calculate the four components of F by defining a diffusion tensor K with a uniform isopycnal and vertical component for the purposes of the initial condition (see the  appendix). We also apply some smoothing to the components of F in the initial condition, since with our variable Sθ grid they can be noisy where ΔS or Δθ is small. We therefore apply a simple 2D boxcar smoothing to these terms. Along with the surface flux term SθUSθsurf, the diffusive fluxes are calculated for each month of our NEMO monthly climatology, and the mean of each set of 12 fields is calculated for input to RTHIM.

f. Optimization

Having established the inputs to Eq. (11), RTHIM can now be solved by minimizing the cost function r2 (in Sv2; 1 SV ≡ 106 m3 s–1) in the following:

 
r2=S,θ(εΔSΔθ)2+w(S,θIadvΔSΔθTnet)2,
formula
(18)

where the term in parentheses on the rhs is a constraint on the net transport through the section, calculated as the difference between the integrated section transport from the RTHIM solution advection term Iadv and the known 1988–2007 time-mean NEMO net section transport including the Bering Strait (Tnet = 0.41 Sv). This net section transport constraint is multiplied by a weighting factor w that adjusts the relative importance of the constraint compared to that of the volume budget (the sensitivity of our results to the factor w is explored in section 4b). When applying RTHIM to observations, Tnet would be set to zero or a known value; since we have an enclosed domain, a zero net transport would be an appropriate assumption in the absence of other information.

An initial control vector is constructed in two sections: the first from the four components of F, which are arrays of values at each S and θ on our grid; and the second from υref. Each section is normalized by the mean of that section’s values to make all the values in the control vector of order of 1; this enables the minimization algorithm to more efficiently search the parameter space. Within MATLAB/fmincon we use the interior point algorithm to carry out the minimization (Byrd et al. 1998). There is more than one option for terminating the optimization on arrival at a satisfactory solution. We explored using the default criterion based on the cost function gradient becoming smaller than a certain value (it should tend to zero at an extremum), but this was susceptible to gradient noise and oscillations in the solution. Instead we used a criterion based on the cost function value itself. When r2 drops below 10−4 Sv2, we accepted the solution (corresponding to a misfit r of 0.01 Sv, chosen to be significantly less than observational uncertainty of hydrographic sections or transbasin mooring arrays).

4. Solution

a. Section transport

In this section we will examine the part of the RTHIM solution that diagnoses the transport through the section bounding the domain. We quantify the method’s skill in reproducing the model truth from NEMO using a number of metrics and test the sensitivity of the solution to a range of reasonable parameters. The purpose is to demonstrate that RTHIM will find a good solution given, for example, an uncertainty on υref that reflects that which would be present when applying the method to observations (i.e., the uncertainty in the surface absolute geostrophic velocities).

RTHIM finds a surface υref by solving Eq. (11), given an initial condition obtained by taking the surface velocities from NEMO and adding some random noise. Given such an initial guess, RTHIM has sufficient skill to converge on the model truth with a high degree of precision (top panel of Fig. 5). The rms error between the RTHIM υref and the model truth for this solution was |υrms| = 0.06 cm s−1. The longitude range shown corresponds to the part of the section indicated by the red line in Fig. 1 between Baffin Island to the west and Norway to the east. The longitude range west of Baffin Island, including Hudson Bay, is not part of our solution as explained in section 3a; the Bering Strait is included in the solution but is not plotted. We then obtain the full section velocities by adding the known model shear to υref (in applying RTHIM to observations the shear will be obtained from thermal wind balance). Even with a high degree of agreement between the RTHIM and NEMO surface velocities, it is still possible to see differences in the full section velocities in the deep part of the section east of the Icelandic Plateau (note that because we use the actual NEMO model shear to construct the section velocities from υref, any differences at depth must be due to differences at the surface). This highlights the importance of obtaining a good solution for υref, particularly where the ocean is deepest. We define an additional skill measure comparing the section transport from the RTHIM solution as follows:

 
Tskill=1Nx=1N[zTsol(x,z)zTmod(x,z)]2,
formula
(19)

where Tsol(x, z) = υsol(x, z)dA(x, z) and Tmod(x, z) = υmod(x,z)dA(x,z) are the volume transports (Sv) through the geographical grid cell (x, z) at the section from RTHIM and NEMO, respectively. The sums over z in Eq. (19) are over the water column, and N is the number of x grid cells along the section. The term Tskill therefore measures the rms error (RTHIM vs NEMO) for the depth-integrated transports at each grid cell in x. For the RTHIM solution shown in Fig. 5, Tskill = 0.03 Sv, which is less than 4% of the rms of the NEMO depth-integrated transports along the section.

Fig. 5.

Velocities corresponding to the section identified by the red lines in Fig. 1. (top) Surface υref from RTHIM solution (red) against the NEMO model truth (black) with the RTHIM initial condition (blue). (middle) Full section velocities from the RTHIM solution obtained by adding the known model shear to the solution υref. (bottom) Full section velocities from NEMO.

Fig. 5.

Velocities corresponding to the section identified by the red lines in Fig. 1. (top) Surface υref from RTHIM solution (red) against the NEMO model truth (black) with the RTHIM initial condition (blue). (middle) Full section velocities from the RTHIM solution obtained by adding the known model shear to the solution υref. (bottom) Full section velocities from NEMO.

The RTHIM solution υref is robust to a range of parameters. If we take an ensemble of initial conditions (blue lines in Fig. 6, top panel) while varying a number of other model parameters, which will be described in the next subsection, the solutions (red lines in Fig. 6, top panel) cluster around the model truth (black line). We can further demonstrate the model skill in solving for υref by examining the solution as it converges (Fig. 6, bottom panel). With a convergence tolerance of r2 = 1, the solution (solid blue line) has deviated from the initial condition (dotted blue line) in only a few places; however, we note that in this initial period of optimization, r2 has reduced from ~106 to <1. This is achieved mainly through adjustments of the elements of F, with some small adjustments to υref required to reduce the net transport constraint in Eq. (18) from ~103 to ~10−5. As the convergence tolerance is decreased, the υref solutions converge toward the model truth, with a tolerance of 10−4 yielding a solution (solid green line) that almost matches the model (dotted black line).

Fig. 6.

(top) The quantity υref from an ensemble of RTHIM solutions (light red lines) with varying initial conditions (light blue lines) and varying model parameters (described in section 4b), compared with the model truth (black line). (bottom) The convergence of υref from its initial condition (blue broken line) toward the model truth (black broken line). Each different colored solid line is an RTHIM solution with a different convergence tolerance applied to the optimization algorithm.

Fig. 6.

(top) The quantity υref from an ensemble of RTHIM solutions (light red lines) with varying initial conditions (light blue lines) and varying model parameters (described in section 4b), compared with the model truth (black line). (bottom) The convergence of υref from its initial condition (blue broken line) toward the model truth (black broken line). Each different colored solid line is an RTHIM solution with a different convergence tolerance applied to the optimization algorithm.

b. Solution sensitivities

We wish to establish some useful metrics for the RTHIM solution appropriate to a wider oceanographic context and assess the sensitivity of the solution to a range of parameters in terms of those metrics. First we define familiar streamfunctions of the section transport diagnosed from the RTHIM solution, binned in terms of density, temperature, and salinity:

 
ψσ(σ*)=σσ*υdA,
formula
(20)

where υdA is the transport through the section, in this case integrated below contours of constant potential density σ2 referenced to 2000 m . The similar streamfunctions Ψθ(θ*) and ΨS(S*) are defined in the same way but integrated below contours of temperature and salinity. We then calculate the meridional overturning circulation (MOC), the meridional heat transport (MHT), and the meridional freshwater transport (MFWT) through the section as follows:

 
MOC=σminσmaxT(σ)forT(σ)>0
formula
(21)
 
MHT=ρ0cp0θminθmaxψθdθ
formula
(22)
 
MFWT=1S0SminSmaxψSdS,
formula
(23)

where in this case T(σ) represents the transport integrated between density contours 0.001 kg m−3 apart, ρ0cp0=4.1×106Jm3K1, and S0 = 35 psu is a reference salinity. Note that we calculate the MOC from T′(σ) as in Li et al. [2017, their Eq. (7)] rather than from the streamfunction.

There are a number of parameters that can affect the RTHIM solution. We have already mentioned adding random noise to the υref initial condition. We also produce an initial condition for the components of F, the other unknown in our inverse problem, based on a uniform diffusivity K made up of along isopycnal and vertical components, K and D, respectively. We test K values of 10, 20, 50, 100, 200, 500, 1000, and 2000 m2 s−1; and D values of 10−5 and 10−4 m2 s−1. There is also the 2D boxcar smoothing applied to the components of F in the initial condition, which has an associated parameter for the number of cells over which to apply the smoothing (we smooth uniformly in both the S and θ directions). We test integer values of 3, 5, 7, and 9 of this parameter, and an additional case where we simply replace each element of FC1C2(S,θ) with the mean of that component over all Sθ space. In section 3b we described the design of an irregular Sθ grid over which we bin all the terms of our volume budget, to focus the available resolution where it is needed by dividing up the transport through the section equally in S or θ coordinates. We can choose the number of segments to divide the ranges into, and test grids of 10 × 10, 13 × 13, 15 × 15, and 17 × 17; at higher resolution than this the number of possible solutions becomes very large and RTHIM has difficulty converging. Finally, the constraint on the net section transport in Eq. (18) has a w that we vary, testing values of 0, 1, 10, 20, 50, 100, 200, 500, and 1000. We summarize the sensitivity to all these parameters of the streamfunctions and metrics defined above in Fig. 7 and Table 1.

Fig. 7.

Streamfunctions (Ψ) in (left) σ, (middle) θ, and (right) S coordinates calculated from RTHIM solution section velocities from ensembles varying across a range of parameters. The same-colored pairs of dashed lines give the upper and lower limits at each point on the y axis of the streamfunctions for an ensemble. The parameters varied were the υref initial condition (red lines), the isopycnal and diapycnal mixing coefficients applied to the F initial condition (purple lines), the Sθ grid used for binning the terms in the volume budget (dark blue lines), the smoothing factor applied to the F initial condition (green lines), and the weighting factor on the net transport constraint (light blue lines). The streamfunctions calculated from the NEMO section velocities are shown (black). Note that the x-axis ranges have been set to assist the reader in discerning details of the contributions of different parameters to the uncertainty, for which it is necessary not to include the full range of the K and D lines. Their full extent can be seen (solid red lines in Fig. 8). For Ψσ, the lower bounds have contributions from RTHIM runs 6, 12, and 13, and are detailed in Table A1; the upper bounds have contributions from runs 2 and 13. For Ψθ, the lower bounds have contributions from runs 2 and 13; the upper bounds have contributions from runs 6, 11, and 13. For ΨS, the lower bound comes from run 13 and the upper bound from run 5.

Fig. 7.

Streamfunctions (Ψ) in (left) σ, (middle) θ, and (right) S coordinates calculated from RTHIM solution section velocities from ensembles varying across a range of parameters. The same-colored pairs of dashed lines give the upper and lower limits at each point on the y axis of the streamfunctions for an ensemble. The parameters varied were the υref initial condition (red lines), the isopycnal and diapycnal mixing coefficients applied to the F initial condition (purple lines), the Sθ grid used for binning the terms in the volume budget (dark blue lines), the smoothing factor applied to the F initial condition (green lines), and the weighting factor on the net transport constraint (light blue lines). The streamfunctions calculated from the NEMO section velocities are shown (black). Note that the x-axis ranges have been set to assist the reader in discerning details of the contributions of different parameters to the uncertainty, for which it is necessary not to include the full range of the K and D lines. Their full extent can be seen (solid red lines in Fig. 8). For Ψσ, the lower bounds have contributions from RTHIM runs 6, 12, and 13, and are detailed in Table A1; the upper bounds have contributions from runs 2 and 13. For Ψθ, the lower bounds have contributions from runs 2 and 13; the upper bounds have contributions from runs 6, 11, and 13. For ΨS, the lower bound comes from run 13 and the upper bound from run 5.

Table 1.

A summary of the sensitivities of RTHIM metrics to various parameters. Each row details an ensemble of runs where a particular parameter (column 1) is varied within a range (column 2). All other parameters are kept constant while exploring a particular parameter. The remaining columns give the metrics of the rms difference between NEMO and RTHIM reference-level velocities and depth-integrated transports υrms and Tskill; and the MOC, MHT, and MFWT as described in section 4. Where range appears in the column head, the lower and upper values of the metric from the ensemble are given. Where a value appears in italics in the column head, these are the NEMO values. The parameters are the number of nodes in S and θ on the RTHIM Sθ grid, the vertical and lateral diffusivities D and K are applied to the diffusive flux initial condition, the w is applied to the section net transport constraint, the random noise is applied to the υref initial condition, and the number of cells of boxcar smoothing is applied to the diffusive flux initial condition. The MOC and MFWT are in Sverdrups. The MHT is in petawatts. The last row of the table summarizes the results, excluding rows 1 and 4. The breakdown of the results used to construct this table is given in Table A1. The “I. C.” with the υref means initial conditions.

A summary of the sensitivities of RTHIM metrics to various parameters. Each row details an ensemble of runs where a particular parameter (column 1) is varied within a range (column 2). All other parameters are kept constant while exploring a particular parameter. The remaining columns give the metrics of the rms difference between NEMO and RTHIM reference-level velocities and depth-integrated transports υrms and Tskill; and the MOC, MHT, and MFWT as described in section 4. Where range appears in the column head, the lower and upper values of the metric from the ensemble are given. Where a value appears in italics in the column head, these are the NEMO values. The parameters are the number of nodes in S and θ on the RTHIM S–θ grid, the vertical and lateral diffusivities D and K are applied to the diffusive flux initial condition, the w is applied to the section net transport constraint, the random noise is applied to the υref initial condition, and the number of cells of boxcar smoothing is applied to the diffusive flux initial condition. The MOC and MFWT are in Sverdrups. The MHT is in petawatts. The last row of the table summarizes the results, excluding rows 1 and 4. The breakdown of the results used to construct this table is given in Table A1. The “I. C.” with the υref means initial conditions.
A summary of the sensitivities of RTHIM metrics to various parameters. Each row details an ensemble of runs where a particular parameter (column 1) is varied within a range (column 2). All other parameters are kept constant while exploring a particular parameter. The remaining columns give the metrics of the rms difference between NEMO and RTHIM reference-level velocities and depth-integrated transports υrms and Tskill; and the MOC, MHT, and MFWT as described in section 4. Where range appears in the column head, the lower and upper values of the metric from the ensemble are given. Where a value appears in italics in the column head, these are the NEMO values. The parameters are the number of nodes in S and θ on the RTHIM S–θ grid, the vertical and lateral diffusivities D and K are applied to the diffusive flux initial condition, the w is applied to the section net transport constraint, the random noise is applied to the υref initial condition, and the number of cells of boxcar smoothing is applied to the diffusive flux initial condition. The MOC and MFWT are in Sverdrups. The MHT is in petawatts. The last row of the table summarizes the results, excluding rows 1 and 4. The breakdown of the results used to construct this table is given in Table A1. The “I. C.” with the υref means initial conditions.

The principle sensitivities of the RTHIM solutions are to the Sθ grid and K applied to the F initial condition. This is evident by the purple and dark blue lines in Fig. 7 and in the first five rows below the headings of Table 1. In particular a value of K = 10 m2 s−1 significantly degrades the solution (see Fig. 8, and upper values in the ranges of row 4, columns 3 and 4 of Table 1). The sensitivity to the other parameters is small; however, we find that if no smoothing is applied to F in the initial condition, then RTHIM will not converge, and if we remove the net transport constraint (i.e., set w = 0), the solution suffers dramatically (see Fig. 8). If we consider that this exercise has informed the appropriate choices of parameters as a minimum Sθ grid size of (13 × 13) and a K of 100–500 m2 s−1, from the last row of Table 1, we can say that the MHT and MFWT are determined by RTHIM with uncertainties of ±0.03 PW, and ±0.03 Sv, respectively, and their mean values agree perfectly with the NEMO model truth. In the case of the MOC, the distribution is slightly skewed, with a range of 10.76–12.91 Sv surrounding a mean of 11.64 Sv; the mean is still in good agreement with the model truth of 11.53 Sv. These uncertainties on the MOC, MHT, and MFWT metrics translate to percentage errors of less than 12%, 10% and 19%, respectively, when compared with the model truth.

Fig. 8.

Two RTHIM runs where the selection of model parameters has significantly degraded the solution. In red is the result of setting K = 10 m2 s−1 in the calculation of the F initial condition (red). In blue is the result of removing the net transport constraint. (top) The quantity υref for the initial condition (dashed lines) and solution (solid lines) for each run. (bottom) The streamfunctions in (left) density, (middle) temperature, and salinity coordinates using the same colored line notations.

Fig. 8.

Two RTHIM runs where the selection of model parameters has significantly degraded the solution. In red is the result of setting K = 10 m2 s−1 in the calculation of the F initial condition (red). In blue is the result of removing the net transport constraint. (top) The quantity υref for the initial condition (dashed lines) and solution (solid lines) for each run. (bottom) The streamfunctions in (left) density, (middle) temperature, and salinity coordinates using the same colored line notations.

A strength of this method is the ability to view the circulation in thermohaline coordinates. Figure 9 (left panel) shows Iadv as it is used within RTHIM: binned in Sθ coordinates. The small residual between the term provided by the RTHIM solution and the term constructed by binning the NEMO section transports on the same grid (right panel) shows the two are very similar. The transport through the section is dominated by a northward flow of warm, salty waters (reds) and a compensating southward flow of cooler, fresher waters (blues). This is what would be expected for this region, since warm, salty North Atlantic water flowing northward is balanced by cooler, fresher water created by boundary processes in the Arctic and subpolar gyre.

Fig. 9.

The advection term Iadv plotted in Sθ coordinates from (left) an RTHIM solution and (right) the residual with the same term based on the NEMO section velocities.

Fig. 9.

The advection term Iadv plotted in Sθ coordinates from (left) an RTHIM solution and (right) the residual with the same term based on the NEMO section velocities.

c. Volume budget

In section 2b we described the volume budget in thermohaline coordinates that RTHIM uses to obtain a solution. We can examine the terms in the volume budget to gain an insight into the contributions of each process to the increase or decrease in volume of water masses in each region of Sθ space (Fig. 10). Mixing processes (Fig. 10a) tend to increase the volume of water toward the center of the Sθ distribution while decreasing the volume of water at its fringes. This is in accordance with what we expect, since mixing acts down tracer gradients to homogenize extremes in S and θ. By contrast, surface fluxes (Fig. 10b) tend to increase the volume of water at the fringes of the Sθ distribution while decreasing the volume toward the center. This again is an expected consequence of surface forcing: in particular, we can see cooling of the warm surface waters as a result of the Arctic winds and salinification of the cooler water as a result of brine rejection in ice formation, both leading to the formation of dense waters. We have already discussed the advection term (Fig. 10c) in terms of the inflow and outflow of different water masses; what we can also see in Fig. 10 is how the residual created by adding the mixing and surface flux terms is largely accounted for by the advection term. The volume trend (Fig. 10d) is comparatively small.

Fig. 10.

Each term from the volume budget of Eq. (11), plotted in Sθ coordinates, from an RTHIM solution. (a) The mixing term, (b) the surface flux term, (c) the advection term, and (d) the volume trend with red colors indicating net positive contribution to the volume of water in that (S, θ) (S, θ) bin by a given process and blue colors a net negative contribution.

Fig. 10.

Each term from the volume budget of Eq. (11), plotted in Sθ coordinates, from an RTHIM solution. (a) The mixing term, (b) the surface flux term, (c) the advection term, and (d) the volume trend with red colors indicating net positive contribution to the volume of water in that (S, θ) (S, θ) bin by a given process and blue colors a net negative contribution.

After adding together the terms in the volume budget, we are left with a residual ε. This reduces by several orders of magnitude between the initial condition and the RTHIM solution (Fig. 11).

Fig. 11.

The residual in the volume budget in Sθ coordinates (left) from the RTHIM initial condition and (right) after optimization. Note the different color scales.

Fig. 11.

The residual in the volume budget in Sθ coordinates (left) from the RTHIM initial condition and (right) after optimization. Note the different color scales.

An alternative way to view the volume budget is by integrating in one direction or the other in thermohaline coordinates. This helps us to see more quantitatively the size of the contributions of each process to the volume budget and their relative roles at each point in the S or θ range. Integrated first through S and then cumulatively in θ (Fig. 12, left panel), we see that at low temperatures the balance is dominated by a competition between mixing and surface fluxes, whereas at higher temperatures the advection is the larger term, balancing out the other two. The y values of the lines in Fig. 12 at a given x represent the transformation in volume below that isotherm (or isohaline). This means that a range in x where the gradient dy/dx is positive corresponds to an increase in the volume of water in that temperature or salinity range by a given process and that a range with negative gradient corresponds to a decrease in volume. Therefore, in the left panel of Fig. 12 we see surface fluxes increasing the volume of the coldest water and decreasing the volume of the warmer water; meanwhile, mixing decreases the volume of the coldest (from −3° to −1.3°C) and warmest waters (>2.5°C), and increases the volume of waters at intermediate temperatures; and finally, advection decreases the volume of water up to 2.5°C (outflow) and increases the volume above that temperature (inflow). The trend term is relatively small. In salinity coordinates (right panel in Fig. 12), the balance is between mixing and surface fluxes between 29.5 and 32.5 psu, and between mixing and advection at higher salinities, with surface fluxes playing a smaller role. Mixing increases the volume between 29.5 and 31.5 psu and decreases the volume between 31.5 and 32.5 psu; meanwhile, surface fluxes have the opposite effect. Mixing then increases the volume between 32.5 and 34.5 psu and decreases the volume of the saltiest water, with both changes being balanced.

Fig. 12.

Terms in the volume budget integrated in Sθ space. (left) Each term has been integrated first through all S values and then cumulatively in θ, and plotted against θ. (right) Each term has been integrated first through all θ values and then cumulatively in S, and plotted against S.

Fig. 12.

Terms in the volume budget integrated in Sθ space. (left) Each term has been integrated first through all S values and then cumulatively in θ, and plotted against θ. (right) Each term has been integrated first through all θ values and then cumulatively in S, and plotted against S.

In section 4b, we validated one of the unknowns in our RTHIM solution, the advection term, by comparing section velocities with a known model truth from NEMO and exploring the solution sensitivities. We now examine the other unknown: the mixing term. As noted in section 2b, the term constrained in the RTHIM solution is Umix, the mixing term from Fig. 10. Using NEMO diagnostics of the tendencies in temperature and salinity resulting mixing, we can construct the components of the transformation vector resulting from mixing as follows:

 
U|Smix-NEMO(S*,θ*)=S0ΔSΔθΠ(S*)Π(θ*)dSdtmixdV.
formula
(24)
 
U|θmix-NEMO(S*,θ*)=1ΔSΔθΠ(S*)Π(θ*)dθdtmixdV,
formula
(25)

where mix/dt (K s−1) and dSmix/dt (s−1) are the NEMO temperature and salinity trends due to mixing. We can then calculate Umix-NEMO, the NEMO equivalent of our RTHIM mixing term, to compare with the RTHIM solutions.

We can compare the mixing terms from the RTHIM solution and NEMO by plotting each on the same Sθ grid (Fig. 13). There are some broadly similar patterns, such as the decrease in volume of waters at low temperatures and an increase in volume toward the center of the Sθ distribution; and there are differences, such as the larger values in NEMO at lower salinities. It is not surprising that RTHIM does not capture the details of the NEMO diffusive transformation. The mixing term is highly differentiated (it is calculated from 2F) and the dominant signals are likely the advection at the boundary and surface fluxes, with mixing being a small residual in many cases. So, a small percentage error in the flow at the boundary may manifest as a large-sized error in the mixing terms.

Fig. 13.

Mixing terms diagnosed directly from (left) NEMO tendencies and (right) an RTHIM solution .Red colors indicate a net positive contribution to the volume of water in that (S, θ) bin and blue indicate a net negative contribution.

Fig. 13.

Mixing terms diagnosed directly from (left) NEMO tendencies and (right) an RTHIM solution .Red colors indicate a net positive contribution to the volume of water in that (S, θ) bin and blue indicate a net negative contribution.

We can also integrate the mixing terms from Fig. 13 in Sθ space, in one dimension or the other (Fig. 14). Although the structure of the NEMO and RTHIM terms differs, the convergence of the blue and black lines at the top end of the ranges in each panel demonstrates that the total diffusive transformation integrated over all Sθ space is almost the same (0.42 and 0.14 Sv in RTHIM and NEMO, respectively). This shows that the bulk diffusive transformation is in the same balance as in the NEMO truth, in equilibrium with advection and surface fluxes when integrated over the control volume.

Fig. 14.

As in Fig. 12, but for mixing terms from an RTHIM solution (blue line) and NEMO (black line).

Fig. 14.

As in Fig. 12, but for mixing terms from an RTHIM solution (blue line) and NEMO (black line).

5. Discussion and conclusions

a. Summary

We have developed “RTHIM,” a new inverse method in thermohaline coordinates, to use temperature, salinity, and surface flux observations to provide estimates of water mass transformation and circulation in an enclosed domain. The method builds on the work of Groeskamp et al. (2014a,b) and others, and allows analysis previously done globally to be applied to a specific region of the ocean, thereby retaining additional useful geographical information about the circulation. Using inputs of time-mean S and θ, and surface fluxes of heat and salt, RTHIM solves for a surface reference-level velocity υref at the section bounding the domain and interior mixing within it while constraining the mixing to be down tracer gradients and υref to be within ±10 cm s−1 of an initial guess. Full section velocities are constructed from υref assuming a known vertical velocity shear. RTHIM was applied to model data from NEMO, and the solution compared to the model truth. The rms error between RTHIM and NEMO for υref was less than 0.63 cm s−1 from an ensemble of solutions with varying parameters in RTHIM. Estimates of the MOC from the ensemble were 10.76–12.91 Sv with a mean of 11.64 Sv, compared to 11.53 Sv in NEMO. The MHT and MFWT from RTHIM were 0.29 ± 0.03 PW and −0.16 ± 0.03 Sv, respectively, in exact agreement with NEMO. The transformation caused by mixing from the RTHIM solutions was also compared with that obtained diagnostically from NEMO, and while there are some similar patterns the agreement is not as strong. However, this is expected because the mixing in NEMO is governed by a parameterization including a specified horizontal eddy diffusivity, whereas in RTHIM we constrain only the diffusive flux tensor F.

b. Previous studies of the region

We have applied RTHIM to a domain including the Arctic and part of the subpolar North Atlantic, since our future aim is to use it to estimate the circulation and mixing from observations in the OSNAP study region. There have been some previous efforts to study the circulation in these regions using a thermohaline-coordinate system. In a purely diagnostic study, Pemberton et al. (2015) calculated the contributions of mixing, surface fluxes, and advection to water mass transformation in the Arctic using NEMO model outputs. They used much higher resolution in Sθ space than we were able to achieve with RTHIM, and they also separated the effects of isopycnal and diapycnal mixing; however, since their method was diagnostic, it could not be applied to observations where the complete fields are not available. In comparing our Fig. 10 to Fig. 16 in Pemberton et al. (2015), there is striking similarity between the advection terms. The surface and interior mixing terms are broadly similar, with formation of water at the fringes of the Sθ distribution by the surface fluxes with destruction toward the center, and vice versa for the interior mixing (note the opposite sign convention for the formation/destruction of water masses). However, there are some differences in the details. We see relatively little transformation of the coldest water between salinities of 29 and 32 psu in the RTHIM solution, whereas Pemberton et al. see significant transformation. On the other hand, we see significant transformation up to ~2°C in the salinity range 32–34 psu, whereas Pemberton et al. find almost all the transformation occurs below 0°C in the same salinity range. The differences are likely due to a combination of the higher resolution in the Pemberton et al. study, the different domains [Pemberton et al.’s Arctic does not include the Greenland–Iceland–Norway (GIN) Seas], and the differences between RTHIM’s and NEMO’s treatment of the interior mixing, as previously discussed. An inverse estimate based on observations was carried out by Tsubouchi et al. (2012) for one month of data in the summer of 2005. Their method relied on current data from moorings and CTD sections, and they calculated boundary fluxes and water mass transformation using a box inverse approach. RTHIM could be applied to observational reanalysis to obtain such estimates over a much longer time period.

c. RTHIM sensitivity

In section 4b, we found that RTHIM solutions were robust, provided some constraints are placed on the parameters involved in the inverse model setup. RTHIM demonstrated good skill for an Sθ grid of at least 13 × 13 cells; however, the remainder of the sensitivity study was carried out with a grid size of 17 × 17, the highest resolution for which RTHIM was able to converge. When applying RTHIM to observations, it would make sense to use the highest possible resolution so as to retain the maximum information in both Sθ space and physical space for the determination of υref. The trade-off would be one of computation time, and this would become a factor if large numbers of RTHIM runs were required. For the diapycnal diffusivity D applied to the initial condition of the diffusive flux tensor, we explored only two values: 10−4 and 10−5 m2 s−1. This is reasonable because the average diapycnal diffusivity over the global ocean is well constrained, although its value can vary locally by several orders of magnitude (see, e.g., Munk and Wunsch 1998). The volume-weighted time and spatial mean of the vertical diffusivity in NEMO for our domain was 5.6 × 10−5 m2 s−1. We explored a wider range of the isopycnal diffusivity K, since this is a more uncertain parameter involved in the parameterization of unresolved eddies in the NEMO model. Our optimal initial condition values of K = 100–500 × 10−5 are also consistent with NEMO, where the volume-weighted mean lateral diffusivity for our domain was 473 m2 s−1. The inability of RTHIM to converge on a solution when F is not smoothed in the initial condition reflects the fact that F is twice differentiated in Sθ space, so any noise on the initial fields may be amplified, causing very large initial values of SθUSθmix. A minimal amount of smoothing of F proved to be enough to overcome this problem. The sensitivity to the other parameters explored was small.

d. Future application to observations

This validation exercise has proven the effectiveness of RTHIM when applied to an internally consistent system, where the ocean has been forced by the prescribed surface fluxes and where there are no gaps in the data input to the inverse model. When applying the model to observations, there will be some additional challenges to overcome.

We can obtain the temperature and salinity fields needed to construct the diffusive fluxes’ initial condition and bin the advective fluxes in Sθ space from observational reanalysis, such as the EN4 dataset (Good et al. 2013). However, the sparseness of observations over much of the Arctic means that these regions are heavily weighted toward model data. It remains to be seen how much of an obstacle this will pose to obtaining good RTHIM solutions; however, since the interior observations are required only for constructing the initial condition, and the validation has proven some resilience to this, it seems reasonable to be optimistic.

In our validation we have also assumed the surface fluxes to be a known quantity. This will not be the case for the application, where again there is quite significant uncertainty in places because of a lack of available observations. This will contribute a degree of uncertainty to the mixing solution, which is affected through the volume budget. Our intended approach is to run RTHIM using a range of surface flux products [e.g., Lindsay et al. (2014) suggest ERA-Interim, CFSR, and MERRA, which are the most consistent with observations in the Arctic] to constrain our estimates. Recent efforts to address the lack of observations, such as the Norwegian Young Sea Ice (N-ICE2015) expedition (Granskog et al. 2016) and the upcoming Multidisciplinary Drifting Observatory for the Study of Arctic Climate (MOSAiC) expedition (http://www.mosaicobservatory.org/index.html), may also provide valuable data for the application of RTHIM to the region.

During the validation we have calculated an initial condition for the reference-level velocity by adding noise to the model truth and have taken the vertical velocity shear as the model shear. When applying to observations, we will obtain our υref initial condition from AVISO surface geostrophic velocities. We have accounted for uncertainties in the AVISO velocities through the noise added to the initial condition, but these errors can be larger at land boundaries. Gourcuff et al. (2011) found that uncertainties in altimetric surface geostrophic velocities near the oceanic western boundary with Greenland led to uncertainties in the East Greenland Irminger Current from an inverse method of ~30%. They also found that in general altimetric velocities were weaker than those obtained from ship acoustic Doppler current profiler observations along the Observatoire de la Variabilité Interannuelle et Décennale en Atlantique Nord (OVIDE) section between Greenland and Portugal. These factors will need to be considered when making use of AVISO surface geostrophic velocities with RTHIM. The assumption of known velocity shear could be met by a variety of observational constraints. For example, current meter moorings could give direct measurements of velocity. Alternatively, below the Ekman layer and away from strong currents, the flow is geostrophic, so thermal wind shear may be derived from observations of S and along the section. For the OSNAP section, both types of observational constraint are present. The region is also well sampled by the Argo program, so reanalyses are well constrained by observations (Good et al. 2013). Furthermore, for the OSNAP region, the Ekman component of the advection is weak, contributing only 1 Sv to the overturning (Mercier et al. 2015).

Finally, we have carried out the current analysis for a 20-yr mean of NEMO data. Ultimately our intention is to use RTHIM to examine temporal variability in the circulation, which would mean shortening the averaging period, perhaps to interannual or shorter variability. This could increase the size of the trend term, which in the validation proved small enough to neglect (although we have not done so); since it will be harder to calculate from uncertain observed interior distributions of S and θ, its importance will place a lower limit on the time scale of variability we can study.

e. Comparison with earlier methods

RTHIM has built on previously developed inverse methods for diagnosing the oceanic circulation from hydographic observations, and it has several advantages over its predecessors. It is a natural extension of the traditional box inverse method for determining the flow through a section by solving for a reference-level velocity, but instead of dividing the section into layers of depth or density, we use contours of S and θ, thereby defining a coordinate system suitable for studying the conservation of salt and heat, respectively. Both the Bernoulli method (Killworth 1986) and beta-spiral method (Stommel and Schott 1977) make use of gradients of tracers, such as potential vorticity on isopycnals, but assume no mixing in their solution. The tracer contour inverse method (TCIM) of Zika et al. (2010) and thermohaline inverse method of Groeskamp et al. (2014a) have mixing as a leading-order process; however, as for the box inverse, the mixing has a fixed spatial structure. In the case of TCIM and the box inverse, they assume one coefficient of isopycnal and diapycnal mixing in each layer. The ability of RTHIM to constrain mixing to be downgradient without imposing such a simplified structure makes it a major advance on earlier methods.

f. Conclusions

In summary, we have presented and carried out a successful validation of a new inverse method, RTHIM. While its application to observations will present some challenges, the method has the potential to provide independent estimates of circulation and interior mixing complementary to the OSNAP array observations, and it could ultimately be applied to any enclosed region of the ocean in future work. A verification test of RTHIM and its assumptions may be derived from a hindcast of independently observed section transport, such as from the OSNAP array. If the hindcast proves to be skillful, that further supports the likely validity of RTHIM solutions for the OSNAP section at other times.

Acknowledgments

NM, CW, and NPH acknowledge funding from the U.K. Natural Environment Research Council under the U.K. OSNAP Large Grant (NE/K010875.1). NPH is additionally supported by the U.K. Natural Environment Research Council’s North Atlantic Climate System Integrated Study program (NE/N018044/1). We would also like to thank George Nurser (NOC) and Adam Blaker (NOC) for their help with the NEMO outputs. The model output used for the inverse model validation is stored at the British Oceanographic Data Centre (under https://doi.org/10.5285/6371e0e7-8197-2959-e053-6c86abc0630d).

APPENDIX

Mixing Constraints, Grid Design, and Sensitivity Study Parameters

a. Symmetric diffusive flux tensor constraint

Here we derive the constraint that the diffusive flux tensor is symmetric as used in RTHIM. Consider temperature and salinity isosurfaces with unit normal vectors ∇θ/|∇θ | and ∇S/|∇S|, respectively. The flux of heat and salt per unit θ and S across these surfaces is given by F with components

 
Fθθ(S*,θ*)=1ΔSΠ(S*)(Kθθ/|θ|)dAθ*,
formula
(A1)
 
FθS(S*,θ*)=1ΔθΠ(θ*)(KθS/|S|)dAS*,
formula
(A2)
 
FSθ(S*,θ*)=1ΔSΠ(S*)(KSθ/|θ|)dAθ*,
formula
(A3)
 
FSS(S*,θ*)=1ΔθΠ(θ*)(KSS/|S|)dAS*,
formula
(A4)

where Π(S*) is a boxcar function over S* ± ΔS/2, K is a 3D and time-dependent diffusion tensor, and dAS* is an area integral over the surface where S = S*. Each component of F, FC1C2, represents the diffusive flux of tracer C1 across and in the direction perpendicular to the isosurface of tracer C2.

Considering specifically the diffusive flux of temperature across the salinity surface and salt across the temperature surface in the limit as Δθ and ΔS become small we have

 
Fθθ(S*,θ*)=1ΔSΔθΠ(θ*)Π(S*)KθθdV
formula
(A5)
 
FθS(S*,θ*)=1ΔSΔθΠ(θ*)Π(S*)KθSdV
formula
(A6)
 
FSθ(S*,θ*)=1ΔSΔθΠ(θ*)Π(S*)KSθdV
formula
(A7)
 
FSS(S*,θ*)=1ΔSΔθΠ(θ*)Π(S*)KSSdV.
formula
(A8)

Typically, coarse-resolution ocean models use diffusion tensors with off-diagonal terms to account for subgrid adiabatic advective effects that arise in Eulerian coordinates (e.g., Stoke’s drift and transient eddy advection). These are not relevant to transport across isosurfaces, which move with such adiabatic advection, and the diffusion tensor can be assumed to have only diagonal elements. We assume the same diffusion tensor applies to both heat and salt (i.e., that turbulent dominates over molecular diffusion). The diffusivity tensor having only diagonal elements and being equivalent for both heat and salt implies that Kθ ⋅ ∇S = KS ⋅ ∇θ and hence that F = FθS.

b. Sθ grid nodes

The grid nodes constructed from the transport functions detailed in Fig. 4, including additional nodes at either end to ensure binning of the whole range of (S, θ) found within the control volume, are as follows:

  • θ: [−4, −3.05, −1.65, −1, −0.75, −0.1, 0.55, 0.95, 1.3, 2, 3.05, 6.1, 7.2, 7.55, 8.8, 13, 15, 14]

  • S: [6, 6.93, 32.15, 32.86, 33.76, 34.12, 34.46, 34.67, 34.75, 34.78, 34.83, 34.88, 34.95, 35.01, 35.05, 35.55, 36]

c. Sensitivity study breakdown

A breakdown of the sensitivity study summarized in Table 1 is shown in Table A1.

Table A1.

Breakdown of the results of the sensitivity study summarized in Table 1. For each group of runs (column 2), all parameters (columns 3–7) are kept constant except the one being tested for the sensitivity of the metrics (columns 8–12) to that parameter. The parameters are the number of nodes on the RTHIM S and θ grids, the vertical and lateral diffusion coefficients are applied to the diffusive flux initial condition D and K (m2 s−1), the w is on the section solution net transport constraint, the random noise is applied to the υref initial condition (I.C., no column), and the smoothing is applied to the diffusive flux initial condition (I.C.). The metrics of υrms (cm s−1), Tskill (Sv), MOC (Sv), MHT (PW), and MFWT (Sv) are as described in section 4. Where metric values are reported with a dash, RTHIM failed to converge on a solution. A 5 cm s−1 random noise is added to the initial υref for all runs, but for each grouping except for the υref I.C. noise grouping, the same noise profile is used for each run within a group.

Breakdown of the results of the sensitivity study summarized in Table 1. For each group of runs (column 2), all parameters (columns 3–7) are kept constant except the one being tested for the sensitivity of the metrics (columns 8–12) to that parameter. The parameters are the number of nodes on the RTHIM S and θ grids, the vertical and lateral diffusion coefficients are applied to the diffusive flux initial condition D and K (m2 s−1), the w is on the section solution net transport constraint, the random noise is applied to the υref initial condition (I.C., no column), and the smoothing is applied to the diffusive flux initial condition (I.C.). The metrics of υrms (cm s−1), Tskill (Sv), MOC (Sv), MHT (PW), and MFWT (Sv) are as described in section 4. Where metric values are reported with a dash, RTHIM failed to converge on a solution. A 5 cm s−1 random noise is added to the initial υref for all runs, but for each grouping except for the υref I.C. noise grouping, the same noise profile is used for each run within a group.
Breakdown of the results of the sensitivity study summarized in Table 1. For each group of runs (column 2), all parameters (columns 3–7) are kept constant except the one being tested for the sensitivity of the metrics (columns 8–12) to that parameter. The parameters are the number of nodes on the RTHIM S and θ grids, the vertical and lateral diffusion coefficients are applied to the diffusive flux initial condition D and K (m2 s−1), the w is on the section solution net transport constraint, the random noise is applied to the υref initial condition (I.C., no column), and the smoothing is applied to the diffusive flux initial condition (I.C.). The metrics of υrms (cm s−1), Tskill (Sv), MOC (Sv), MHT (PW), and MFWT (Sv) are as described in section 4. Where metric values are reported with a dash, RTHIM failed to converge on a solution. A 5 cm s−1 random noise is added to the initial υref for all runs, but for each grouping except for the υref I.C. noise grouping, the same noise profile is used for each run within a group.

REFERENCES

REFERENCES
Bower
,
A. S.
,
M. S.
Lozier
,
S. F.
Gary
, and
C. W.
Böning
,
2009
:
Interior pathways of the North Atlantic meridional overturning circulation
.
Nature
,
459
,
243
247
, https://doi.org/10.1038/nature07979.
Byrd
,
R. H.
,
M. E.
Hribar
, and
J.
Nocedal
,
1998
:
An interior point algorithm for large-scale nonlinear programming
.
SIAM J. Optim.
,
9
,
877
900
, https://doi.org/10.1137/S1052623497325107.
Good
,
S. A.
,
M. J.
Martin
, and
N. A.
Rayner
,
2013
:
EN4: Quality controlled ocean temperature and salinity profiles and monthly objective analyses with uncertainty estimates
.
J. Geophys. Res. Oceans
,
118
,
6704
6716
, https://doi.org/10.1002/2013JC009067.
Gourcuff
,
C.
,
P.
Lherminier
,
H.
Mercier
, and
P. Y.
Le Traon
,
2011
:
Altimetry combined with hydrography for ocean transport estimation
.
J. Atmos. Oceanic Technol.
,
28
,
1324
1337
, https://doi.org/10.1175/2011JTECHO818.1.
Granskog
,
M.
,
P.
Assmy
,
S.
Gerland
,
G.
Spreen
,
H.
Steen
, and
L.
Smedsrud
,
2016
:
Arctic research on thin ice: Consequences of Arctic sea ice loss
.
Eos
,
97
, https://doi.org/10.1029/2016EO044097.
Groeskamp
,
S.
,
J. D.
Zika
,
T. J.
McDougall
,
B. M.
Sloyan
, and
F.
Laliberté
,
2014a
:
The representation of ocean circulation and variability in thermodynamic coordinates
.
J. Phys. Oceanogr.
,
44
,
1735
1750
, https://doi.org/10.1175/JPO-D-13-0213.1.
Groeskamp
,
S.
,
J. D.
Zika
,
B. M.
Sloyan
,
T. J.
McDougall
, and
P. C.
McIntosh
,
2014b
:
A thermohaline inverse method for estimating diathermohaline circulation and mixing
.
J. Phys. Oceanogr.
,
44
,
2681
2697
, https://doi.org/10.1175/JPO-D-14-0039.1.
Groeskamp
,
S.
,
B. M.
Sloyan
,
J. D.
Zika
, and
T. J.
McDougall
,
2017
:
Mixing inferred from an ocean climatology and surface fluxes
.
J. Phys. Oceanogr.
,
47
,
667
687
, https://doi.org/10.1175/JPO-D-16-0125.1.
Hieronymus
,
M.
,
J.
Nilsson
, and
J.
Nycander
,
2014
:
Water mass transformation in salinity–temperature space
.
J. Phys. Oceanogr.
,
44
,
2547
2568
, https://doi.org/10.1175/JPO-D-13-0257.1.
Killworth
,
P. D.
,
1986
:
A Bernoulli inverse method for determining the ocean circulation
.
J. Phys. Oceanogr.
,
16
,
2031
2051
, https://doi.org/10.1175/1520-0485(1986)016<2031:ABIMFD>2.0.CO;2.
Large
,
W. G.
, and
S. G.
Yeager
,
2004
: Diurnal to decadal global forcing for ocean and sea-models: The data sets and flux climatologies. NCAR Tech. Note NCAR/TN-460+STR, 105 pp., https://doi.org/10.5065/D6KK98Q6.
Li
,
F.
,
M. S.
Lozier
, and
W. E.
Johns
,
2017
:
Calculating the meridional volume, heat, and freshwater transports from an observing system in the subpolar North Atlantic: Observing system simulation experiment
.
J. Atmos. Oceanic Technol.
,
34
,
1483
1500
, https://doi.org/10.1175/JTECH-D-16-0247.1.
Lindsay
,
R.
,
M.
Wensnahan
,
A.
Schweiger
, and
J.
Zhang
,
2014
:
Evaluation of seven different atmospheric reanalysis products in the Arctic
.
J. Climate
,
2
,
2588
2606
, https://doi.org/10.1175/JCLI-D-13-00014.1.
Lozier
,
M. S.
,
2012
:
Overturning in the North Atlantic
.
Annu. Rev. Mar. Sci.
,
4
,
291
315
, https://doi.org/10.1146/annurev-marine-120710-100740.
Lozier
,
M. S.
, and Coauthors
,
2017
:
Overturning in the subpolar North Atlantic program: A new international ocean observing system
.
Bull. Amer. Meteor. Soc.
,
98
,
737
752
, https://doi.org/10.1175/BAMS-D-16-0057.1.
Madec
,
G.
,
2008
: NEMO ocean engine. Version 3.6, Note du Pôle de Modélisation, Institut Pierre-Simon Laplace 27, 386 pp.
Mercier
,
H.
, and Coauthors
,
2015
:
Variability of the meridional overturning circulation at the Greenland–Portugal OVIDE section from 1993 to 2010
.
Prog. Oceanogr.
,
132
,
250
261
, https://doi.org/10.1016/j.pocean.2013.11.001.
Munk
,
W.
, and
C.
Wunsch
,
1998
:
Abyssal recipes II: Energetics of tidal and wind mixing
.
Deep-Sea Res. I
,
45
,
1977
2010
, https://doi.org/10.1016/S0967-0637(98)00070-3.
Naveira Garabato
,
A. C.
,
J. D.
Zika
,
L.
Jullion
,
P. J.
Brown
,
P. R.
Holland
,
M. P.
Meredith
, and
S.
Bacon
,
2016
:
The thermodynamic balance of the Weddell Gyre
.
Geophys. Res. Lett.
,
43
,
317
325
, https://doi.org/10.1002/2015GL066658.
Pemberton
,
P.
,
J.
Nilsson
,
M.
Hieronymus
, and
H. E. M.
Meier
,
2015
:
Arctic Ocean water mass transformation in ST coordinates
.
J. Phys. Oceanogr.
,
45
,
1025
1050
, https://doi.org/10.1175/JPO-D-14-0197.1.
Stommel
,
H.
, and
F.
Schott
,
1977
:
The beta spiral and the determination of the absolute velocity field from hydrographic station data
.
Deep-Sea Res.
,
24
,
325
329
, https://doi.org/10.1016/0146-6291(77)93000-4.
Timmermann
,
R.
,
H.
Goosse
,
G.
Madec
,
T.
Fichefet
,
C.
Ethe
, and
V.
Dulière
,
2005
:
On the representation of high latitude processes in the ORCA-LIM global coupled sea ice–ocean model
.
Ocean Modell.
,
8
,
175
201
, https://doi.org/10.1016/j.ocemod.2003.12.009.
Tsubouchi
,
T.
, and Coauthors
,
2012
:
The Arctic Ocean in summer: A quasi-synoptic inverse estimate of boundary fluxes and water mass transformation
.
J. Geophys. Res.
,
117
,
C01024
, https://doi.org/10.1029/2011JC007174.
Walin
,
G.
,
1982
:
On the relation between sea-surface heat flow and thermal circulation in the ocean
.
Tellus
,
34
,
187
195
, https://doi.org/10.3402/tellusa.v34i2.10801.
Wunsch
,
C.
,
1978
:
The North Atlantic general circulation west of 50°W determined by inverse methods
.
Rev. Geophys.
,
16
,
583
620
, https://doi.org/10.1029/RG016i004p00583.
Zika
,
J. D.
,
T. J.
McDougall
, and
B. M.
Sloyan
,
2010
:
A tracer-contour inverse method for estimating ocean circulation and mixing
.
J. Phys. Oceanogr.
,
40
,
26
47
, https://doi.org/10.1175/2009JPO4208.1.
Zika
,
J. D.
,
M. H.
England
, and
W. P.
Sijp
,
2012
:
The ocean circulation in thermohaline coordinates
.
J. Phys. Oceanogr.
,
42
,
708
724
, https://doi.org/10.1175/JPO-D-11-0139.1.
Zika
,
J. D.
,
N.
Skliris
,
A. J.
George Nurser
,
S. A.
Josey
,
L.
Mudryk
,
F.
Laliberté
, and
R.
Marsh
,
2015
:
Maintenance and broadening of the ocean’s salinity distribution by the water cycle
.
J. Climate
,
28
,
9550
9560
, https://doi.org/10.1175/JCLI-D-15-0273.1.

Footnotes

© 2018 American Meteorological Society.