On the observability of turbulent transport rates by Argo: supporting evidence from an inversion experiment

. Although estimation of turbulent transport parameters using inverse methods is not new, there is little evaluation of the method in the literature. Here, it is shown that extended observation of the broad-scale hydrography by Argo provides a path to improved estimates of regional turbulent transport rates. Results from a 20-year ocean state estimate produced with the ECCO v4 (Estimating the Circulation and Climate of the Ocean, version 4) non-linear inverse modeling framework provide supporting evidence. Turbulent transport parameter maps are estimated under the constraints of ﬁtting the extensive collection of Argo proﬁles collected through 2011. The adjusted parameters dramatically reduce misﬁts to in situ proﬁles as compared with earlier ECCO solutions. They also yield a clear reduction in the model drift away from observations over multi-century-long simulations, both for assimilated variables (temperature and salinity) and independent variables (biogeochemical tracers). Despite the minimal constraints imposed speciﬁcally on the estimated parameters, their geography is physically plausible and ex-hibits close connections with the upper-ocean stratiﬁcation as observed by Argo. The estimated parameter adjustments furthermore have ﬁrst-order impacts on upper-ocean stratiﬁcation and mixed layer depths over 20 years. These results identify the constraint of ﬁtting Argo proﬁles as an


Introduction
Direct observational estimates of vertical and lateral turbulent transport rates are largely limited to studies of surface drifter dispersion (e.g., Krauß and Böning, 1987), surface eddy fluxes estimated from satellite data (e.g., Abernathey and Marshall, 2013), occasionally released tracer dispersion (e.g., Ledwell et al., 1993), and rare micro-structure measurements (see Waterhouse et al., 2014). However, the vast collection of in situ profiles provided by Argo (Roemmich et al., 1999(Roemmich et al., , 2009 may offer new opportunities to infer turbulent transport rates from the sea surface to 2000 m depth. Inferences of diffusivities is possible through the analysis of Argo variance fields combined with conceptual turbulence models (Wu et al., 2011;Whalen et al., 2012;Cole et al., 2015). The extensive observation of the broad-scale hydrography characteristics by Argo may also provide a basis for the inversion of turbulent transport rates.
The idea that turbulent transports can be inferred from observed characteristics of the broad-scale hydrography goes back to Iselin (1936Iselin ( , 1939 and likely even further. It is for example the basis of the Munk (1966) estimate of diapycnal diffusivity from temperature profile curvatures below 1000 m (following upon Wyrtki, 1962). Many subsequent studies have pursued comparable inferences of turbulent transport rates (or parameters) from observed broadscale hydrography characteristics using conceptual models as well as general circulation models and adjoint techniques (e.g., Schott and Zantopp, 1980;Walin, 1982;McDougall, 1984;Olbers et al., 1985;Tziperman, 1986;Tomczak and Large, 1989;Ganachaud, 2003;Stammer, 2005;Ferreira et al., 2005;Lumpkin and Speer, 2007;Zika et al., 2009; Published by Copernicus Publications on behalf of the European Geosciences Union. Liu et al., 2012). The large influence of vertical and lateral turbulence in setting the large-scale characteristics of the ocean state in numerical models (e.g., Danabasoglu and Marshall, 2007;Eden et al., 2009;Griffies et al., 2010;Jochum et al., 2013;Gnanadesikan et al., 2014;Melet et al., 2014) underscores that the large-scale ocean state observation carries a wealth of information on turbulent transport rates.
Despite this long history and strong modeling evidence, it is not clear that inverse parameter estimates are useful, robust and/or physically meaningful. Munk's estimate of diapycnal diffusivity provides an example of the possible downfall of the method (Munk, 1966; see also Munk and Wunsch, 1998). Munk's estimate of about 10 −4 m s −1 that was originally based upon a one-dimensional model of the Pacific led to the search for the "missing mixing" since direct observations of diapycnal mixing in the thermocline were an order of magnitude smaller than Munk's value (Ledwell et al., 1998(Ledwell et al., , 2000. However, it is now commonly accepted that neglected physics (e.g., the adiabatic upwelling of North Atlantic Deep Water in the Southern Ocean; Toggweiler and Samuels, 1998;Webb and Suginohara, 2001) could explain the discrepancy.
Importantly, this example illustrates the difficulty in interpreting inverse parameter estimates and the possible confusion brought on by assumptions built in the fitted model. In light of the well-documented heterogeneity in ocean mixing rates and processes (reviewed in MacKinnon et al., 2013;Fox-Kemper et al., 2013), the notion that ocean mixing as a whole can be cast into a one-dimensional model now seems incongruous. Models of extreme simplification, while they can illuminate individual mechanisms, also discount the composite and complex essence of ocean observations. General circulation models used within a least-squares framework a priori provide a suitable framework to avoid misinterpreting observations (e.g., aliased small-scale signals) while taking advantage of complementary data sets (e.g., altimetry) and constraints (e.g., atmospheric reanalyses) to infer large-scale ocean balances and diagnose ocean variability (see Wunsch and Heimbach, 2013, for a review). However, the few publications that followed this approach to infer turbulent transport parameters (Stammer, 2005;Ferreira et al., 2005;Liu et al., 2012) provide little, if any, evaluation of the method and its results. At best, they point to the usefulness of the estimated parameters in an ad hoc manner -for example Ferreira et al. (2005) derived a parameterization of eddy diffusivity as a function of stratification that improved the unconstrained ocean model solution.
How many and which degrees of freedom may be well constrained by available observations is just one of the rather technical, yet crucial questions that have never been tackled. Thus, one could still argue that parameter inversion (along with forcing fields adjustment) through general circulation models is nothing more than an objective and practical way of tuning these models (i.e., a marginally useful approach). Whether the estimated parameters have any intrinsic value beyond the chosen estimation framework and settings (assimilated observations, period of assimilation, numerical model, etc.) remains unclear just as in the case of conceptual models. Further assessment of data constraints and parameter estimates is clearly needed and equally justified given the obvious importance of the subject matter.
In this paper we aim to demonstrate that extended observation of the broad-scale hydrography by Argo (see, e.g., Forget and Wunsch, 2007;Forget, 2010;Speer and Forget, 2013) should translate into improved estimates of upperocean turbulent transport rates. In other words, we seek to assert that turbulent transport rates are "observable" by means of Argo's collection of temperature and salinity profiles. This is clearly not a trivial proposition given the variety and heterogeneity of oceanic processes intermingled within observational data -especially in the upper ocean. Overarching oceanographic questions regarding the observability of turbulent transport rates by means of broad-scale measurements largely remain to be answered; for example, which specific features of the ocean state are informative of which turbulent transport rates? How precise may inverse estimates of turbulent transport rates be depending on the limited data availability? How are inverse problems best formulated to take full advantage of available observations?
The presented analysis reaches preliminary answers to these open questions, while providing clear supporting evidence that observations of the large-scale ocean stratification by Argo yield useful constraints to estimate turbulent transport parameters, and establishing a frame of reference for further research into their observability. It focuses on the  ocean state estimate over 1992-2011 covering the Argo era and mainly on the 0-2000 m oceanic layer that Argo observes extensively. In this framework, vertical and lateral turbulent transport parameters are estimated by fitting the simulated large-scale ocean state to observations (notably Argo profiles of temperature and salinity). Even though minimal constraints are imposed on the parameters themselves, their geography and impacts are found to be physically meaningful and shown to be useful beyond the estimation procedure.
Section 2 summarizes the estimation method, establishes that the estimated parameters are broadly consistent with the observed large-scale ocean state, and shows that they reduce spurious model drifts using independent biogeochemistry data. Section 3 demonstrates the high sensitivity of the observed upper-ocean stratification to the estimated parameters. This result provides clear supporting evidence that Argo profile collections yield a useful observational constraint of regional turbulent transport rates. The estimated turbulent parameters themselves and their relationship to the large-scale ocean state are discussed in Sect. 4. Section 5 provides a discussion of the uncertainties in the approach. The findings are summarized and perspectives are drawn in Sect. 6. Further details on the solutions and these misfits can be found in .

Reduced model errors
The Estimating the Circulation and Climate of the Ocean version 4 (ECCO v4) estimate of the evolving ocean state over the period 1992-2011 as well as the associated model and estimation settings are presented in detail in . In summary, the 20-year solution of the global 1 • resolution model is fitted to a suite of data constraints, including the vast collection of Argo profiles of temperature (T ) and salinity (S), through iterative adjustments of turbulent transport parameters (time-invariant), atmospheric forcing fields (biweekly) and initial conditions (on 1 January 1992).  show that the turbulent transport parameter adjustments are particularly important to the close fit of ECCO v4 to observed in situ profiles. They allow for a clear reduction in widespread misfits in ECCO v4 (Fig. 1, top left panel) as compared with (1) earlier ECCO solutions that optimized surface forcing fields but not turbulent transport parameters (top right corner panels) and (2) ECCO2 eddying model simulations (bottom three panels). The turbulent transport parameters being estimated, and the focus of this paper, are time-invariant threedimensional maps of bolus velocity coefficient K gm (Gent and Mcwilliams, 1990), isopycnal diffusivity K σ (Redi, 1982), and background diapycnal diffusivity K d (aside from mixed layer parameterizations). The tendency equation for a tracer φ in the ocean interior can thus be written in simplified form as ∂φ ∂t where K z includes K d plus contributions from mixed layer parameterizations (Gaspar et al., 1990;Duffy et al., 1999, plus simple convective adjustment), ∇ σ φ is the lateral tracer gradient on isopycnal surfaces (Redi, 1982), u is the Eulerian velocity, and u * is the parameterized bolus velocity representing the advective (adiabatic) effect of meso-scale eddies. After Gent and Mcwilliams (1990) the non-divergent bolus velocity is evaluated as u * = −∇ × * , where * = (K gm S y , K gm S x , 0) is the bolus streamfunction, and S x and S y are the isopycnal slopes in the zonal and meridional directions, respectively.
The first guess values are the constants K gm = 10 3 , K σ = 10 3 and K d = 10 −5 m 2 s −1 with respective uncertainties set to u gm = 500, u σ = 500 and u d = 10 −4 m 2 s −1 . As part of ECCO v4, the specification of error covariances for K gm , K σ and K d is limited to imposing smoothness at the scale of three grid points, thus allowing regional adjustments to emerge simply from observational constraints under the dynamical model constraint. The respective ranges of permitted adjustment are 10 2 <K gm <10 4 , 10 2 <K σ <10 4 , and 10 −6 <K d <5 × 10 −4 m 2 s −1 . The specified ranges span values found (directly or indirectly) in observations as well as values typically used in general circulation models. Note that the constraints imposed specifically on the turbulent parameters (smoothness and range) are minimal on purpose. One could envision imposing further constraints, for example on the vertical structure or the energetics of the turbulent parameters. Here, instead, the turbulent parameters can adjust freely within the specified ranges that reflect their large a priori uncertainty at the scale of a few grid points. Accordingly, the adjusted parameters spread out over the specified ranges as a result of observational constraints being imposed (Fig. 2). The primary objective of this paper is to reveal the observational basis to these adjustments being provided by the rapidly growing in situ data base (Sects. 3,4).
To test the estimated turbulent parameters beyond the 20 years of the estimation window, two model integrations are carried out for 500 years that perpetually loop over the 1992-2011 forcing of the ECCO v4 state estimate. One integration uses the estimated parameters while the other uses the constant first guess parameters (which, we should underscore, are common values used in ocean models of comparable horizontal resolution). The turbulent transport parameter adjustments yield a clear reduction in the model's tendency to drift away from observations over multiple centuries (Fig. 3). This is true not only for the physical variables that were directly constrained by Argo profiles (left panels) but also for biogeochemical variables that were not (middle and right panels).
The biogeochemistry result is most remarkable as it provides independent evidence that the turbulent transport parameter adjustments as estimated from observations of physical variables reduce internal model error more generally. In itself this novel result demonstrates strong constraints imposed (directly or indirectly) by observations of the physical state on biogeochemistry. The resulting improvement in biogeochemistry ( Fig. 3) implies that regional turbulent transport rates are at least partly observable by available observations of physical variables. It further motivates the assessment of the estimated turbulent transport parameters and of their observational basis presented below.
Before returning the focus to the 20-year estimation period, it is worth illuminating the consequences of K gm , K σ and K d over longer timescales (Fig. 4). The fact that the benefits of the ocean state estimation procedure extend much beyond its 20 years (Fig. 3) and through the abyss over centennial timescales (not shown) is indeed clearly not a trivial result. It implies dramatic changes in the formation and ventilation of oceanic water masses, as demonstrated by zonal mean oxygen concentrations after 500 years (Fig. 4).
The fact that the ECCO v4 turbulent transport parameter adjustments ( and independent climatological data for biogeochemistry (middle and right panels) in simulations with (blue) the K gm , K σ and K d estimated parameter maps or (red) the K gm , K σ and K d constant first guess parameters. Each plotted value is a cost function d 2 /σ 2 where d is a model-data difference, σ is an uncertainty estimate, and the overbar denotes averaging over all data points. Values of one would indicate model-data differences that are on average exactly at the estimated level of uncertainty. Argo cost function details are reported in . The biogeochemistry model is from Dutkiewicz et al. (2005) with settings provided by H. Song (personal communication, 2015). The corresponding cost functions compare annual mean biogeochemistry model fields at 300 m with the climatologies that were used to initialize the model. The 300 m climatology standard deviation is further used as an ad hoc uncertainty estimate to form a cost function. Top middle: alkalinity (Key et al., 2004); bottom middle: phosphate (Garcia et al., 2010); top right: dissolved inorganic carbon (Key et al., 2004); bottom right: dissolved oxygen (Garcia et al., 2010).  Fig. 3 caption for details) using the K gm , K σ , and K d estimate (middle) or the K gm , K σ , and K d first guess (top).
tic Bottom Water (also evident near Antarctica at all depths). The improved maintenance of a high oxygen content in the Arctic and of the mid-latitude oxygen minimum region is equally remarkable (Fig. 4). While the oxygen minimum exists due to oxygen consumption by biogeochemistry, it is also largely shaped by physical processes (Wyrtki, 1962). Figure 4 provides evidence that it is effectively constrained by Argo's collection of T and S profiles via their constraint of turbulent transport parameters revealed in Sects. 3 and 4. Whether the constraint of the oxygen minimum results directly from K gm − K gm , K σ − K σ and K d − K d appearing in the biogeochemical tracer equations or rather indirectly from improvements in other physical variables (T , S, etc.) impacting biochemistry is an interesting question left for further investigation.

Parametric controls of ocean stratification
The new constraints on turbulent transport parameter inversions brought by data collected over the past decade can most immediately be gauged using linear adjoint sensitivities ( Fig. 5 for K gm ). The results show that the constraint provided by one decade of Argo data collection (bottom) generally exceeds that provided by two decades of altimetry (middle). Furthermore, the large sensitivity increase seen upon moving from the top panel (10 years of altimetry) to the middle panel (20 years of altimetry) illustrates how the second decade of Argo data collection should even further solidify the observational basis to turbulent transport parameter inversions. These considerations motivate the detailed investigation of the existing Argo data constraint presented below.
The geography and values of ocean stratification ( ∂ρ ∂z and mixed layer depths) are a priori particularly relevant to turbulent transport rate inferences. Stratification is indeed a prime candidate amongst observational constraints as it is intimately related to potential vor- is computed with the adjoint model. The adjoint model and J Argo are documented in details in . For altimetry, J uses the large-scale formulation of . In all cases, all turbulent transport parameters are reset to their unadjusted values, so that the adjoint computation is representative of the starting point of the estimation process. The results are displayed non-dimensionally as log 10 (||u gm ∂J ∂K gm || 2 ) where || . || 2 denotes the zonal mean squared norm, and u gm = 500 m 2 s −1 . ticity, water mass formation and ventilation, which in turn strongly constrain the general ocean circulation (see e.g., Walin, 1982;Luyten et al., 1983). Argo's collection of T and S profiles yields an extensive observation of ocean stratification (Fig. 6). Its importance as an observational constraint of K gm , K σ and K d is established below.
At 300 m outside of the tropics, lows in ∂ρ ∂z (Fig. 6, bottom left) delineate regions of deep winter convection (Fig. 6, bottom right) and mode water formation (see Speer and Forget, 2013). In the tropics at 300 m, lows in ∂ρ ∂z denote thermoclines shallower than 300 m. Conversely, subtropical and tropical thermoclines intersecting 300 m are marked by highs in ∂ρ ∂z (Fig. 6, bottom left). The state estimate's stratification geography and values (Fig. 6, middle left) generally are in very good agreement with Argo observations (Fig. 6, bottom  left). Thus, the estimated K gm , K σ and K d are consistent with the observed stratification.
To evaluate the full non-linear sensitivity of ocean stratification to each of the three parameters, three forward perturbation experiments are conducted. In each experiment, one parameter is individually reset to its constant first guess value. It is worth recalling that atmospheric forcing field adjustments (in temperature, specific humidity, downward radiation, precipitation, and wind stress) have a lesser impact on the subsurface hydrography than turbulent transport parameter adjustments do in ECCO v4 . In particular it is readily evident in Fig. 6 that the realism of the state estimate stratification is noticeably improved by the turbulent transport parameter adjustments (compare middle and top panels). Therefore, the adjusted atmospheric forcing fields are retained in all perturbation experiments presented here to focus on turbulent transport parameters and their observability.
The mean squared deviation of ∂ρ ∂z between each perturbed solution and the state estimate over the period 2008-2010, normalized by the corresponding ∂ρ ∂z variance, is shown in logarithmic scale as a function of latitude and depth (Fig. 7). Logarithmic values above 0 indicate ∂ρ ∂z contrasts between solutions being as large as the zonal and seasonal mean ∂ρ ∂z variance in the state estimate. Such large impact of turbulent transport parameter adjustments on ∂ρ ∂z occurs at all latitudes in the upper 2000 m. Figure 7 thus provides decisive supporting evidence that inversion of regional turbulent transport rates is effectively guided by the constraint of fitting Argo profiles. The collocation of large ocean stratification changes (Fig. 7) and turbulent transport parameter adjustments (Fig. 2) further confirms that they are tied to each other in ECCO v4 (see Sect. 4 for more detail).
Generally, in the upper 2000 m and over the whole water column at high latitudes, ocean stratification is found to be more sensitive to K gm (Fig. 7, top) than to K σ (middle) or K d (bottom). The predominant impact of K gm suggests that ocean stratification profiles may most efficiently constrain (i.e., observe) the rates of bolus advection, albeit with several noteworthy exceptions. High sensitivity to K σ (Fig. 7, middle) is found at high latitudes. The impact of K σ on ∂ρ ∂z is indirect since ∇ σ φ in Eq. (1) vanishes for φ = σ . It does however have dynamical impacts and it is not surprising that they are magnified at high latitudes where dense layers outcrop and interact with mixed layer and sea-ice processes (e.g., see England, 1993;Zika et al., 2009). Furthermore, haline density variations overcome thermal density variations at those latitudes (e.g., see Forget and Wunsch, 2007; and this transition provides an environment conducive to dynamical impacts of K σ . High sensitivity to K d is found in the tropics and in the Arctic in the 0-2000 m layer, as well as near the sea surface at all latitudes. This result implies that observation of the stratification's broad structure informs us of diapycnal diffusion rates in these regions. Winter mixed layers modulate water mass transformations and the penetration of surface buoyancy fluxes. As such, they are a key factor in ocean stratification and deserve special attention. The state estimate closely agrees with the geography Ocean Sci., 11, 839-853, 2015 www.ocean-sci.net/11/839/2015/ Stratification at 300 m depth 90th percentile mixed layer depth Figure 6. (Left) median stratification at 300 m depth (shown as log 10 ( ∂ρ ∂z ) in kg m −4 ) and (right) 90th percentile mixed layer depth (shown as log 10 (mld) in m) for (bottom) in situ profiles, (middle) the corresponding state estimate profiles, and (top) model profiles generated by resetting K gm , K σ , and K d to the K gm , K σ , and K d first guess while retaining all other settings of the state estimate. Percentiles are computed from the distribution of individual profile values within grid boxes and mapped. and magnitude of observed mixed layers (Fig. 6, right panels). The sensitivity of mixed layer depths to turbulent transport parameters (Fig. 8) varies on a regional basis and mostly reflects the underlying stratification sensitivity (Fig. 7). Logarithmic values above 0 in Fig. 8 indicate that changes in mixed layer depths that result from the parameter adjustments are as large as the seasonal contrasts seen in the state estimate. Such differences in mixed layer depths provide additional evidence of efficient constraints imposed by Argo on turbulent transport rates.
Several dynamical regimes can be distinguished in Figs. 7 and 8 that are in broad agreement with theoretical views and worth commenting upon before proceeding with further assessment of the geography (Sect. 4) and uncertainty (Sect. 5) of K gm , K σ and K d . The thermocline in the upper 2000 m (between 60 • S and 60 • N) is primarily adjusted through adiabatic circulation controls (Fig. 7, top panel) rather than di-apycnal mixing changes (Fig. 7, bottom panel). This behavior is prompted by data constraints leading to K gm , K σ and K d rather than by a priori assumptions. Thus, Fig. 7 provides objective observational support for predominantly adiabatic, ventilated thermocline theories (see Vallis, 2006).
The sizable impact of the increased K d at the base of the ventilated thermocline further supports the notion that it is embedded in an internal thermocline (Samelson and Vallis, 1997) where diffusion balances advection (see Figs. 2 and 7). Another regime appears near the sea surface in the tropics, where shallow mixed layers and stratification also strongly respond to background diapycnal mixing (Figs. 7, 8). Finally, high latitudes where deep isopycnals outcrop and experience large isopycnal variations in T and S expectedly show a sizable impact of K σ , along with the predominant impact of K gm (Figs. 7, 8) Figure 7. Sensitivity of ocean stratification to K gm (top), K σ (middle), and K d (bottom). Mean squared deviations are computed between perturbation experiments and the state estimate for the 2008-2010 monthly climatology of ∂ρ ∂z , and then normalized by the corresponding state estimate variance. The color shading shows the log 10 of this ratio. Thus, orange denotes >10 % differences, while dark red denotes >100 % differences. The perturbation experiment where the K gm (resp. K σ , K d ) estimate is reset to the constant K σ (resp. K σ , K d ) first guess is shown at the top (resp. middle, bottom). Overlaid contours: zonal mean potential density from the OCCA atlas.
cal regimes points to the importance of not overly simplifying the inversion problem and the underlying model.

Estimated turbulent transport parameters
The results from Sects. 2 and 3 suggest that the estimated K gm , K σ and K d may have oceanographic value beyond the ECCO v4 state estimate itself -a notion supported by Fig. 3 in particular. The main goal of this section is to further assert the relationship between the observed ocean stratification and estimated turbulent transport parameters, and to assess whether the estimated parameter maps are physically meaningful.
Each of the estimated parameters varies by orders of magnitude on a regional basis and thus shows a great degree of heterogeneity (Fig. 9) -a physically reasonable proposition (see, e.g., Eden, 2006;Eden et al., 2007;Fox-Kemper et al., 2013). Highs and lows often alternate at the same latitudes over extended regions. K gm , K σ and K d however also show Figure 8. Mixed layer depth sensitivity (color shading) to K gm (top), K σ (middle), and K d (bottom). Computational and plotting details are similar to Fig. 7. Overlaid blue contours (resp. magenta contours) denote the 60, 70, 80, 90th (resp. 10, 20, 30, 40th) percentiles of the observed mixed layer depth map (mld(x, y) from the bottom right panel of Fig. 6). extended regions of limited departure from K gm , K σ and K d (regions in green). This behavior ought to be expected in regions where turbulent transports may be weak regardless of the corresponding coefficients.
Patterns of large parameter adjustments generally align with contours of observed stratification (Fig. 9). Stratification contrasts are indeed expected to delineate between regimes of turbulence, where one or another physical process may become predominant. Parameterized turbulent transports furthermore act along or across isopycnals, and stratification contrasts are indicators of the isopycnals' geography.
Parameterized advection of tracers by meso-scale eddies (see Eq. 1) is controlled by K gm . Reduction in K gm is most distinctly seen at the equatorward flank of the subtropical thermocline bowl in all oceanic basins (Fig. 9, top). Low K gm values appear consistent with the observed tropical thermocline characteristics (Figs. 1, 6). The tropics are indeed known to show relatively low levels of meso-scale eddy activity and a predominance of baroclinic wave trains instead (e.g., Tulloch et al., 2009) Each overlaid contour corresponds to a percentile of the observed stratification map ∂ρ ∂z (x, y) that is depicted at 300 m in Fig. 6 (bottom left). Black contours (resp. magenta contours) denote the 60th, 70th, 80th, and 90th (resp. 10th, 20th, 30th, and 40th) percentiles of ∂ρ ∂z (x, y).
K gm acts to increase the equatorward slant of isopycnals per Eq. (1). Instances of K gm ≥ K gm are common at higher latitudes where large isopycnal slopes are found. K gm >K gm is found at subpolar latitudes and not necessarily at the cores of the Antarctic Circumpolar Current (ACC) and subtropical jets. We refer the reader to Ferreira et al. (2005) and Liu et al. (2012) for discussions of comparable results. Following the various ACC branches and meanders, instances of both K gm <K gm and K gm >K gm are found depending on longitude, as should be expected from Eden (2006) and subsequent studies. Furthermore, vertical contrasts in K gm are readily apparent in Fig. 9, such as the appearance of K gm >K gm patches in the ACC between 300 and 900 m depth. They become even more predominant at greater depths (see Sect. 5).
Isopycnal tracer diffusion due to lateral turbulence is controlled by K σ (see Eq. 1). Similar adjustment patterns can be seen in K σ (Fig. 9, middle) and K gm (Fig. 9, top). Notably, K σ shows minima in the tropics at 300m and below (in qualitative agreement with Cole et al., 2015) that correspond to K gm minima (see also Fig. 2). This is not entirely surprising since the two parameterized processes, while generally distinct, are associated with lateral turbulence, so that K σ and K gm may covary (see, e.g., . It should be noted however that the estimation settings did not impose covariance between the two parameters, but it is allowed to emerge from data constraints if adequate. Differences between K σ and K gm could conversely indicate rewww.ocean-sci.net/11/839/2015/ Ocean Sci., 11, 839-853, 2015 gions where the ocean state is not equally sensitive to the diffusive and advective effects of lateral turbulence. One notable feature in K σ (Fig. 9, middle right) that is not in K gm (Fig. 9, top right) is its reduction from K σ near 900 m on the southern side of the ACC and near Antarctica (see also Fig. 2). It is physically reasonable that isopycnal diffusion is a particularly sensitive process in regions showing large tracer gradients along isopycnals (see Gnanadesikan et al., 2014). Also notable are two contrasting situations in the North Atlantic. The subpolar gyre shows increased K gm and K σ . The eastern Atlantic subtropics at 900 m (near the Mediterranean outflow depth) show reduced K gm but increased K σ . As emphasized by , K gm and K σ are not expected to be equal except under specific limits. Gnanadesikan et al. (2014) further illustrate that tying K σ strictly to K gm (as is often done in ocean modeling) can be problematic. The contrasted cases seen in Figs. 2 and 9 provide observational evidence that K σ and K gm are not systematically tied to each other.
Background diapycnal diffusivity K d is increased from K d at 300 m near 30 • latitude (Fig. 9, bottom left). This geography qualitatively agrees with the fine structure observed by Argo ( Fig. 1 in Whalen et al., 2012) and theoretical predictions of intensified parametric subharmonic instability (e.g., MacKinnon and Winters, 2005). However, interleaving of weak and strong mixing layers is a common feature of the K d inverse estimate. Hence, K d at 900 m (Fig. 9, bottom right) is rather reduced near 30 • latitude.
A pronounced interleaving is also seen in K d near the Equator (Fig. 10). Increased K d in the upper 100 m is consistent with the analysis of (amongst others) Moum et al. (2009) and found in all basins. A secondary K d maximum can be seen in the Pacific and Atlantic immediately underneath the Equatorial Undercurrent. Shear instability is a good candidate mechanism also in this case. The deepest Pacific maximum (near 1000 m in Fig. 10) is a tropical rather than equatorial feature (Fig. 9, bottom right). Its dynamical origin is likely very different from the two upper maxima, and associated with the internal thermocline (see Sect. 3; Fig. 7).
In summary, this section shows that the regional contrasts in K gm , K σ and K d are tied to regional stratification contrasts and that these estimates generally are physically plausible. It is not implied however that turbulent transport rates are uniquely dependent on stratification alone and everywhere. The interplay of instability processes that in turn govern turbulent transports and stratification on a regional basis certainly remains a subject of active research of great importance beyond this paper. The presented analysis simply supports the notion that turbulent transport rates are effectively constrained (i.e., observable) by means of Argo's extensive collection of stratification profiles.

Assessment of uncertainties
A comparison of ECCO v4 parameter estimates (Figs. 9, 11) with earlier inversion results that did not cover the Argo era (Ferreira et al., 2005;Stammer, 2005;Liu et al., 2012) should be indicative of overall observational uncertainty levels. None of the listed estimates is provided with a formal error estimate, which more generally remains a major caveat of ocean modeling and data synthesis. In such context, intercomparison of solutions is a commonly accepted, practical method to assess uncertainties even though its results are often difficult to interpret precisely (see Danabasoglu et al., 2014;Balmaseda et al., 2014).
In particular, attributing point by point differences amongst inverse estimates to specific causes would be a perilous exercise, due to various differences in model and estimation settings and therefore no such attempt is made. It is worth highlighting, however, that ECCO v4 benefits from many innovations (e.g., updated numerics, the addition of a sea ice model and of the Arctic, and increased vertical resolution) as compared with previous generation model setups used by Ferreira et al. (2005), Stammer (2005, and Liu et al. (2012). A more exhaustive list of innovations is provided by .
The most meaningful comparison may be between the ECCO v4 and Liu et al. (2012) results since their respective experimental settings are most comparable. Importantly, estimated parameter adjustments are larger in ECCO v4 than in Liu et al. (2012) (compare Fig. 11 with Figs. 7c, 12b and13b in Liu et al., 2012). It should be noted that the impact of the K gm , K σ and K d adjustments in ECCO v4 was shown to generally exceed the impact of model errors unrelated to turbulent transport parameterizations . Therefore, the increased parameter adjustment amplitude (in ECCO v4 as compared with Liu et al., 2012) is thought to primarily reflect the extensive data constraints added over 2002-2011) rather than differences in model settings.
Regardless of the noted caveats, it is encouraging that ECCO v4 parameter estimates (Figs. 9, 11) bear some resemblance to the Ferreira et al. (2005), Stammer (2005 and Liu et al. (2012) results, which may indicate robust oceanic features. All estimates are rich in regional adjustment patterns aligned with contours of the large-scale hydrography. The three K d estimates (Stammer, 2005;Liu et al., 2012 and ECCO v4) show elevated mixing near 30 • latitude at 300 m, and interleaving of high and low mixing in the tropics. The K d map of Liu et al. (2012) at 300 m (their Fig. 7a) more generally is in a good qualitative agreement with Fig. 9 (bottom left). As a final example, the three K gm estimates (Ferreira et al., 2005;Liu et al., 2012 and ECCO v4) show maxima associated with the ACC, minima in the tropics, and maxima in the North Atlantic subpolar gyre.
There are, however, also many differences amongst inverse parameter estimates. Firstly, ECCO v4 shows intensified diapycnal mixing in the tropical thermocline below 500 m (Fig. 11, bottom panel) as compared with Stammer (2005) and Liu et al. (2012). Secondly, the three K gm estimates largely differ in their magnitude and vertical structure in the ACC and at 40 • N. Ferreira et al. (2005) show increased K gm in the upper 1000 m in the ACC, with a nearsurface maximum, whereas Liu et al. (2012) and ECCO v4 (Fig. 11, top panel) show increased K gm between 1000 and 3000 m. In the Gulf Stream, Ferreira et al. (2005) show a K gm maximum in the upper 1000 m, whereas Liu et al. (2012) and ECCO v4 show K gm minima at these depths. Also, neither Ferreira et al. (2005) nor Liu et al. (2012) show the increased K gm seen near the coast of Antarctica in ECCO v4. Finally, K σ in ECCO v4 (Fig. 11, middle panel) hints at a steering level effect (Green, 1970;Abernathey et al., 2010;Ferrari and Nikurashin, 2010) below 500 m in the ACC that is not seen in Liu et al. (2012). Such differences are indicative of large overall uncertainty in inverse parameter estimates.
Furthermore, regions where parameters remain virtually unadjusted (Figs. 2,9;regions in green) likely denote large uncertainties reflecting that available data constraints are insufficiently sensitive to prompt sizable parameter adjustments (as noted by Liu et al., 2012). The relatively weak values of K d − K d in the ACC (as compared with, e.g., Whalen et al., 2012, but not with Liu et al., 2012 may be one example. Similarly, the fact that K d − K d is generally muted in the abyss (albeit with notable exceptions in the Southern Ocean) is not surprising and does not imply that K d is a precise first guess. Indeed the equilibration of the abyssal stratification (or the lack thereof) and the recycling of abyssal water masses are dominated by very long timescales, and abyssal observations are very sparse as compared with upper-ocean data constraints.
It is in fact encouraging that K d shows even marginal increases near the sea floor (Fig. 12, right; Fig. 2, bottom right) as it is often expected to result from the interaction of barotropic tides (amongst others) and bottom topography (Polzin et al., 1997;Ledwell et al., 2000;Naveira Garabato et al., 2004;Sloyan, 2005). It is intriguing that low K d values are also found near the bottom -mostly in the Southern Ocean along deep canyon margins (Fig. 12, right; Fig. 2, bottom left). However, reductions in diapycnal diffusivity (from K d to K d ) may compensate for an unknown amount of numerical diffusion implied by the advection schemes, and bottom boundary layers are notoriously difficult to simulate adequately in ocean models. For these reasons, and since they are not seen in Stammer (2005) or Liu et al. (2012), the K d contrasts in Fig. 12 should be interpreted most cautiously.
The deep Southern Ocean can adjust relatively fast due to the proximity of bottom water formation sites and the presence of a deep wind-and eddy-driven thermocline (see Karsten and Marshall, 2002) that result in a strong coupling of superficial and deep layers in the ACC region (see Fig. 7). Hence the large values of K gm at 3000 m estimated over 20 years in the Southern Ocean may be physically meaningful (Fig. 12, left). Maxima in K gm are located along the ACC path just downstream of Kerguelen and at Drake Passage, as well as in the Brazil-Malvinas confluence region, in the Ross Sea and in the Weddell Sea. These regions are indeed characterized by relatively large isopycnal slopes, and eddying numerical models show sizable meso-scale eddy activity even at 3000 m along the ACC (Ponte, 2012 at 3000 m depth at the sea floor to see whether the K gm maxima seen in Fig. 12 are confirmed (or otherwise) in inversion experiments conducted once deep Argo profiles are available.

Summary and perspectives
This study asserts that extended observation of the broadscale hydrography by Argo should translate into improved inverse estimates of regional turbulent transport rates in the upper 2000 m over the global ocean. Time-invariant threedimensional maps of turbulent transport parameters (K gm , K σ and K d ) are estimated as part of ECCO v4  and under the observational constraint of Argo T -S profiles collected through 2011. The presented exploration of the method of turbulent transport parameter inversion, while still incomplete, fills a major gap in the oceanographic literature. The observability of turbulent transport rates is asserted by focusing on ocean stratification ( ∂ρ ∂z and mixed layer depths). Argo indeed readily observes these key oceanographic variables with unprecedented data coverage. The estimated turbulent transport parameters are consistent with the observed stratification by construction, since it is well reproduced by the ECCO v4 ocean state estimate. It is shown that ocean stratification over the upper 2000 m is highly sensitive to the estimated parameter adjustments, and that the geography of estimated parameter adjustments is aligned with contours of observed ocean stratification. Thus, the constraint of fitting Argo profiles is identified as an effective observational basis for the inversion of turbulent transport parameters.
As part of the inversion method evaluation it is shown that the estimated parameters reduce spurious model drifts (i.e., accumulating model biases) of the physical ocean state in multi-century simulations. They also lead to remarkable improvements in simulating biogeochemistry variables that were not involved in the parameter optimization. The estimated parameter adjustments themselves and the resulting adjustments in ocean stratification are physically plausible despite the minimal constraints that were built in the optimization. These results demonstrate that the estimated parameters have intrinsic value beyond the optimized solution of the 20-year evolving ocean physical state.
The asserted "observability" of turbulent transport rates by Argo does not, however, require or imply that the present parameter estimates are very precise or accurate. Given the noted contrasts amongst published inversion results, and given that vast regions show negligible parameter adjustments, it is unlikely that available T -S profiles suffice to determine regional turbulent transport rates uniquely and everywhere. The overall weakness of estimated parameter adjustments in the abyss is also revealing of current limitationsextensive Argo data collection has only reached 2000 m and 20 years is too short a period to fully resolve (im)balances of the abyss. Additional observational constraints (e.g., passive and biogeochemical tracer observations), statistical constraints (e.g., observed fine-scale and meso-scale statistics) and dynamical constraints (e.g., longer timescales, energetics) ought to complement the constraint of T -S profiles in future inversion experiments.
The lack of a practical technology to associate estimated turbulent transport parameters with formal error estimates is arguably the main caveat here as well as in Stammer (2005), Ferreira et al. (2005), and Liu et al. (2012). It is the reason why observability of turbulent transport rates by means of Argo and other global data sets largely remains to be quantified. The presented inter-comparison of inversion results, however, provides clues into overall levels of uncertainty. The closest agreement appears to be between K d estimates at 300 m. A much more contrasted picture emerges for e.g., K gm in the main mid-latitude jets (ACC, Gulf Stream, Kuroshio). The three inverse estimates of K gm (Ferreira et al., 2005;Liu et al., 2012; show both minima and maxima along or near the jets but with differences in their geographic and vertical distributions. Aside from the need for formal error estimates, and additional constraints, much remains to be done to refine inverse estimates of turbulent transport rates and our understanding Ocean Sci., 11, 839-853, 2015 www.ocean-sci.net/11/839/2015/ of their observability. One practical approach would consist in conducting twin estimation experiments (see, e.g., Forget et al., 2008a) focusing on turbulent transport parameter inversions. Another one would consist in conducting dedicated estimation experiments where real data sets are withheld or added one at a time to further compare the constraints that they respectively provide (see, e.g., Forget et al., 2008b). And the presented results should eventually be re-evaluated on the basis of additional estimation experiments that would differ from  in regard of parameter ranges, first guess values, and error covariance specifications (all of which are uncertain estimation settings).
Ongoing research aiming to diagnose and alleviate numerical diffusion (e.g., Hill et al., 2012) and structural model uncertainty  is of direct relevance to turbulent transport parameter inversions. Numerical diffusivity of advection schemes or e.g., momentum equation biases could very well contaminate turbulent transport parameter estimates.  show that turbulent transport parameter adjustments exceed what may be expected to compensate for model errors due to advection and momentum scheme choices. However, turbulent transport parameter inversions will need to be conducted with a variety of numerical models before one can reach more definitive conclusions in this regard. Initially, it will be interesting to see which beneficial impacts of the presented K gm , K σ and K d carry over to different ocean and climate models.