Tracer and observationally-derived constraints on diapycnal diffusivities in an ocean state estimate

Use of an ocean parameter and state estimation framework–such as the Estimating the Circulation & Climate of the Ocean (ECCO) framework–could provide an opportunity to learn about the spatial distribution of the diapycnal diffusivity parameter (κρ) that observations alone cannot due to gaps in coverage. However, we show that the inclusion of misfits to observed physical variables–such as in situ temperature, salinity, and pressure–currently accounted for in ECCO is not sufficient, as κρ from ECCO does not agree closely with any observationally-derived product. These observationally-derived κρ 5 products were inferred from microstructure measurements, derived from Argo and CTD data using a strain-based parameterization of finescale hydrographic structure, or calculated from climatological and seafloor data using a parameterization of tidal mixing. The κρ products are in close agreement with one another, but have both measurement and structural uncertainties, whereas tracers can have relatively small measurement uncertainties. With the ultimate goal being to jointly improve the ECCO state estimate and representation of κρ in ECCO, we investigate whether adjustments in κρ due to inclusion of misfits to 10 a tracer–dissolved oxygen concentrations from an annual climatology–would be similar to those due to inclusion of misfits to observationally-derived κρ products. We do this by performing sensitivity analyses with ECCO. We compare multiple adjoint sensitivity calculations: one configuration that uses misfits to observationally-derived κρ and the other uses misfits to observed dissolved oxygen concentrations. We show that adjoint sensitivities of dissolved oxygen concentration misfits to the state estimate’s control space typically direct κρ to improve relative to the observationally-derived values. These results suggest that 15 the inclusion of oxygen in ECCO’s misfits will improve κρ in ECCO, particularly in (sub)tropical regions.


Introduction
We consider the challenges with using observational data in the context of a parameter and state estimation framework to infer the global distribution of ocean mixing. Ocean models must parameterize the unresolved, turbulent diffusion of oceanic tracers. further distinguish "observationally-inferred" values, which are from the currently accepted method of observing a quantity such as κ ρ but are not measured, and "observationally-derived" values because the latter data depend on a method that requires additional assumptions. These terms only apply to values calculated making use of observations.) κ ρ,micro are based on an expression for the isotropic turbulence field, which is proportional to the viscosity of water and the velocity shear resolved to 125 dissipative scales (Thorpe (2007); and references therein). The depth ranges of the data collected by Waterhouse et al. (2014) go from the upper several hundred meters to the full water column. The profiles are seasonally biased at higher latitudes and span decades. There are thousands of vertical profiles from 24 different campaigns that comprise this data set, with samples being taken in North Pacific Ocean, North Atlantic Ocean, tropical Pacific, near Drake Passage, near the Kerguelen Plateau, and in the South Atlantic Ocean. Many of the profiles were taken in regions with both smooth and rough bottom topography. To 130 compare the microstructure profiles with model output, the nearest neighbors to each model's grid are selected. A geometric average is taken for each profile because this is more representative than an arithmetic average for a small sample size and when the data are not normally distributed (Manikandan, 2011), like the log-normal distribution of κ ρ (Whalen, 2021).
We make use of multiple data sets for κ ρ derived from observations. Two of these data sets-listed in Table 1 κ ρ values are derived from finestructure observations of temperature, salinity, and pressure using a strain-based finescale parameterization, which has been developed and implemented in different ways (Henyey et al., 1986;Gregg, 1989;Polzin et al., 1995Polzin et al., , 2014, but typically assumes a mixing efficiency of 0.2 (St. Laurent and Schmitt, 1999;Gregg et al., 2018). The finescale parameterization assumes that 1) the production of turbulent energy at small scales is due to an energy transfer driven by water areas-which skew its global average (standard deviation) below the mixed layer to about 1.2 × 10 −4 (3 × 10 −3 ) s −2 .
The vertical gradients in N 2 are typically between −10 −5 and 10 −5 m −1 s −2 and have an average value (standard deviation) 190 of about −10 −7 (4 × 10 −6 ) m −1 s −2 , with its largest magnitudes on continental shelves-high-latitude ones in particular-and in the eastern equatorial Pacific Ocean. The spatial correlation between the annual mean vertical gradients in oxygen (Figs. 1a,c,e) and the annual mean vertical gradients in N 2 is about 0.25, which suggests that stratification is one candidate factor in explaining why oxygen concentrations are correlated with κ ρ and . However, we do not test this with model experiments that incorporate information about N 2 -which directly compare N 2 from our model and observations-because the vertical 195 resolution of our ocean model is so much coarser than that of observations. Instead, we run a set of model experiments that compare oxygen concentrations, κ ρ , or . We perform these model experiments to further explore the potential information that oxygen concentrations provide about κ ρ and indirectly infer-via the Osborn (1980) relation-the possible role of stratification.

Modeling system
We use the Estimating the Circulation & Climate of the Ocean (ECCO) framework in our analysis. ECCO uses a time-invariant 200 but spatially varying background κ ρ field, calculated with a parameter and state estimation procedure, where κ ρ associated with temperature and salinity are assumed to be identical. Details about the model simulations we perform are summarized in Table 2.

ECCO
The modeling system used here is ECCOv4r3 (Fukumori et al., 2017). The underlying ocean-sea ice model is based on the  (Adcroft et al., 1997). The MITgcm uses a dynamic/thermodynamic sea ice component (Menemenlis et al., 2005;Losch et al., 2010;Heimbach et al., 2010) and a nonlinear free surface with freshwater flux boundary 210 conditions (Campin et al., 2004). The wind speed and wind stress are specified as 6-hourly varying input fields over 24 years . Average adjustments to the wind stress, wind speed, specific humidity, shortwave downwelling radiation, and surface air temperature are re-estimated and then applied over 14-day periods. These adjustments are based on estimated prior uncertainties for the chosen atmospheric reanalysis (Chaudhuri et al., 2013), which is ERA-Interim (Dee et al., 2011). The net heat flux is then computed via a bulk formula (Large and Yeager, 2009). The ocean variables, on the other hand, do not get pe-215 riodically adjusted. A parameterization of the effects of geostrophic eddies (Gent and McWilliams, 1990) is used. Mixing along isopycnals is accounted for according to the framework provided by Redi (1982). Vertical mixing-diapycnal plus the vertical component of the along-isopycnal tensor-is determined according to the Gaspar et al. (1990) mixed layer turbulence closure, simple convective adjustment, and estimated background κ ρ . Here, κ ρ represents a combination of processes, including-but potentially not limited to-internal wave-induced mixing. κ ρ , the Redi coefficient, and the Gent-McWilliams coefficient are 220 time-independent because of the under-determined problem of inverting for initial conditions and model parameters would be even more under-determined if they were allowed to vary in time-explained below. In order to simulate oxygen concentrations, tracers are carried using Biogeochemistry with Light, Iron, Nutrients and Gases (BLING) model (Galbraith et al., 2015).

BLING is an intermediate complexity biogeochemistry model that uses eight prognostic tracers and parameterized, implicit
representations of iron, macronutrients, and light limitation and photoadaptation. BLING has been shown to compare well with 225 the Geophysical Fluid Dynamics Laboratory's full-complexity biogeochemical model, TOPAZ (Galbraith et al., 2015), and has been adapted for use in the MITgcm with its adjoint (Verdy and Mazloff , 2017). Initial conditions and model parameters for the runs performed here are from ECCOv4r3. The least-squares problem solved by the ECCO model uses the method of Lagrange multipliers through iterative improvement, which relies upon a quasi-Newton gradient search (Nocedal, 1980;Gilbert and Lemarechal, 1989). Algorithmic (or automatic) differentiation tools (Griewank,230 1992; Giering and Kaminski, 1998) have allowed for the practical use of Lagrange multipliers in a time-varying non-linear inverse problem such as ocean modeling, eliminating the need for discretized adjoint equations to be explicitly hand-coded.
Contributions of observations to the model-data misfit function are weighted by best-available estimated data and model representation error variance (Wunsch and Heimbach, 2007). The observational data included in the ECCO state estimation pro-sitivity of one variable to another, computed by making use of the model's adjoint. Formally, an adjoint sensitivity is ∂J/∂X, where the cost function J is a sum of weighted misfits to observations and a control variable X is a variable that the model estimates by making use of its adjoint and observations-see Section 2.2.2. The adjoint sensitivities provide information about which directions the model's input parameters should change X in order to minimize J. Masuda and Osafune (2021) showed 260 some examples of adjoint sensitivities of several model parameters in their ocean state estimate to a vertical mixing parameter (slightly different from κ ρ ). We also compute adjoint sensitivities in the present study, but using ECCO with respect to X = κ ρ .
The following is a summary of the ECCO experiments we run (Table 2): -E-CTRL -a forward ECCOv4 simulation that uses the parameters from ECCOv4r3; this simulation can be referred to as a "re-run" 265 -E O -an adjoint sensitivity (with respect to X = κ ρ ) experiment in which only oxygen concentrations from the World Ocean Atlas (2013) climatology are included in the misfit function J in our ECCO run using BLING are shown in Fig. 3b. The differences are largest in the Arctic Ocean, northeastern Pacific 290 Ocean, and near the coasts, particularly on the eastern side of the American continent, the southwestern side of the African continent, around the Kuroshio/Sea of Japan region, along almost every coastline of Oceania, and in the Mediterranean Sea ( Fig. 3b). Point-wise differences between the initial conditions for oxygen concentrations in ECCO and the World Ocean Atlas (2013) product are shown in Fig. 3c, which suggests that there is strong agreement between the two fields. Where there are disagreements, the initial conditions for oxygen concentrations in ECCO are more often too small (particularly in the Atlantic

295
Ocean, as shown in Fig. 3b) than too large. These differences are likely due to the deficiencies in model resolution, the sparse observations in regions such as the Arctic Ocean, the locations of sea ice (Bigdeli et al., 2017), and the parameterization of the tracer air-sea fluxes (e.g., Atamanchuk et al. (2020)). We need to consider the spatial patterns shown in Fig. 3b when interpreting the signs of the adjoint sensitivities.

ECCO adjoint sensitivity analyses 300
Short of including a particular data set (e.g., dissolved oxygen concentrations) in the misfits of a new optimization run of ECCO, we assess whether the inclusion of a particular data set in the model's misfits could lead to a more accurate estimate of a control variable that can be observed (e.g., κ ρ ). In order to understand whether κ ρ could be estimated more accurately through the inclusion of oxygen concentrations in the model's misfit, we need to further explain the details of our adjoint sensitivity experiments with ECCO. We define the objective (or cost) function here to more formally explain what the adjoint sensitivity 305 is. ECCO calculates the cost function to be minimized, J, (Stammer et al., 2002) as-focusing here only on the observational misfit terms while omitting regularization terms for the control variables: where t f is the final time step,x is the model-based estimate of the state vector x, S is the observation matrix that relates the model state vector to observed variables y (such that Sx is the model-based estimate of the observables y), and W is the 310 weight (inverse square of approximate uncertainties accounting for measurement and representation errors) of the observations.
In each of our adjoint sensitivity experiments, the data vector y only contains the data set specific to the experiment (see Table   2) so we emphasize here that J is different for each of our experiments. The uncertainties in κ ρ,ECCO in E κ are set to be three times the values of the observationally-derived κ ρ because of the level of agreement between the κ ρ,Argo and κ ρ,micro (Whalen et al., 2015). The uncertainties in oxygen concentrations in E O are set to be 2% of the values of the measured dissolved oxygen 315 concentrations.
The adjoint sensitivities computed in this study are the derivatives of J in Eq. 1 with respect to κ ρ . We consider evaluating directions in the control space in which to improve κ ρ , given the control vector from the ECCOv4r3 solution. While the adjoint sensitivities of J to the control space in experiment E O must be computed online, those in E κ can either be computed online or offline using an analytical equation (see below). The adjoint sensitivity run with κ ρ included in the misfit calculation of 320 experiment E κ can be calculated offline using output from the E-CTRL run instead of being calculated online as follows: (2) Here, X = κ ρ is the control variable, X obs is the observationally-derived value of X described in the previous section, X model is the value that ECCO estimates for X, and σ X is taken to be 3X obs /1.96 (or the base-10 logarithm of this in the case of κ ρ ) due to the factor of 3 uncertainty corresponding to an approximate 95% confidence interval in Whalen et al. (2015). For X model ,

325
we use the offline values calculated from the E-CTRL run following Eq. 2. While this assumes a diagonal W and minimal impact of the smoothing operator applied over a decorrelation length scale diameter of ∼ 100 km, the offline Eq. 2 and online sensitivities have been verified to be in agreement.
Because the observations of κ ρ here are not direct measurements, we first need to show that observationally-derived κ ρ has a smaller bias with respect to independent observations than the model's estimate of κ ρ . We devote the first portion of our 330 study to determining whether |κ ρ,Argo − κ ρ,micro | < |κ ρ,ECCO − κ ρ,micro | (and, by extension, κ ρ,CT D in place of κ ρ,Argo ) is true. We do this because κ ρ,micro is limited in its spatial coverage compared to κ ρ,Argo , κ ρ,CT D , and κ ρ,tides . Also, κ ρ,Argo and κ ρ,CT D are still limited spatial coverage relative to dissolved oxygen concentrations. While κ ρ,tides has global spatial coverage, its measurement plus structural uncertainties are not well-known compared to dissolved oxygen concentrations. The data product with higher accuracy (dissolved oxygen concentrations) will have larger weights (W in Eq. 1) and thus will exert 335 more influence in constraining κ ρ,ECCO -bringing it closer to microstructure values. So if we can show that the adjustments to κ ρ in ECCO are similar, whether we provide information from observationally-derived κ ρ or a measured tracer with relatively small uncertainties (dissolved oxygen concentrations), then we would include the tracer in the misfits.
One problem with doing a direct comparison of the adjustments is that the uncertainties in observationally-derived κ ρ products are large, so we first quantify the extent to which the adjoint sensitivities from two runs (here, E κ and E O ) have 340 the same sign at each location and depth. Specifically, we inspect whether ∂J/∂κ ρ has the same sign in E κ and E O where |κ ρ,Argo − κ ρ,ECCO | is significantly different from zero (i.e., κ ρ,Argo is more than a factor of three greater or less than a factor of three smaller than κ ρ,ECCO ). We are interested in regions where κ ρ is significantly erroneous and where the errors in oxygen are due to errors in the physics (e.g., κ ρ ), not initial conditions; hence, these choices. We perform these comparisons in regions where the difference between the observationally-derived κ ρ products and κ ρ,ECCO exceeds three times the observational 345 products' magnitudes (i.e., statistically distinguishable from zero). Because model errors unrelated to κ ρ can confound the correlations between the adjoint sensitivities from E κ and E O , we additionally look at regions where the difference between oxygen concentrations from the model and the World Ocean Atlas (2013) is relatively small to determine whether oxygen concentrations guide the state estimate's control space to improve the magnitude of κ ρ . In this subset of regions, we calculate the correlations between the adjoint sensitivities from E κ and E O -despite the difficulty with determining their significance.

350
To investigate whether the results are sensitive to our assumptions about the signal-to-noise ratio of our data-through W in Eq. 1-we additionally perform Monte Carlo simulations for the adjoint sensitivities from E κ -using three different data sets for κ ρ : κ ρ,micro , κ ρ,Argo together with κ ρ,CT D , and κ ρ,tides . In the Monte Carlo simulations, at each location and depth, we randomly sample κ ρ values within its uncertainty-σ κ -simultaneously with randomly sampled values of σ κ between a factor of 2-3 of κ ρ . This accounts for the uncertainties in the κ ρ products in both the numerator-corresponding to uncertainties in the 355 observationally-derived estimates-and denominator of Eq (2)-corresponding to the weights. With each of the 10,000 samples of κ ρ and σ κ , we recompute the adjoint sensitivity for E κ and then its correlation with that for E O . With these Monte Carlo simulations, we report the maximum possible correlation between each experiment's adjoint sensitivities.

Results
We first show that the disagreements between κ ρ from ECCO and κ ρ from various observations are larger than the observations' 360 approximate 95% confidence intervals. Then we analyze results from pairs of adjoint sensitivity runs: one with misfits to observed κ ρ derived from the finescale parameterization and the other with misfits to observed O 2 . We use these results to investigate the potential to use O 2 as a constraint for improving κ ρ,ECCO in a future optimization. We then compare the results of the adjoint sensitivity runs using misfits in κ ρ with ones using misfits in to infer a potential role of stratification in any information that O 2 provides about κ ρ . In the Appendix, we show there is general agreement between κ ρ from observations and 365 a free-running earth system model that calculates a physically-motivated parameterization for κ ρ , but poor agreement between κ ρ from observations and a sequential ocean data assimilation system based on the same earth system model ( Fig. A1a).

Model-inverted vs microstructure-inferred κ ρ comparisons
We compare the average κ ρ,micro profile that is comprised of 24 campaigns' worth of data (Waterhouse et al., 2014) (see their The average κ ρ,ECCO,0 profile-i.e., the initial guess of κ ρ,ECCO -is typically smaller than the microstructure profile, particularly at 1000 m where the difference is approximately an order of magnitude ( Fig. 4a). At iteration 59 (which is the EC-

375
COv4r3 solution), the difference between κ ρ,ECCO and κ ρ,micro decreases. However, agreement between the average profiles of κ ρ,ECCO and κ ρ,micro is still worse than the agreement between κ ρ,Argo and κ ρ,micro . The agreement between κ ρ,Argo and κ ρ,micro at each of the three depth bins is well within a factor of three (dotted black curves in Fig. 4a) and the spatial standard deviation of κ ρ,micro (dashed black curves in Fig. 4a). The agreement between the average profiles of κ ρ,ECCO and κ ρ,tides is poor, with κ ρ,ECCO typically too small, notably so at deeper depths ( Fig. 4b). κ ρ,ECCO includes internal-wave-induced mix-380 ing as well as potentially numerical diffusion. However, numerical diffusion cannot explain the errors in κ ρ,ECCO where κ ρ is too small in the model relative to the observationally-derived products because numerical diffusion would increase κ ρ,ECCO .
In these regions, one likely explanation is that errors in other model parameters (e.g., the Redi coefficients) compensate for the errors in κ ρ,ECCO .
We also compare κ ρ,ECCO and κ ρ,ECCO,0 profiles with κ ρ,micro from 16 example campaigns in Fig. 5. In some regions, 385 the κ ρ,ECCO and κ ρ,ECCO,0 profiles are constant (10 −5 m 2 s −1 , the default background value) because ECCO does not sufficiently resolve the bathymetry so we exclude those from Fig. 5. We also exclude some others, for example, in the subpolar North Atlantic because temporal variations in κ ρ can be large there (Fig. A1b). The κ ρ,ECCO profiles (blue curves in Fig. 5) and κ ρ,ECCO,0 profiles (grey curves in Fig. 5) are often within the approximate (factor of three) uncertainties in the κ ρ,micro profiles (dashed black curves in Fig. 5), but not always. Without taking an average over all of the campaigns, there can be 390 large regional disagreements between the model and observations. Also, the κ ρ,ECCO profiles are not always closer to the κ ρ,micro profiles than the κ ρ,ECCO,0 profiles. This suggests that performing more iterations of the optimization of ECCO is not necessarily going to lead to more accurate representation of κ ρ with the current data constraints.

Model-inverted vs finescale parameterization-derived κ ρ comparisons
We next show κ ρ,Argo and κ ρ,tides as well as how they contrast with κ ρ,ECCO because this highlights the spatial patterns of many under-determined parameters in ECCO could also lead to erroneous κ ρ,ECCO . These are some reasons why values of κ ρ,ECCO do not agree well with κ ρ from observations-κ ρ,micro , κ ρ,Argo , κ ρ,CT D , or κ ρ,tides .

Adjoint sensitivities in ECCO
Because the data that are currently included in ECCO's misfits are insufficient for κ ρ,ECCO to match κ ρ,micro , κ ρ,Argo , κ ρ,CT D , or κ ρ,tides , including additional variables controlled by mixing in the model's misfits may assist in further improving 420 the modeled mixing parameters. Oxygen is a candidate since its distribution is, in part, determined by the local κ ρ . To test this, we run multiple adjoint sensitivity experiments in which either observationally-derived κ ρ or oxygen is included in the misfit calculation to guide constraints on κ ρ . We expect that the signs of sensitivities agree most in regions away from where air-sea fluxes and transport of oxygen-e.g., by intensified jets-are large. One of these regions is the subtropical North Atlantic Ocean, away from the Gulf Stream Extension. Further, we expect to find more agreement between the signs of sensitivities in tropical 425 OMZs and other (sub)tropical regions because of the known importance of diapycnal mixing in these regions.
We show the adjoint sensitivity calculations using Eq. 2 for κ ρ misfits (experiment E κ in Table 2 6b,d,f-aren't representative of the ocean where Argo measurements were taken. 440 We next compare ∂J/∂ κ ρ from E κ using κ ρ,Argo and κ ρ,

450
We also show the adjoint sensitivity calculations using Eq. 2 for κ ρ misfits (experiment E κ in Table 2) in Fig. 10 using κ ρ,tides and compare ∂J/∂ κ ρ from E κ using κ ρ,tides with ∂J/∂ κ ρ from E O . The signs of ∂J/∂ κ ρ using κ ρ,tides ( where κ ρ,ECCO needs to be increased become dominant closer to the seafloor, particularly in the Atlantic Ocean. The signs 10a,c and Figs. 10b,d) using κ ρ,tides instead of κ ρ,Argo and κ ρ,CT D . The regions with agreement in signs of sensitivities using κ ρ,tides account for just over half of the ocean's volume where they can be compared globally ( Fig. 11; Table 4); this is also true for the subtropics. The equatorial regions have 460 the highest percent volume of agreement in signs of sensitivities using κ ρ,tides over all depths, but the North Atlantic also has fairly high agreement (Fig. 11a). There is high agreement in the regions of the Arctic Ocean that are north of the Greenland and Barents Seas too. Compared with shallower depths, regions below 3000 meters depth tend to be derived from Antarctic Bottom Water (Marshall and Speer, 2012 -see their Figure 1) and therefore have different oxygen concentration characteristics such as weaker vertical gradients, have differences between κ ρ,ECCO and observationally-derived κ ρ that are more commonly 465 statistically indistinguishable, and have less overwhelmingly positive adjoint sensitivities from E κ using κ ρ,tides (Fig. 10c).
As a result, all of the depths with the highest levels of agreement in signs of sensitivities using κ ρ,tides are between the mixed layer depth and 3000 meters depth (Fig. 11b). Most differences in the spatial distribution of agreements between the signs of sensitivities across different observationally-derived products (Figs. 9 and 11) are at least partially due to their spatial coverage-Argo versus global-and the < 100% overlap in processes accounted for by the various κ ρ products. Thus, the level 470 of agreement in signs of sensitivities from E κ and E O is high over many regions, and is qualitatively consistent in its spatial distribution across the different observationally-derived products for κ ρ .
We need to address whether any of the agreement in signs of sensitivities is random-as their correlation is due to the large uncertainties in observationally-derived κ ρ -or underpinned by physical reasons. We first focus on the locations with statistically indistinguishable errors in κ ρ,ECCO . These regions and those where there can be significant differences between 475 oxygen concentrations in ECCO and the World Ocean Atlas (2013) product correspond to the white regions in Fig. 9 that are red or blue in Fig. 8-likewise for Fig. 11 versus Fig. 10. The vast majority of the locations where disagreements occur in sensitivity signs are in places with statistically indistinguishable differences between κ ρ,ECCO and observationally-derived κ ρ . The regions with statistically indistinguishable differences in κ ρ account for 56.2% (24.1%) of the volume of the ocean where the adjoint sensitivities from E O and E κ can be compared using κ ρ,Argo and κ ρ,CT D (κ ρ,tides ). Thus, we exclude a 480 large portion of the ocean from the remainder of our analysis because we cannot determine whether agreements in signs of sensitivities are by random chance in these regions.
We next inspect the sensitivity sign patterns in regions with statistically significant κ ρ misfits. The regions where the signs of ∂J/∂ κ ρ agree from the two experiments and have large differences between κ ρ,ECCO and the combined κ ρ,Argo and κ ρ,CT D product tend to have relatively small oxygen concentration misfits (Fig. 3b). This is also true when using the κ ρ,tides product.

485
For example, when only regions with less than one standard deviations above the average oxygen concentration misfits are selected, the signs of the adjoint sensitivities agree between E O and E κ over 60.8% (50.5%) of the volume with sufficient data using κ ρ,Argo and κ ρ,CT D (κ ρ,tides ). However, the larger the oxygen concentration misfits, the more often the signs of the sensitivities agree. When only regions with more than three standard deviations above the average oxygen concentration misfits are selected, the signs of the sensitivities always agree. Thus, the regions with the largest disagreements in oxygen 490 concentrations can always decrease their oxygen misfits by changing κ ρ,ECCO with a sign consistent with decreasing its disagreement with observationally-derived κ ρ , wherever differences in κ ρ are detectable.
Where there are statistically significant differences in κ ρ , we still need to determine whether there is a physical basis for the agreements in signs of sensitivities. We next show results that are consistent with our hypothesis that the tropical OMZs and other ( We further address whether the potential information dissolved oxygen concentrations provide about κ ρ is due to the information oxygen contains about stratification. To determine whether oxygen provides information about stratification-and through stratification, about κ ρ -we use the adjoint sensitivity results obtained from experiment E with observationally-derived dissipation rates, = N 2 κ ρ /0.2 (e.g., Figs. 1b,d,f) instead of κ ρ (e.g., Figs. 6b,d,f), in the misfit function via Eq (2) and multi-515 ply the adjoint sensitivity of E O by 0.2/N 2 so that their sensitivities are each taken with respect to (parentheses in Table 4).
We find approximately equal agreement between the signs of the adjoint sensitivities from E O (scaled by 0.2/N 2 ) and E as we do between those from E O and E κ in every region, regardless of which observationally-derived product we use. Because is related to κ ρ through the stratification, we suggest that the information oxygen concentrations provide about κ ρ is likely independent of the stratification field.

520
Lastly, given that the general agreement in signs of sensitivities between E κ and E O are likely underpinned by physical reasons unrelated to stratification, we pursue whether there is a statistically significant relationship between the adjoint sensitivities from E κ and E O . We again focus on the regions where the difference between κ ρ,ECCO and observational κ ρ products (from tides, Argo/CTD, and microstructure) is statistically significant (greater than a factor of three), but also filter out the adjoint sensitivities where the differences between oxygen concentrations from ECCO and those from the World Ocean Atlas (2013) 525 are statistically significant. The simple correlations between the adjoint sensitivities from E κ and E O in the remaining regions tend to be small but positive (Fig. 12). In addition to taking simple correlations, we've performed that the information oxygen provides about κ ρ is not conditional on the stratification. This suggests that κ ρ,ECCO may be con-535 strained by the information provided by oxygen concentrations. That is, oxygen concentrations inform adjoint sensitivities that typically direct κ ρ,ECCO to improve relative to observationally-derived κ ρ . However, inclusion of accurately known oxygen concentrations in the model's misfits is not a perfect substitute for the inclusion of accurately known κ ρ itself in the model's misfits.

Discussion
This study evaluated the potential to improve the diapycnal diffusivities (κ ρ ) in the ECCOv4 ocean parameter and state estimation framework. We assessed the fidelity of the inverted field of κ ρ,ECCO by first comparing the average inverted vertical profiles of κ ρ,ECCO with those inferred from microstructure. The comparison was not favorable. κ ρ,ECCO is inverted for within the ECCO framework through constraints of vertical profiles of temperature and salinity-e.g., from Argo profiles.

545
Model choices-e.g., the initial guess of κ ρ = 10 −5 m 2 s −1 everywhere-can lead to errors in κ ρ,ECCO even in the presence of globally complete hydrographic observations (see Section 4.2), but we investigated whether κ ρ,ECCO can benefit from new information.
We then investigated which additional observations can be used as new constraints to improve the fidelity of the inverted κ ρ,ECCO . The products we used were observationally-derived κ ρ based on Argo and ship-based CTD hydrographic data, 550 observationally-derived κ ρ based on climatological and seafloor data, and climatological oxygen concentrations. To justify the use of the observationally-derived κ ρ products, we also evaluated them by comparing them with the microstructure-inferred product. κ ρ,Argo and κ ρ,CT D have better agreement with the microstructure-inferred data than κ ρ,ECCO does.
We inspected the misfit of the model parameter κ ρ,ECCO with respect to κ ρ,Argo and κ ρ,CT D as well as to κ ρ,tides and motivated use of dissolved oxygen concentration data as a potential constraint in ECCO. One drawback of the observationally-555 derived data products for κ ρ is that they have large uncertainties-here, approximated by a factor of three. Observed oxygen concentrations, on the other hand, have a relatively small uncertainties. More importantly, we showed that vertical oxygen gradients have similar geographical patterns to energy dissipation rates. We therefore performed an additional adjoint sensitivity experiment with oxygen concentration data as the only data in the misfit function. Adjoint sensitivities results were compared between the experiment with measured oxygen in the misfit function and observationally-derived κ ρ in the misfit function.

Caveats and future directions
Many factors-including a significant absence of independent observations for assessment, a combination of measurement 570 and structural errors, numerical diffusion in our simulations, and unconstrained parameters in the biogeochemical moduleshave stymied progress in state estimation of ocean subgrid-scale transport and mixing parameters. First, the observational measurement errors used here are only approximate. We assumed uncertainties equal to a factor of three of the observationallyderived κ ρ and 2% of the oxygen concentrations. These do not account for interpolation/averaging errors that entered the data prior to our calculations, but are conservative estimates nonetheless. The observational uncertainties affect the weights given in 575 the misfits that enter the adjoint sensitivity calculations and our Monte Carlo simulations of the correlations between the adjoint sensitivities account for the possibility that these weights are misspecified. Second, only one ocean subgrid-scale transport or mixing parameter-namely, κ ρ -has been compared with independent observational data-microstructure. This is the primary reason why we focused on κ ρ in our study. Third, the ECCO-estimated κ ρ accounts for other model errors-e.g., structural ones suggested by Polzin et al. (2014)-which explains some of the model biases relative to microstructure observations. For 580 instance, the ECCO-estimated κ ρ should be time-dependent as well as spatially-varying, but it is only spatially-varying. In the presence of other estimated model parameters and initial conditions, some parameters could be compensating for errors in κ ρ .
The ECCO-estimated κ ρ can also be sensitive to the a-priori estimate of κ ρ and we showed how one particular initial guess-10 −5 m 2 s −1 everywhere-can evolve from the first optimization iteration to the final one. Additionally, there is numerical diffusion in the model, which could confound some physical inferences about the model-e.g., regarding how sensitive the 585 model's state is to κ ρ relative to along-isopycnal diffusion. Numerical errors could remain and result in the primary source of error in the ocean state estimate even if additional constraints are placed on κ ρ in ECCO. Retaining a parameter which absorbs structural model errors that are not expressed in the model's functional form may be necessary to improving the ECCO state estimate itself, but minimizing numerical errors would be beneficial to improving the ECCO-estimated κ ρ . Lastly, there are several unconstrained parameters in biogeochemical modules used to calculate biogeochemical tracers (Verdy and Mazloff ,590 2017), so some of the disagreements in signs of the adjoint sensitivities found here could be associated with other inaccurate parameters.
These challenges can continue to be overcome by allowing models and observations to inform each other. First, the observationallyderived κ ρ from the finescale parameterization could be further scrutinized using ship-based CTD data taken concurrently with microstructure velocity shear data. A preliminary analysis suggests that the percent difference between the full depth-averaged 595 CTD-derived κ ρ from the finescale parameterization and the microstructure-inferred κ ρ at the same locations is 1.68%, which is indistinguishable from zero, but the quality of the the CTD data taken concomitantly with microstructure has not been fully assessed. Second, we would need to account for the time-dependence of κ ρ in a future ocean state estimate. The underdetermined nature of the parameter estimation procedure makes this difficult. These efforts would also benefit from minimizing spurious mixing due to numerical diffusion (e.g., Holmes et al. (2021) informative of κ ρ than in other locations. We did not pursue this in the present study because our adjoint runs use a global misfit. If we perform an ensemble of adjoint sensitivity runs with a single observation in each run, then we could calculate 615 the effective proxy potential at each of these observation locations. Lastly, the (imperfectly-known) initial conditions of biogeochemical tracers will also need to be included in the input control vector during optimization of the ocean state estimate.
If biogeochemical tracers are included in the misfit calculation in an optimization run, their impact on variables such as κ ρ would depend upon how they are weighted relative to the physical variables-e.g., temperature, salinity, and pressure. A more complete representation and understanding of κ ρ is possible through these analyses and methods.

620
Appendix A: Model with a sequential data assimilation framework

A1 GEOS-5 and the GMAO S2S Ocean Analysis
To demonstrate a problem with κ ρ in a sequential data assimilation framework, we present example output from a reanalysis product and output from an identical ocean model hindcast without any data assimilation. GEOS-5 includes a global, finite volume atmospheric general circulation model that is used for numerical weather prediction, seasonal-to-decadal forecasts,

A2 Steric sea level budget framework
In order to examine whether the analysis increments can dynamically impact κ ρ , we analyze a model's buoyancy budget, which is broken down into heat and salt budgets and used to calculate the steric sea level budget. The tracer tendency equation terms required for the heat and salt budgets were computed online and saved as the reanalysis was running. The tracer equations can 690 be broken down into individual contributions (Palter et al., 2014), where d/dt = ∂/∂t + (v + v * ) · ∇ is the material derivative, v is the resolved velocity field, v * is the eddy-induced or quasi-Stokes velocity field that represents parameterized motions, Θ is the potential temperature, S is the salinity, ρ is the locally B3 Relationships between steric sea level budget terms There are distortions in temperature and salinity fields from applying analysis increments, violating conservation principles and potentially causing the model to undergo baroclinic adjustment (Stammer et al., 2016). Thus, we examine whether the velocity field itself changes because of the analysis increments. To do this, we show the relationship between the analysis increments and resolved advection terms in the steric sea level budget for the GMAO S2S Ocean Analysis in Fig. A3a.     ) (tides) (tides) Figure 12. Shown are scatterplots between the adjoint sensitivities from EO and Eκ where they are both negative (panels a, c, and e) and where they are both positive (panels b, d, and f), where Eκ has its adjoint sensitivities calculated with: the tidal mixing product κ ρ,tides (panels a-b), Argo-derived κρ (κρ,Argo and κρ,CT D ; panels c-d) or microstructure-inferred κρ (panels e-f). Only the adjoint sensitivities where the differences between κρ,ECCO and observational κρ products are statistically significant (greater than a factor of three) and where the differences between oxygen concentrations from ECCO and those from the World Ocean Atlas (2013) Table 2. Listed are the ECCO simulations performed and analyzed in the present study as well as the observationally-derived data or measured data included in each simulation. Either observationally-derived data or measured data are included in the experiments through its misfit calculation (Eq. 1). Here, κ ρ,obs denotes an observationally-derived κρ product derived from a parameterization (κρ,Argo and κρ,CT D or κ ρ,tides ) or inferred from microstructure (κρ,micro), = κρN 2 /0.2 indicates an observationally-derived dissipation rate (N 2 is the stratification from the World Ocean Atlas or WOA (2013)), and O2 is the climatology of measured oxygen concentrations from WOA  Table 3. The control variables that ECCO inverts for and optimizes. Some of these control variables are initial conditions only (indicated with the "initial condition" column). Other control variables are time-varying (indicated with the "time-varying" column). The rest are not initial conditions, but also are time-independent. Also noted is whether the control variable's field is two-dimensional or three-dimensional-there are no control variables that vary in both time and over all locations and depths of the ocean. Redi coefficients (Redi, 1982) no no 3 Gent-McWilliams coefficients (Gent and McWilliams, 1990) no no 3 background κρ no no 3 surface forcing fields no yes 3 Table 4. Listed are the percent volumes where the signs of the adjoint sensitivities agree between Eκ-for both the Whalen et al. (2015) and where sufficient observations are available to derive κρ and where the difference between the model-calculated and observationally-derived κρ is greater than the uncertainty (i.e., three times the observationally-derived κρ). In parentheses are the same, except for the dissipation rates, ρ = N 2 κρ/0.2, where N 2 is the stratification and 0.2 is an empirical coefficient (see, e.g., Gregg et al. (2018)).