## 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.

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 × 10^{9} m^{2} 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 × 10^{16} m^{3} is contained within 10% of the occupied bins.

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 *S*_{A} 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 *I*_{adv} = *I*_{adv}(*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:

where ∂(*V*)/∂*t* is the volume tendency per unit (*S*, *θ*), ∇_{Sθ} ⋅ $US\theta dia$ is the diathermohaline volume transport divergence per unit (*S*, *θ*), *ρ*_{0} = 1024 kg m^{−3} is a reference density of seawater, *f*_{m} is the surface boundary mass flux per unit (*S*, *θ*), and *I*_{adv} is the section volume advection per unit (*S*, *θ*);

where ∂(*VS*)/∂*t* is the salt tendency per unit (*S*, *θ*), ∇_{Sθ} ⋅ $\u2061(SUS\theta dia)$ is the diathermohaline salt transport divergence per unit (*S*, *θ*), *m*_{S} is the salt convergence as a result of interior diffusion per unit (*S*, *θ*), and *I*_{adv}*S* is the section salt advection per unit (*S*, *θ*); and

where ∂(*Vθ*)/∂*t* is the heat tendency per unit (*S*, *θ*), ∇_{Sθ} ⋅ $\u2061(\theta US\theta dia)$ is the diathermohaline heat transport divergence per unit (*S*, *θ*), $cp0$ is the specific heat capacity of seawater (4000 J kg^{−1} K^{−1}), *f*_{h} 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 *I*_{adv}*θ* is the section heat advection per unit (*S*, *θ*). Note that the superscript dia denotes a transport through a surface of constant temperature or salinity.

We begin the separation into inputs and outputs by identifying the component terms of the time-mean diathermohaline volume transport $US\theta dia=(U|Sdia,U|\theta dia)$. Here $U|Sdia$ is the transport across an *S* surface between two *θ* surfaces per unit *θ*; and $U|\theta 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θ} ⋅ [⋅] = (∂/∂*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 *I*_{adv} 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):

where

are due to surface processes;

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

are related to the local response, $lS\xaf=\rho 0V\u2202S/\u2202t\xaf$ and $l\theta \xaf=\rho 0cp0V\u2202\theta /\u2202t\xaf$. 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):

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 *f*_{m}, 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\theta mix$, as the divergence of a tensor **F**, whose properties we may constrain to ensure net downgradient mixing:

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,

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 *I*_{adv} and the elements $FC1C2$ of **F** using constrained optimization (MATLAB’s fmincon), with constraints *F*_{SS}, *F*_{θθ} ≥ 0 (net downgradient mixing), and *F*_{Sθ} = *F*_{θS} (symmetric diffusive flux tensor). Note that because Eq. (11) involves the $\u2207S\theta 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 $\u2207S\theta 2N=0$ and **F**′ is the modified **F** —yet to still obtain the same volume balance: $\u2207S\theta 2F=\u2207S\theta 2F\u2032$. 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 (Δ*x*^{2} + Δ*y*^{2}). Therefore, we emphasize that the interpretation of the RTHIM solution should be focused solely on the role of $\u2207S\theta \u22c5US\theta 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:

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.

### c. Advective flux

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

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

(and similarly for *I*_{ref}), where $\upsilon rel\u2061(x,z)=\upsilon \u2061(x,z)\u2212\upsilon surf\u2061(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 *I*_{ref} = **A***υ*_{ref}, where *I*_{ref} and *υ*_{ref} are (*nS* × *n**θ*)- and (*nx*)-element column vectors, respectively, and where *nS*, *nθ*, and *nx* are the number of bins used to discretise the *S*, *θ*, and *x* grids. Then *I*_{adv}(*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 $\u2207S\theta \u22c5US\theta surf$ for input to RTHIM. The contribution of the surface fluxes to the volume budget, $\u2207S\theta \u22c5US\theta 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, Δ*z*_{surf}, obtaining *f*_{w} (s^{−1}), and $U|Ssurf$ is then calculated as follows:

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 *I*_{adv} ensures that ∇_{Sθ} ⋅ $US\theta surf$ will coincide with *I*_{adv} on the discrete grid, since $U|Ssurf$ is calculated between *S* nodes and on *θ* nodes, while $U|\theta surf$ (below) is calculated on *S* nodes and between *θ* nodes; therefore, $\u2207S\theta \u22c5US\theta surf$ = (∂/∂*S*, ∂/∂*θ*) ⋅ $\u2061(U|Ssurf,U|\theta surf)$ will be on *S* and *θ* nodes, as is *I*_{adv}. The heat flux (W m^{−2}) (Fig. 1, right panel) is divided by $\u2061(\rho 0cp0\Delta zsurf=4.1\xd7106\Delta zsurf)$ to obtain *f*_{θ} (K s^{−1}), and $U|\theta surf$ is calculated as follows:

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 $\u2207S\theta \u22c5US\theta 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 *r*^{2} (in Sv^{2}; 1 SV ≡ 10^{6} m^{3} s^{–1}) in the following:

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 *I*_{adv} and the known 1988–2007 time-mean NEMO net section transport including the Bering Strait (*T*_{net} = 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, *T*_{net} 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 *r*^{2} drops below 10^{−4} Sv^{2}, 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:

where *T*_{sol}(*x*, *z*) = *υ*_{sol}(*x*, *z*)*dA*(*x*, *z*) and *T*_{mod}(*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 *T*_{skill} 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, *T*_{skill} = 0.03 Sv, which is less than 4% of the rms of the NEMO depth-integrated transports along the section.

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 *r*^{2} = 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, *r*^{2} has reduced from ~10^{6} 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 ~10^{3} 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).

### 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:

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:

where in this case *T*(*σ*) represents the transport integrated between density contours 0.001 kg m^{−3} apart, $\rho 0cp0=4.1\xd7106Jm\u22123K\u22121$, and *S*_{0} = 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 m^{2} s^{−1}; and *D* values of 10^{−5} and 10^{−4} m^{2} 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\u2061(S,\theta )$ 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.

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 m^{2} 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 m^{2} 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.

A strength of this method is the ability to view the circulation in thermohaline coordinates. Figure 9 (left panel) shows *I*_{adv} 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.

### 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.

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).

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.

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 $\u2207S\theta \u22c5US\theta mix$, 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:

where *dθ*^{mix}/*dt* (K s^{−1}) and *dS*^{mix}/*dt* (s^{−1}) are the NEMO temperature and salinity trends due to mixing. We can then calculate $\u2207S\theta \u22c5US\theta mix-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 $\u2207S\theta 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.

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.

## 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} m^{2} 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} m^{2} 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 m^{2} 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 $\u2207S\theta \u22c5US\theta 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

where Π(*S**) is a boxcar function over *S** ± Δ*S*/2, **K** is a 3D and time-dependent diffusion tensor, and $\u222b\u222bdAS*$ is an area integral over the surface where *S* = *S**. Each component of **F**, $FC1C2$, represents the diffusive flux of tracer *C*_{1} across and in the direction perpendicular to the isosurface of tracer *C*_{2}.

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

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* = **K**∇*S* ⋅ ∇*θ* and hence that *F*_{Sθ} = *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

## REFERENCES

## Footnotes

© 2018 American Meteorological Society.

This article is licensed under a Creative Commons Attribution 4.0 license (http://creativecommons.org/licenses/by/4.0/).