The Observation-Based Application of a Regional Thermohaline Inverse Method to Diagnose the Formation and Transformation of Water Masses North of the OSNAP Array from 2013 to 2015

The strength of the meridional overturning circulation (MOC) in the North Atlantic is dependent upon the formation of dense waters that occurs at high northernlatitudes. Wintertime deep convection in the Labrador and Irminger Seas forms the intermediate water mass known as Labrador Sea Water (LSW). Changes in the rate of formation and subsequent export of LSW are thought to play a role in MOC variability, but formation rates are uncertain and the link between formation and export is complex. We present the ﬁrst observation- based application of a recently developed regional thermohaline inverse method (RTHIM) to a region encompassing the Arctic and part of the North Atlantic subpolar gyre for the years 2013, 2014, and 2015. RTHIM is a novel method that can diagnose the formation and export rates of water masses such as the LSW identiﬁed by their temperature and salinity, apportioning the formation rates into contributions from surface ﬂuxes and interior mixing. We ﬁnd LSW formation rates of up to 12Sv (1Sv [ 10 6 m 3 s 2 1 ) during 2014–15, a period of strong wintertime convection, and around half that value during 2013 when convection was weak. We also show that the newly convected water is not exported directly, but instead is mixed isopycnally with warm, salty waters that have been advected into the region, before the products are then exported. RTHIM solutions for 2015 volume, heat, and freshwater transports are compared with observations from a mooring array deployed for the Overturning in the Subpolar North Atlantic Program (OSNAP) and show good agreement, lending validity to our results.


a. The subpolar North Atlantic
The meridional overturning circulation (MOC) of the ocean is characterized in the North Atlantic by a northward flow of warm, salty waters in the upper layers and a compensating southward flow of cooler, fresher waters at depth (see Fig. 1). The transformation into denser waters occurs at high latitudes, where heat is lost to the atmosphere and freshwater added, and it has been proposed that the strength of the MOC is linked to the rate of production of these dense waters (Lozier 2012). Observations suggest that rates of dense water formation are fairly constant in the Nordic Seas, but more variable in the Labrador and Irminger Seas (Smeed et al. 2014), and climate models have indicated a link between changes in Labrador and Irminger Sea convection and the MOC (e.g., Zhang 2010; Danabasoglu et al. 2012); however, direct observational evidence for this link is lacking (Lozier et al. 2017). Paleoceanographers have also linked deep convection in the Labrador Sea with the strength of the AMOC over the last 1500 years using Denotes content that is immediately available upon publication as open access.
Supplemental information related to this paper is available at the Journals Online website: https://doi.org/10.1175/JPO-D-19-0188.s1. proxies (Thornalley et al. 2018). The formation of dense waters in the Labrador and Irminger basins is therefore a subject of interest.
The Overturning in the Subpolar North Atlantic Program (OSNAP) has been continuously monitoring the flow through the section shown in Fig. 1 since August 2014. OSNAP aims to quantify the strength and variability of the MOC and heat and freshwater transports through the section, using a combination of moored instruments and glider surveys described in Lozier et al. (2017). OSNAP results reported by Lozier et al. (2019 suggested that despite the deployment coinciding with a period of strong deep convection in the Labrador Sea, the MOC variability was dominated by changes in the Irminger and Iceland basins. The present study complements the OSNAP array observations with new information about the processes that transform the waters north of the section, while also providing independent estimates of the section transports. A key test of our results is that they are validated by the OSNAP observations.

b. Labrador Sea Water
Labrador Sea Water (LSW) is an intermediate water mass formed when convection creates mixed layers of cold, fresh water as deep as 1500 m in winter (Yashayaev 2007. It is formed predominantly in the Labrador Sea but also occasionally in the Irminger Sea (Pickart et al. 2003;Fröb et al. 2016) and can be found at middepths throughout the North Atlantic north of 408N and farther south along the western boundary. It intersects with high-salinity lower Mediterranean water and mixes isopycnally to produce the upper North Atlantic Deep Water (NADW), thereby contributing to the upper cell of the MOC (Talley and McCartney 1982).
In a review of 45 years of observations, Y2007 described the LSW and the processes behind its formation and transformation over its life cycle. Over a period of years, repeated convection events in winter form a class of LSW identifiable by its temperature T and salinity S properties (see, e.g., Y2007, Figs. 6, 7), which then evolves with time. The T and S of a given LSW class depend on the conditions that led to its formation, and the atmospheric forcing both throughout the year of formation and in previous years play a role in preconditioning the ocean for convection. For example, heating and freshwater flux in summer increases surface stratification, which works against convection the following winter; on the other hand, deep convection one winter homogenizes the water column, meaning that it occurs more easily during the next. The role of preconditioning in convection was confirmed by Loder (2017, hereafter YL2017). During convection and after it finishes, the newly formed LSW mixes isopycnally with warm saline intermediate waters arriving into the Labrador Sea in the boundary current that FIG. 1. The circulation of the North Atlantic Subpolar Gyre. Warm, salty water (red arrows) flows northward, loses heat and joins with colder, fresher waters (blue arrows) flowing southward. Also marked is the OSNAP and the OVIDE section (Daniault et al. 2016). Figure from Holliday et al. (2018). flows westward around the southern tip of Greenland (Fig. 1). Gradually the LSW drains away to other parts of the North Atlantic, being replaced by lighter waters as they flow into the Labrador Sea. Much of the water also recirculates within the Labrador and Irminger Seas, steadily mixing with warmer and saltier waters, thereby maintaining a steady transfer of heat and salt into the LSW. As it spreads eastward across the Irminger Sea, modified LSW can also come into contact with Northeast Atlantic Deep Water (NEADW) spreading west. The NEADW is modified Iceland-Scotland Overflow Water (ISOW) that mixes with warm/salty older NEADW, with Denmark Strait Overflow Water (DSOW), and with the eastward spreading LSW before it is advected into the Labrador Sea (Yashayaev and Dickson 2008).
The formation rate of LSW is variable interanually and estimates of its mean also vary: Pickart and Spall (2007) compiled estimates of 1-8.5 Sv (1 Sv [ 10 6 m 3 s 21 ) from studies between 1972 and 2003; Y2007 reported an average annual formation rate of 2 Sv for 1970-95, with a higher rate of 4.5 Sv for 1987-94; Haine et al. (2008) summarized estimates of 1.3-12.7 Sv from a range of studies; and LeBel et al. (2008) suggested a long-term formation rate of 11.9 Sv based on CFC-11 inventories. There is also conflicting evidence of the link between deep convection (and subsequent LSW formation) and the rate of LSW export from the subpolar gyre. Schott et al. (2004) found similar export rates in two different time periods with significantly different amounts of convection; on the other hand Loder (2016, hereafter YL2016) found LSW export rates between 2002 and 2015 were larger in strong convection years. The main export pathway for LSW has been assumed to be in the Deep Western Boundary Current, but recent studies using real (Bower et al. 2009(Bower et al. , 2011 and simulated Zou and Lozier 2016) Lagrangian floats have found that much, if not most LSW is exported via interior pathways, and that there is significant recirculation within the subpolar gyre. In fact the exchange of LSW between the Labrador and Irminger Seas was noted in Y2007 as a ''broad wellestablished communication'' between the two basins.

c. The ocean in thermohaline coordinates
Water masses such as the LSW may be defined in various ways, for example, as water within a T-S class or density class. Transformation of water masses occurs when water moves between the relevant defined classes, via a flux of volume. The transformation may lead to water mass formation (destruction), when there is a net increase (decrease) of volume in a particular class. Figure 2 shows a volumetric distribution in T-S space for the Arctic/subpolar North Atlantic Ocean volume used in this study, with the definitions of various important water masses used in our analysis indicated by black boxes. The warm, saline water mass at the top right of the plot is the North Atlantic Water (NAW). Just cooler than the NAW are the Labrador Sea Water (LSW) and Overflow Water (OW), the latter of which includes the ISOW and DSOW; these water masses collectively occupy a small region of T-S space but considerable volume. Just saltier than these is Arctic/Atlantic Water (AAW; see, e.g., Polyakov et al. 2004;Shu et al. 2019). The bottom of the plot has the cold, fresh Arctic Surface Water (ASW) and the cold, salty Arctic Deep Water (ADW); the latter being the densest and most voluminous water mass.
There are a number of processes which transform water masses, moving them around in T-S space. These processes and their effects are summarized in Fig. 3. The tendency of surface fluxes to increase the spread of water masses in T-S space is counteracted by the tendency of mixing to bring them together. In this work, we utilize this competition to determine the relative roles of the different processes involved in water mass formation and transformation.
The regional thermohaline inverse method (RTHIM) was developed to investigate water mass transformation north of the OSNAP section, and to diagnose the relative roles of the transforming processes for each water mass. Using inputs of surface fluxes of heat and freshwater, Conservative Temperature Q, and Absolute Salinity S A along the section, and initial estimates of interior mixing and section velocities, RTHIM solves for the section flow and mixing within the control volume bounded by the section. The method was validated by Mackay et al. (2018 using model data and a control volume bounded by a simple section on the model grid. In this study we define a volume-bounding section that follows the OSNAP section and includes Bering Strait and apply RTHIM to gridded observations from 2013 to 2015. These years overlap with the OSNAP observations, allowing independent validation of the section velocities component of our solution for 2015. The period also coincides with the development of a new LSW class between 2012 and 2016 reported by YL2016, which was one of the deepest and thickest LSW layers observed since the 1990s. There are contrasts in atmospheric forcing over the three years, with 2013 seeing relatively weak convection and 2014-15 relatively strong convection in the Labrador Sea. RTHIM will enable us to diagnose the formation (or destruction) and export rates of LSW in each year, and partition that formation (destruction) into contributions from surface fluxes and interior mixing. In the results presented, rates of formation or destruction of a water mass will be referred to generically as ''formation rates,'' where negative values are associated with destruction.
Previous estimates of LSW formation rates have tended to be from inverse box models, inferred from airsea fluxes or tracer inventories, or using the thickness of an LSW layer identified using hydrographic sections as a proxy for its volume (note that YL2016 were more sophisticated, combining ship and float data to analyze the evolution of LSW week by week). By combining RTHIM with observational products that use a wealth of recent observations from remote sensing and autonomous profiling floats, we can improve upon these estimates. RTHIM combines a number of the best aspects of previous methods. It conserves volume like a box model, but also implicitly conserves heat and salt. It uses the available surface fluxes while ensuring consistency with the conservation of volume, heat and salt. And finally it apportions the water mass transformation into contributions from the surface flux and interior mixing while imposing physically realistic constraints on that mixing (see section 2 and appendix D, or M2018 for full details). We will capture the episodic nature of the LSW formation and export by analyzing three individual years of data. In addition, Y2007 suggest that their LSW production rates are likely to be underestimates as they do not account for the loss by mixing and entrainment during and shortly after convection; our solutions account for these losses.
The paper is organized as follows: in section 2 we summarize the regional thermohaline inverse method and detail how we have applied it to observation-based data. In section 3 we compare the section transports from the RTHIM solutions with observations from the OSNAP array. In section 4 we examine the formation rates for important water masses and the relative roles of the processes effecting this formation. Section 5 discusses our results and section 6 presents a summary and conclusions.

a. Method theory
The full details of RTHIM and its validation are laid out in M2018. Here we summarize the method and FIG. 2. Volumetric distribution of water masses calculated using the 2015 time mean of the EN4 objective analysis (Good et al. 2013). The colors represent the volume of water between pairs of isotherms and isohalines, summed over the ocean volume bounded by the section shown in Fig. 4 (note the logarithmic color scale). The bin sizes are DQ 5 0.18C in Conservative Temperature and DS A 5 0.02 g kg 21 in Absolute Salinity. Overlaid are boxes defining water masses in T-S space: North Atlantic Water (NAW), Arctic Atlantic Water (AAW), Arctic Deep Water (ADW), Arctic Surface Water (ASW), Labrador Sea Water (LSW), and overflow water (OW). Potential density contours s 0 are also overlaid. describe how we have applied it to observations in this study. We begin by defining an ocean control volume consisting of the Arctic and part of the NASPG, bounded to the south by the circumpolar section shown in Fig. 4a. This section is chosen because it coincides with the OSNAP array, giving us the opportunity of comparing our inverse model solutions with the OSNAP observations. We can subdivide the section into areas enclosed by isotherms and isohalines (Fig. 4b), which project into the volume and may outcrop. We then consider a volume element V between a pair of isotherms and a pair of isohalines, the volume of which will change when these isosurfaces move. The rate of volume change ›V/›t of the element is governed by the flow I adv through the section in between the isosurfaces, surface fluxes of heat and freshwater in between the same isosurfaces where they outcrop, and interior mixing across the element isosurfaces within our control volume. This volumetric balance can be expressed in thermohaline coordinates as where = 2 Su F is the mixing term, = Su Á U surf Su is the surface flux term, and « is the error in the inverse solution which we minimize. The operator = Su is the divergence operator in thermohaline coordinates. The full derivation of this equation is given in M2018, and some more details of the terms can be found in appendix D; we note here only that the diffusive flux tensor F is constrained such that it is symmetric and that mixing acts down tracer gradients. The fact that no spatial structure is imposed on the interior mixing is a unique advantage of RTHIM over other inverse methods such as the tracer contour inverse method (Zika et al. 2010), which has uniform isopycnal K and diapycnal D diffusivities on each neutral density surface, or the thermohaline inverse method (Groeskamp et al. 2014b), which has uniform K and imposes a horizontally uniform vertical structure on D.
In RTHIM, the advection term I adv is obtained similarly to a box inverse method (e.g., Wunsch 1978), by solving a reference level velocity y ref (we use the surface) and imposing a velocity shear to reconstruct the velocities for the whole section (see section 2d). The volume trend and surface flux terms are also imposed, calculated from the time-evolving Q, S A , and surface flux fields, also described in the next section. For full details of the calculation of the terms in Eq. (1), the reader is referred to M2018.

b. Datasets
We have applied RTHIM to two distinct datasets as a means of exploring the uncertainty in our inverse solutions that is due to uncertainties in the input data. The additional uncertainty due to choices of inverse model parameters is also explored (see section 3b). We apply RTHIM to each dataset for the years 2013, 2014, and 2015, giving six sets of results.
The first dataset is ECCO v4r3, a dynamically consistent time-evolving ocean state estimate with known heat and buoyancy forcing (Forget et al. 2015;Fukumori et al. 2017), which includes all of the variables required by RTHIM: surface heat and freshwater fluxes, time-evolving T and S for the control volume, and sea surface height (SSH) for the calculation of an initial geostrophic surface y ref . ECCO v4r3 has not been constrained by any observations collected for OSNAP. The second dataset uses a combination of products. The T and S fields are obtained from EN4 version 4.2.1, an objective analysis of quality controlled T and S profiles (Good et al. 2013); surface fluxes of heat and freshwater are from the ERA-Interim reanalysis (Dee et al. 2011); and SSH are from AVISO satellite altimetry. 1 In what follows, ''EN4 solution'' or ''EN4 dataset'' in the context of our RTHIM results refers to this combination of products. Of course, since ECCO uses the available observations to calculate its state estimate, the two datasets are not independent; however they are different enough to help in the exploration of uncertainty. The ECCO dataset acts as a good test bed for the application of RTHIM to real-world observations, because its time-evolving fields must obey the dynamical equations of the general circulation model on which it is based, including the conservation of volume, heat and salt consistent with its surface forcing. This strength is also a weakness however when it comes to making inferences about the real ocean from the application of RTHIM to the ECCO fields, since the ability of the model to represent the real ocean is constrained by its resolution and by any limitations in its subgrid-scale parameterizations. It is therefore expected that the model will depart from available observations to some degree in order that the fields remain dynamically consistent [see Forget et al. (2015) for an explanation of how they minimize the model-data misfit and Carton et al. 2019 for a comparison of ECCOv4r3 with temperature and salinity observations]. By contrast, we expect the combined dataset of EN4/ERA-Interim/AVISO to most closely match the available observations; however where observations are absent EN4 is relaxed to climatology, and consequently Good et al. (2013) urge caution when using the dataset to diagnose trends. Given the sparsity of observations under the polar ice cap in particular, this must be taken into consideration when interpreting our RTHIM solutions. Away from Argo float coverage, EN4 relies on hydrographic section data which are necessarily sparse in time, therefore the seasonal evolution of the T and S fields is unlikely to be resolved in these regions. In addition, while in the depth range of the Argo floats there are plentiful observations contributing to EN4, their number is considerably reduced below 2000 m. Finally there are some known issues with the EN4 fields high in the Arctic which result in some unphysical spatial distributions of T and S where it seems that sparse point observations of the North Pole Environmental Observatory (NPEO) may have been too heavily weighted (A. T. Blaker 2017, personal communication). We explore the sensitivity of our solutions to uncertainties in the EN4 dataset in sections 3 and 4.

c. T-S grid design
We have identified a number of specific water masses on which we focus our analysis, defined according to their T-S properties. Taking the volumetric T-S distribution of water masses in the control volume from each dataset, we draw boxes around each water mass in T-S space as shown in Fig. 2 for EN4 and in Fig. A1 for ECCO (see Table A1 in appendix A for water mass definitions). We are restricted by our inverse model grid to the use of rectangles in T-S space; this has some implications at the boundaries between the Labrador Sea Water and Overflow Waters that will be discussed later. The total volume of the water in all our defined water masses represents 95% of the whole control volume for EN4, and 98% for ECCO. Having established our water mass boundaries in T-S space, we then construct T-S grids that allow us to determine formation rates for those water masses by integrating terms in the RTHIM solutions over the relevant parts of T-S space. Details of these grids are in appendix B.
We have chosen a simple approach to defining water masses based on their volumetric T-S distributions for ease of comparison across different datasets when evaluating the formation rates (section 4). In the case of the LSW, the focus of this paper, there is some variation in the volumetric peak between datasets and over different years (Fig. A2). Our LSW definition includes some water that is colder than the temperature range suggested by the volumetric peaks for 2013 and 2014 but allows for a consistent definition across all three years and both datasets. We have carried out some investigations using alternative definitions of the LSW and found that regions of T-S space with low volume do not contribute significantly to the formation rates, and consequently do not impact our results. We also note that our LSW definition is a general one for the whole of our control volume, and therefore the volumetric peaks in subregions of the volume (such as the Labrador Sea, for example) are likely to be different. The use of a thermohaline coordinate system necessitates this approach but has the advantage that we make no assumptions about where in our control volume a water mass is located.

d. Further model adaptations
To build the section shown in Fig. 4a, we construct a series of subsections consisting of approximately evenly spaced points joining the lat/lon coordinates of the OSNAP mooring arrays, plus additional subsections where needed to make the section circumpolar. The dataset T and S fields are then bilinearly interpolated onto these points in the lateral plane on each depth level, giving a 2D section of T and S in the along-section/depth plane, coinciding with the OSNAP array. We then calculate relative velocities through the section using the thermal wind relation: where g is the acceleration due to gravity 5 9.81 m s 22 , f is the Coriolis parameter, r is the density of water at that point along the section calculated from the T and S fields, and z and x are the depth and distance along the section in meters. The relative velocities y geos are then added to surface reference level velocities y ref from the RTHIM solution to construct the full section velocities. The gridded EN4 and ECCO datasets have resolutions of 18 and 0.58, respectively. In order that the section definitions from the two datasets match as well as possible, we bilinearly interpolate the EN4 fields onto a 1/48 grid before they are input to RTHIM. We used a 1/48 interpolation in order that the EN4 grid points can line up with the ECCO grid points which fall every X.258 and X.758, while at the same time avoiding the loss of any information from EN4. We then apply a land mask constructed using 1/128 ETOPO5 bathymetry data 2 to both datasets. An initial condition for the surface reference level velocity y ref is calculated from geostrophy using the annual-mean SSH h: ( To reduce near-gridscale noise in the initial condition (originating mainly from the observation error of the sea surface height, the observation and representation error of the geoid, and errors in the optimal interpolation in the gridded product), we explore smoothing y ref using a boxcar filter of different widths. In the case of the AVISO h fields, we find that some smoothing is required to produce a realistic y ref , which we apply using a moving average. We explore a range of parameters for the smoothing in the RTHIM ensembles (see section 3b). The ECCO SSH anomalies are quite smooth, and so do not require this step.
To more closely match the transport calculation method in L2019 for our comparison to the OSNAP observations, we introduce an additional constraint to RTHIM. It is known from long-term observations that the net transport through Davis Strait has a long-term mean value of 21.6 Sv (Curry et al. 2014), so we add this constraint, with a weighting factor, to the net transport through the OSNAP West part of the circumpolar section. This is consistent with the zero total net transport constraint for the whole section that was used in M2018 and L2019, and we explore the sensitivity of the RTHIM solution to the weights on both constraints in section 3. We do not impose an additional constraint on the transport through Bering Strait because this would require a third weighting factor and would therefore increase the size of the parameter space to explore for sensitivity of the RTHIM solutions.

a. Qualitative comparison with observations
In this section we compare the section transports obtained from RTHIM solutions with those derived from OSNAP observations. The OSNAP velocity fields were produced using a combination of mooring data, Argo float profiles, glider data, CTD sections, and the World Ocean Atlas climatology (for full details the reader is referred to L2019). First, we have made a qualitative comparison in Fig. 5 between 2015 annual mean section velocities from RTHIM solutions and those derived from OSNAP observations. Key features such as the Labrador Current, East and West Greenland Currents, Irminger Current, a branch of the North Atlantic Current, and the southward flowing East Reykjanes Ridge Current are common between the RTHIM solutions and the observations. There is a southwardflowing current just to the west of the West Greenland Current in the EN4 solution but not the ECCO solution. Overturning streamfunctions calculated from ensembles of RTHIM solutions for 2015 using the EN4 and ECCO datasets are compared with those calculated using the OSNAP observations for the same period in Fig. 6. The ensembles contain RTHIM solutions for a range of model parameters (see section 3b and Table 1). The OSNAP section velocity, temperature and salinity fields were constructed as described in L2019. The monthly mean fields have been further averaged to obtain a 2015 mean, and the section densities calculated from the mean T and S using TEOS-10 (McDougall and Barker 2011). We use the time-mean T and S here so that the OSNAP streamfunctions are comparable to those from RTHIM, for which densities calculated from 2015 mean section T and S from each dataset have been combined with the 2015 solution section velocity. Overturning streamfunctions are calculated according to where y dA is the transport through the section (or part section) and is integrated under contours of constant density s*. For the full OSNAP section the EN4 ''best fit'' solution (see section 3b) gives the better fit of the two to the observations, and the observations fall within the envelope of the ensemble for most of the density range. The maximum of the EN4 streamfunction just above 27.7 kg m 23 is around the right density, although too large in magnitude. The ECCO streamfunction maximum is both too large and too deep, and the observed streamfunction is outside the envelope of the ECCO ensemble here. Both RTHIM solutions correctly show the majority of the overturning occurring across OSNAP East, as reported by L2019, but neither matches the observed structure across OSNAP West, with the observations going outside the ensemble envelopes for much of the density range. The largest discrepancy is the strong peak in the ECCO streamfunction for OSNAP West around FIG. 6. Overturning streamfunctions for the 2015 mean transport through the OSNAP section. Transports are integrated from lower to higher density so that the x intercept is the total transport through each section part. Individual best-fit solutions from the RTHIM ensembles using the EN4 (solid red lines) and ECCO (solid blue lines) datasets are shown along with the OSNAP observations (solid black lines). The RTHIM ensemble members are also plotted as dotted lines in their respective colors for each dataset.
27.75 kg m 23 , which can be explained by looking at the slope of the isopycnals across OSNAP West on the bottom panel of Fig. 5. The 27.7 and 27.75 kg m 23 contours are close together on the western side of the basin, where the Labrador Current flows southward, and farther apart on the eastern side where the West Greenland Current flows northward. The resulting large net northward transport in this density range gives the peak in the ECCO OSNAP West streamfunction in Fig. 6. By contrast, the same isopycnals are symmetrical across the basin in EN4 (top panel of Fig. 5) so the section transports due to the boundary current on either side of the basin in this density range largely cancel out. The density structure in EN4 is much more similar to the observations here. The observed streamfunction across OSNAP West between 27.4 and 27.9 kg m 23 is characterized by a northward flow in the lighter waters, a southward flow in the intermediate waters, and northward flow in the densest waters, resulting in a weak net overturning as measured by the maximum of the streamfunction. It may be that the overestimate in the peak from ECCO indicates an overproduction of dense water in the Labrador Sea: Li et al. (2019) suggested that such an overproduction causes a bias in the MOC in some models. We must therefore bear in mind this discrepancy when interpreting our inverse solutions based on the ECCO dataset.
The MOC derived from the full ECCO velocity fields (as opposed to the geostrophic velocity estimate we have used with RTHIM) can be seen in appendix C (Fig. C1). It is similar to the RTHIM solution using the ECCO dataset over much of the density range, but with a less pronounced spike at the density of the maximum of the streamfunction, and the maximum also appears at a lighter density. The same figure shows the adjustment that has occurred between an initial condition calculated from geostrophic surface velocities [Eq. (3)] and thermal wind shear [Eq. (2)] and the solution for the two datasets; this is most significant in the EN4 case where there is significant net transport through the section in the initial condition.

b. Quantitative comparison with observations
We now make a quantitative comparison between section transports derived from the OSNAP array observations and our RTHIM solutions. To do so we define nine metrics: the MOC, meridional heat transport (MHT), and meridional freshwater transport (MFT) for the whole OSNAP section, for OSNAP West, and for OSNAP East. The MOC is the maximum of the streamfunction in Eq. (4); the MHT and MFT are as follows: where r 0 c 0 p 5 4:1 3 10 6 J m 23 K 21 , S A0 5 35 g kg 21 is a reference salinity, and c Q and c S A are calculated analogously to c s in Eq. (4), but integrated under contours of constant Q and S A , respectively.
To determine the uncertainties on our metrics for the RTHIM solutions, we have carried out ensembles of RTHIM runs on both datasets where we explore the sensitivity of the solutions to a range of model parameters, summarized in Table 1. To test the sensitivity for each parameter, we vary one at a time while fixing the other parameters, as was done in M2018. This method TABLE 1. Summary of model parameters varied in RTHIM ensembles. Note that y ref smoothing, surface flux error, and T and S error only apply to runs using the EN4 dataset. The ERA-Interim surface flux errors were calculated by first creating arrays of random numbers between 10.5 and 20.5 with the same dimensions as the surface flux arrays. These were then multiplied by a scalar 40 W m 2 in the case of the heat flux and 8 3 10 29 m s 21 in the case of the freshwater flux, creating arrays of random errors of 6#20 W m 2 and 6#4 3 10 29 m s 21 with zero mean that were then added to their respective surface fluxes. The EN4 T and S errors were calculated in a similar way, except that the multipliers for the random number arrays were arrays with the same dimensions as the T and S fields, composed of the uncertainties on T and S included with the EN4 product. Each ensemble using the EN4 dataset contains one run where errors have been added to the surface fluxes and one where errors have been added to the T and S fields as described.

Parameter
Description Values y ref smoothing Boxcar smoothing over 2n 1 1 grid cells of y ref initial condition Integers 1-12 inclusive F smoothing 2D boxcar smoothing of diffusive flux tensor initial condition Integers 1-4 inclusive W1 Weighting factor applied to total net transport constraint 1, 10, 50, 100 W2 Weighting factor on transport through OSNAP West constraint 0.0005, 0.001, 0.005, 0.01, 0.05, 0.1 D Vertical diffusivity used in calculating F initial condition 10 24 , 10 25 m 2 s 21 K Lateral diffusivity used in calculating F initial condition 100, 200, 500 m 2 s 21 T-S grid size Number of nodes in the T and S ranges constructing the grid 15 3 15, 16 3 16, 17 3 17 Trend term Inclusion of the trend term in the volume budget Yes/no Surface flux error Random errors within bounds added to surface fluxes Heat 6#20 W m 2 , freshwater 6#4 3 10 29 m s 21 T and S error Random errors within bounds added to interior T and S Uncertainty limits from EN4 dataset is a compromise between the desire to examine fully the uncertainty on our solutions and the available computation time, since there are too many permutations for an exploration of the full parameter space to be feasible. For each ensemble, the result presented is an individual ensemble member chosen to most closely resemble the ensemble mean in terms of our metrics. We do this rather than presenting the ensemble mean of the metrics because individual solutions obey the balance based on volume, heat and salt conservation from Eq.
(1), and as such are more physically realistic than the ensemble mean. The ensemble member most closely resembling the ensemble mean is established using a function C: The sum is over our nine transport metrics TP i ; TP and TP are the individual solution and ensemble mean transport metrics, respectively; and d TP is the ensemble standard deviation. The ensemble member taken as having the best fit to the ensemble mean is that which gives the smallest value of C in Eq. (6). The best-fit solutions from the 2015 EN4 and ECCO RTHIM ensembles are shown by the colored bars in Fig. 7, along with their observational equivalents. The RTHIM solutions' error bounds show the ensemble range for each metric; note that the colored bars do not fall in the center of the error bounds because they correspond to an individual solution rather than the ensemble mean. The observational error bounds are calculated according to DMOC 5 d MOC = ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 365/(2 3 16 days) p (and similarly for DMHT and DMFT), where d MOC is the standard deviation of the MOCs calculated from the 12 individual months of observations, and 16 days is the integral time scale calculated from the autocorrelation function of the daily MOC time series (see section 1e of L2019's supplementary material). This is reported in L2019 as being a close estimate to the observational error they obtain from a more sophisticated Monte Carlo technique.
The agreement between the transport metrics from the RTHIM ensembles and the observations is generally good, with the MOC across OSNAP West from the ECCO runs the notable exception (Fig. 7). The MHT and MFT from both ensembles mostly agree with the observations within their uncertainties, with the exception being the OSNAP West MFT where both slightly underestimate the (southward) freshwater transport. The reasonable agreement between the RTHIM solution overall section transports and the observations gives us some confidence in the inferred formation rates presented in the next section. We can also diagnose the same metrics for the other two years of RTHIM ensembles as we have done for 2015; these are summarized in Table 2.

Formation rates
We now examine the volume budgets from our RTHIM solutions; this is the part of the solution from which we calculate the transformations of individual water masses, and the contributions to these transformations from different physical processes. The volume budget for an RTHIM solution from the EN4 2015 dataset is shown in Fig. 8. The solution is the ensemble member that best fits the ensemble mean, identified as described in section 3b. In the surface flux term, we see the formation of the Arctic Water by cooling around 08C and the destruction of the warm, salty North Atlantic Water. There is also a formation signal in the warmer waters between 30 and 34 g kg 21 , which falls in a region of T-S space with little to no volume on the time average (Fig. 2). Examining the monthly surface T-S distribution (see animation in the online supplemental material) reveals that the isotherms and isohalines bounding this region of T-S space outcrop in the summer months, which is responsible for this signal in the surface flux term. These waters are then mixed and transformed into other waters at the same average rate as they are  formed, as seen in the mixing term. The mixing term is in competition with the surface flux term, having the opposite sign over much of T-S space, with the strongest formation around 08C and around 38C, 34.5 g kg 21 .
In the advection term we can see the northward flow of warm salty waters and southward flow of cooler, fresher waters. The volume trend is generally smaller than the other terms, with some net formation in the main water masses, and destruction of the saltiest Arctic Water. Figure 9 shows the volume budget for the EN4 2015 solution, zoomed in to show the LSW and OW water masses. There is a formation signal in the surface flux just cooler and fresher than the box defining the LSW, and at the same density. There is also a northward flow in the advection term just warmer and saltier than the LSW, which presumably comes from the Irminger Current as it enters the Labrador Sea through OSNAP West (Fig. 1). Meanwhile, the mixing term shows destruction of these adjacent water masses while showing a net formation within the LSW box itself. This suggests that the LSW has been transformed through mixing, predominantly along isopycnals, of the cold, freshwater mass formed by deep convection and the warm, salty water mass advected in. The overflow waters also are formed predominantly by mixing and then advected southward, with the bulk of the formation in the colder OW box. The water masses mixing to form the OW are likely in the two adjacent blue areas (blue representing net destruction): the AAW slightly saltier and around FIG. 8. Volume budget in T-S coordinates for the EN4 2015 RTHIM solution. The four panels show the four terms described by Eq. (1), with regions of T-S space where there is net formation by a given process shown in red, and regions where there is net destruction shown in blue. The boundaries of our important water masses as defined in Fig. 2 are shown by the black boxes, with water mass labels shown on the volume trend term (NAW 5 North Atlantic Water, AAW 5 Arctic/Atlantic Water, ADW 5 Arctic Deep Water, ASW 5 Arctic Surface Water). Gray contours are s 0 potential density surfaces.
the same temperature, and the unlabeled water mass fresher and around 28C. We now compare the formation rates from RTHIM solutions applied to our six datasets: EN4 2013/14/15 and ECCO 2013/14/15 (Fig. 10). We have integrated the formation rates from the RTHIM solutions over the water mass definitions from Fig. 2, grouped as follows: NAW, LSW, OW, and Arctic Water (AW 5 ASW 1 ADW 1 AAW). For each of the six ensembles the colored bars show the solution that most closely resembles the ensemble mean as established using Eq. (6).
NAW is advected into the region and then transformed by a combination of mixing and surface flux, with mixing the larger term. The advection is consistent between the EN4 and ECCO solutions and constant over the three years. The picture is more complex for the mixing and surface flux terms, as mixing in EN4 reduces between 2013 (weak convection) and 2014-15 (strong convection), whereas ECCO shows mixing increasing and then decreasing while the surface flux increases steadily. The difference in this pattern is explained by the differences in the trend term, which is diagnosed directly from the T and S fields and must also be the sum of the other three terms. In the EN4 solution for 2013 the trend term is negative, meaning that the net effect of the inflow and transformation by mixing and surface flux is to reduce the volume of NAW within the region. In ECCO for 2013 the trend is positive, that is, the volume of NAW is increasing. In 2014 the opposite is true: the ECCO solution shows decreasing NAW volume while in EN4 there is a small increase. To close the volume budget in each case, RTHIM has attributed differing amounts of transformation to mixing in each solution for each year.
In the LSW there is larger variation between the solutions from the different datasets. Part of this may be attributable to the large volumes of water in the small surface flux. These waters are then exported (negative advection term), and the ECCO solution shows a decrease in volume (negative trend term). As the years progress both sets of solutions show a steady increase in the rate of production of LSW by mixing, but the export rates follow an upward trend in the EN4 solutions versus a downward trend in ECCO. The surface flux and trend terms are very small for the EN4 solutions but play a significant role in the balance for ECCO; in the latter the increased production of LSW in 2014-15 results in storage (positive trend term), rather than export. Both the EN4 and ECCO solutions show a net destruction of the LSW by the surface flux in 2015, our strongest deep convection year. This perhaps surprising result can be understood by referring back to Fig. 9: the surface fluxes are cooling our LSW and some fresher waters adjacent to it while forming a cooler and fresher water mass, which meanwhile gets mixed back to form the LSW.
In the OW, both sets of solutions show similar rates of export in all three years, although the rates are higher in ECCO. In five out of the six ensembles mixing forms the OW, the exception being the ECCO solution in 2015 where the surface flux has a major role. A closer inspection of the volume budget for the ECCO 2015 RTHIM solution (not shown) reveals that the large formation signal is at the warm boundary of the cooler OW box, and likely represents formation due to cooling by convection in the Labrador Sea. This signal appears in the EN4 solution as discussed above, but in fresher waters that are outside our OW definition.
In the AW we see an even balance between formation by surface flux and destruction by mixing, with almost no signal in the advection since these water masses do not generally flow through the section. The EN4 and ECCO solutions disagree on the magnitudes of the rates of formation/destruction of AW, and neither show an obvious trend over the 3 years. The effect of random errors on the surface fluxes has been accounted for in our EN4 ensembles, but if there is a bias in the surface flux for either dataset then this would explain the disagreement. For example, if the magnitude of the (negative) heat flux from ERA-Interim were increased this would increase the rate of production of AW in the surface flux term, and RTHIM would increase the rate of destruction by mixing to maintain the volumetric balance, reducing the disagreement with the ECCO solution.

a. Dataset comparison
We have presented RTHIM solutions based on two distinct datasets: one composed of outputs from the ECCO state estimate and one using a combination of EN4 and ERA-Interim reanalyses and AVISO satellite altimetry. Both datasets produced solutions for the advection through the OSNAP section which agree reasonably well with observations from the OSNAP array, but the EN4 dataset was the closer of the two. On examining the formation rates, the two datasets generally agree on the sign, but not the magnitude, of contributions to the volume budget from advection, mixing, surface fluxes and the volume trend. Where there is disagreement between the solutions derived from each dataset, it is difficult to know which is more realistic, since EN4 and ECCO include different types of observational influence. Where observations are always missing, EN4's persistence-based forecast form of objective analysis adjusts the solution toward climatology Good et al. (2013). However, ECCO's 4DVar state estimate is effectively a free-running forward model solution with forcing/parameters chosen so that its state most closely fits the observations. This means that information from near-surface observations can, in principle, continue to propagate into the deep ocean in the case of ECCO, provided there are dynamical connections in the forward model (Forget et al. 2015). We also look at how the formation rates differ between 2013 and 2015, three years with varying rates of convection in the Labrador Sea. Both datasets show an increase in the production rate of LSW by mixing as convection increased between 2013 and 2015, but EN4 and ECCO show opposite trends in export rates, while the ECCO solutions also show a larger role for surface fluxes and the volume trend in the balance.
In the high Arctic, the interior T and S of our control volume are poorly constrained by observations, leading to issues with EN4 described in section 2b. This affects two aspects of RTHIM: the initial condition for the mixing term which is based on T-S gradients, and the fixed trend term which is calculated from changes in the volumetric distribution of T and S. We have explored the effects of variations in the mixing term initial condition through our ensembles, and their ranges are plotted in the error bars in Figs. 7 and 10. For the EN4 ensembles these ranges include runs where we have added random errors to the T and S fields. We also briefly explore the uncertainty on the trend term in both datasets by including an ensemble member where it is set to zero (Table 1). This can be seen in Fig. 10 in the fact that the trend term uncertainties all have zero as either an upper or lower bound.
There are differences in the magnitudes of the formation rates due to surface flux for the two datasets, and with only two to compare it is difficult to be confident in a preference for one over the other. In most cases the differences are not large enough to change the nature of the volumetric balance, but there are a few exceptions in the Labrador Sea Water and Overflow Waters. It is the large volumes associated with these water masses (combined they occupy 12.5% and 12.2% of the total control volume for the EN4 and ECCO datasets, respectively) and their close proximity in T-S space which make them difficult to distinguish in this coordinate system.

b. Labrador Sea Water
The period of 2013-15 addressed in this study coincides with the development of a class of LSW defined in YL2017 as LSW 2012-16 . On their Fig. 2 this water mass has potential temperature around 3.28-3.78C and practical salinity around 34.80-34.90 psu. In our coordinates of (Q, S A ) the temperature is equivalent to our degree of precision, and the salinity converts to 34.97-35.07 g kg 21 . Our LSW definition based on the ECCO and EN4 2015 volumetric T-S distributions of Q 5 3.38-48C and S A 5 35-35.1 g kg 21 is slightly warmer and slightly saltier than this, but fits with the traditional LSW density range of s 0 5 27.68-27.80 kg m 3 quoted by Li et al. (2019) as shown in our Fig. 9. Perhaps the upper end of our LSW salinity range is a little too high and includes some of the OW for the EN4 dataset; the position of the boundary was a compromise reached to try to keep a consistent definition for both EN4 and ECCO. We may also have excluded some LSW with our lower temperature bound of 3.38C: most of the LSW in the YL2017 figure is in the temperature range 3.28-3.48C. It is possible that the warmer waters we have identified from the volumetric T-S distributions of our Figs. 2 and 6 are legitimate LSW, either formed in the Irminger Sea or formed in the Labrador Sea and recirculated. Equally it is possible that our boundary of 48C between the LSW and the NAW is a little too high. It is also worth noting that the sections from YL2017 on which these comparisons are based are snapshots from a survey done in May 2016, whereas our definitions are constructed using a time-mean of the 2015 volumetric T-S distribution.
We can interpret our results in the context of forcing over the seasonal cycle. The formation rates diagnosed due to each process are annual means, but the formation itself is likely to have taken place over shorter time periods and at different times of the year. For example, the formation of the water mass just cooler and fresher than the LSW in Fig. 9 will have occurred during convection in the winter months; meanwhile the mixing forming the LSW itself probably began during convection but continued for some time after. The isopycnal mixing in the Labrador Sea of newly convected cold, freshwater with warm, salty water brought in by advection is consistent with the description by Y2007 of the evolution of LSW which we introduced in section 1b. We also see little export of the newly convected water in the advection term of Fig. 9; instead the waters are mixed into the LSW box before being exported. These findings fit with those of Pickart and Spall (2007), who suggest that LSW is generated at the boundaries of the Labrador Sea via adiabatic eddies, and of Georgiou et al. (2019), who show that dense water formed in the interior of the Labrador Sea is laterally advected into the boundary current by eddies before it can be exported.
The steady warming, salinification and recirculation of LSW over the years following its initial formation described by Y2007 may also explain the fact that much of the water identified in our volumetric T-S distributions as LSW is somewhat warmer and saltier than that seen in the sections of YL2016 and YL2017. It is likely that much of what we see in the volumetric census of 2015 is recirculated recently formed LSW that has had more time to mix with other waters. The suggestion of Yashayaev et al. (2007) that LSW spreads to the Irminger Sea 1-2 years after its formation is consistent with this idea.
Our RTHIM solutions using the EN4 dataset give LSW formation rates from a combination of mixing and surface fluxes of 6.2, 8.3, and 11.3 Sv for 20136.2, 8.3, and 11.3 Sv for , 20146.2, 8.3, and 11.3 Sv for , and 2015, respectively, with mixing dominating. Using the ECCO dataset the formation rates are 5.6, 11.9, and 6.0 Sv, with the surface fluxes making a significant contribution to the formation; in particular the effect of the surface flux is to contribute to the formation of LSW in 2014 but to counteract it in 2015 by cooling waters that are already in the LSW class. These formation rates are in the same ballpark as the historical estimates discussed in section 1b and are also of the same order as the export rates of 3.2 and 8.9 Sv in weak convection and strong convection years, respectively, reported by YL2016. Our results may in fact be closer to the real formation rates as we have explicitly accounted for the contributions of surface fluxes and mixing to formation. As discussed above, the upper bound on our temperature range for defining the LSW may be a little high. If it is reduced to 3.78C, RTHIM solutions using the EN4 dataset give slightly lower formation rates of 4.2, 3.6, and 9.0 Sv for 2013, 2014, and 2015, respectively; the latter value in close agreement with YL2016 for strong convection years.
The question of the link between formation and export rates [i.e., respectively 2= Su Á U surf Su 2 = 2 Su F and 2I adv from Eq. (1)] remains unresolved due to the differences between our solutions from EN4 and ECCO. However, the ECCO solution offers an illustration of the disconnect because while the formation rates were higher in the strong convection years, the export rates as seen in the advection term were lower. The difference is taken up in the trend term: the volume of LSW increased in 2014 following convection (positive trend term in Fig. 10), with a further increase in 2015 partially counteracted by convection creating more source waters for the formation of new LSW by mixing (negative surface flux term, large positive mixing term). YL2016 reported that during the Argo era the average winter LSW volume was about 70% larger in strong convection years than in weak ones; however the reduction in volume from winter to autumn was 180% larger, giving a factor of 2.8 difference in potential LSW export rates in strong convection years. It will be interesting to see what is revealed by the next set of OSNAP observations in the context of the recently observed deep convection. The previously described observation that dense water formation rate is less variable in the Nordic Seas than in the Labrador and Irminger Seas is supported by our results: the AW and OW formation rates have ranges of ;2.5 and ;3.5 Sv, respectively, while the LSW formation rate has a range of ;4.5 Sv.
It is also difficult to draw firm conclusions from this study about the role of the LSW in the MOC variability. On the one hand, the MOC across the whole OSNAP section was larger in 2014-15 when we had higher rates of LSW formation; however, the relationship between the whole section MOC and that diagnosed across its western and eastern parts is unclear (see Table 2). If we consider only the EN4 solutions on the grounds of the large discrepancy between the observed MOC across OSNAP West and that derived from the ECCO solutions, the full-section MOC increases steadily while the contributions from both parts of the sections fluctuate. It may be necessary to analyze more years of data in order to establish a possible link between LSW formation and the MOC.

Conclusions
We have applied a regional thermohaline inverse method (RTHIM) to diagnose the water mass transformation in an enclosed region of the Arctic and Subpolar North Atlantic Ocean. Six sets of results were obtained by applying RTHIM to three separate years of data from two different datasets; the year 2013 where the convection in the Labrador Sea was relatively weak and 2014-15 where it was strong. For each solution we obtain transports through the section bounding the volume, and formation rates due to water mass transformation within the volume. Comparisons between inverse solution section transports and those derived using independent observations from the OSNAP mooring array were good, giving confidence in the formation rates. The latter were summarized in terms of the contributions of different processes to the formation of important water masses, with a particular focus on the LSW, and the results from the three years and two datasets compared.
Annual mean formation rates for LSW ranged from a low of ;6 Sv in 2013 when convection was weak to highs of ;12 Sv in either 2014 or 2015 (depending on which dataset was used) when convection was strong, with interior mixing playing a leading role in the formation. The effect of winter convection was to create a water mass slightly colder and fresher than the resident LSW class, but this water was not exported directly. Instead the newly convected water was mixed isopycnally with warm, salty waters carried in by advection. The product, the intermediate temperature and salinity LSW, was then exported.
This was the first application of the recently validated regional thermohaline inverse method to observationbased data. Its success in diagnosing a MOC for the time period coinciding with OSNAP that is consistent with observations indicates its potential for further analysis of the circulation in the region. By applying RTHIM to other years of the observational data products used in this study, or to similar products, yet further context could be provided for the OSNAP observations. It may be possible to look at interannual variability in the meridional overturning circulation as measured across the OSNAP section and explain this in terms of changes in the water mass transformation in the region.
volumetric distributions in the region of T-S space occupied by the LSW for the EN4 and ECCO datasets for 2013, 2014, and 2015. Table A1 shows the water mass definitions that are plotted in Figs. 2 and A1.

T-S Grid Details
The T-S grids used in this study are each defined using two row vectors: one in the T dimension and one in the S dimension. These row vectors are detailed in Table B1. The grid points in each vector define the midpoints of T-S bins for which the terms in Eq. (1) are calculated. The grids are designed such that T-S bin boundaries, which are at the midpoints between row vector points, correspond to defined water mass boundaries. The vectors cover the full range of T and S within the control volume for their respective datasets, plus an additional ''halo'' of grid points, which is required due to the inverse model discretization. The water mass boundaries plotted in Figs. 2 and A1 used the 15 3 15 EN4 2015 and ECCO 2015 grids, respectively.

APPENDIX C
RTHIM Transport Adjustment Figure C1 shows overturning streamfunctions for RTHIM solutions from ECCO and EN4 compared to their initial conditions.  Fig. 2, but using the 2015 time mean of the ECCO state estimate (Forget et al. 2015;Fukumori et al. 2017). Note that the water mass boundaries differ from those in Fig. 2 due to the different T-S grid (see appendix B).

Details of the Inverse Method
Here we describe some additional details of Eq. (1) and how it is solved to obtain our RTHIM solutions [for a full explanation including the derivation of Eq. (1) the reader is referred to M2018]. The first term, I adv , describes the transport through the section surrounding our control volume (the section shown as a red line in Fig. 4a). In each (S, u) bin of our RTHIM T-S grid, I adv (S, u) is the transport per unit S and u perpendicular to the section (positive northward, into the control volume) between the pairs of isohalines and isotherms defining that bin. The transport is the section velocity y 5 y ref 1 y geos (see section 2d) integrated over the area between the isohalines and isotherms. The I adv term, as the other terms in Eq. (1), has units of m 3 s 21 g 21 kg K 21 .
The mixing term = 2 Su F is the formation rate per unit S and u within the control volume due to mixing. It is calculated by applying the operator = 2 Su to the diffusive flux tensor F, where = Su [] 5 (›/›S, ›/›u)[] is the divergence operator in thermohaline coordinates. The diffusive flux tensor has four components (so that when applying the divergence operator twice we obtain a scalar): F uu , F Su , F uS , and F SS , where F C 1 C 2 represents the diffusive flux of tracer C 1 across and in the direction perpendicular to the isosurface of tracer C 2 . The RTHIM solution is constrained such that F SS , F uu $ 0, and F Su 5 F uS .
The volume trend term ›V/›t(S, u) is the rate of change of the volume of water in each (S, u) bin, that is, that contained within pairs of isohalines and isotherms defining each bin. This is divided by the width of the bins in T-S space, DSDu, so that it has the same units of m 3 s 21 g 21 kg K 21 . It is calculated by taking a volumetric census of the water masses at the start and end of an averaging period (e.g., the year 2015) using the T and S data from our datasets (EN4 or ECCO). The surface flux term = Su Á U surf Su (S, u) is the divergence in thermohaline coordinates of the vector U surf Su (S, u), which has components of [U surf S (S, u), U surf u (S, u)]. The components are calculated by integrating the surface fluxes of heat and freshwater from our datasets between the isohalines and isotherms defining each (S, u) bin where they outcrop.
In the RTHIM inverse calculation, the volume trend and surface flux terms are prescribed, and the advection and mixing terms are solved for. In the case of the mixing term, we prescribe the relative velocities y geos and solve for the surface reference velocity y ref , using Eq. (3) as the initial condition. In the case of the mixing term, we calculate an initial condition for the diffusive flux tensor F from gradients in our dataset temperature and salinity fields, and solve for the components of F given the constraints outlined above. During the optimization used to obtain our RTHIM solution, the sum over all T-S space of («DuDS) 2 is minimized, and this sum has final values of ,10 24 Sv 2 .