Using the German Estimating the Circulation and Climate of the Ocean (GECCO) synthesis framework, four separate eddy tracer mixing coefficients are adjusted jointly with external forcing fields, such as to reduce a global misfit between the model simulations and ocean observations over a single 10-yr period and weighted by uncertainties. The suite of the adjusted eddy tracer mixing coefficients includes the vertical diffusivity kz, the along-isopycnal surface diffusivity kredi, the isopycnal layer thickness diffusivity kgm, and the along-iso-thickness advection coefficient kgmskew. Large and geographically varying adjustments are found in all four parameters, which all together lead to an additional 10% reduction of the total cost function, as compared to using only surface flux parameters. However, their relative contribution to the cost reduction varies from 1% to 50% among the four coefficients, with the adjusted kgm contributing most. Regionally, the estimated kgm ranges from less than −800 to about 2500 m2 s−1. Largest adjustments in kgm reside in the vicinity of large isopycnal slopes and support a mixing length hypothesis; they also likewise support the hypothesis of a critical layer enhancement and high potential density gradient suppression. In a few occasions, resulting negative net kgm values can be found in the core of main currents, suggesting the potential for an inverse energy cascade transfer there. Large adjustments of kredi and kgmskew are found in the vicinity of isopycnal slopes. The adjustments of kz in the tropical thermoclines suggest deficiencies of the mixed layer parameterization.
Water mass properties can be greatly influenced by diffusion and redistribution of water masses induced by eddies. Physically most appropriate is to think of eddy-induced tracer mixings as processes occurring along (Redi 1982) and across isopycnal surfaces. In coarse-resolution models these processes are not explicitly resolved; instead, the associated subgrid-scale eddy-induced mixing process needs to be parameterized. Gent and McWilliams (1990) and Gent et al. (1995) presented a parameterization scheme (hereafter GM) suggesting that eddy mixing can be regarded as diffusion of isopycnal layer thickness, which can be treated as an extra advection of tracers across the isopycnal slopes. More recently, Eden et al. (2007) suggested an additional skew part to GM, another advection process of tracer along rather than across the isopycnal slopes. Following the above parameterization schemes, the along-isopycnal diffusivity kredi, the across-isopycnal diffusivity kz, the thickness diffusivity kgm, and its skew counterpart kgmskew have to be specified. With those terms included, the equation for any tracer T can be written as
Here, sources or sinks were neglected, and kz is the coefficient of the vertical diffusion, which can be considered as diapycnal mixing. Here, u and u* are Euler and eddy (or bolus) velocities, respectively, and is given as
which are the parameterizations specified by GM and Eden et al. (2007), respectively, and are perpendicular to each other. They are parameterized by the isopycnal slopes and advect mass properties. The vector denotes the vector of isopycnal slopes, and the tensor
projects gradients of T onto the isopycnal surface. In these parameterization schemes, the small-slope approximation |s| ≪ 1 is invoked because a typical slope in the ocean interior is only of the order of 10−4. However, slopes in regions such as deep convection or the surface mixed layer can be (infinitely) large, which leads to large variations in the tensor elements and hence large bolus velocities, resulting in instability of numerical modeling. This was solved by artificially reducing the slopes with “tapering” functions (e.g., Large et al. 1997).
Generally, kgm should be positive so that (on global average) the mean potential energy is converted to eddy kinetic energy. However, this requirement does not need to hold locally as suggested by eddy-resolving models, which showed energy conversion in the opposite direction in western boundary currents (WBCs; Eden et al. 2007). Evidence that there is an inverse cascade eddy transfer, from eddy kinetic energy to large scales and up-gradient potential vorticity (PV) fluxes, both leading to negative kgm coefficients, was also provided by Venaille et al. (2011) and Waterman and Jayne (2011).
Estimates from both observations and numerical studies have found a wide range of eddy diffusivity amplitudes and an evident spatially varying pattern (Rix and Willebrand 1996; Stammer 1998; Bryan et al. 1999; Treguier 1999; Nakamura and Chao 2000; Roberts and Marshall 2000; Drijfhout and Hazeleger 2001; Solovev et al. 2002; Zhurbas and Oh 2004), which is consistent with Gent and McWilliams’s (1990) argument that it may be important to make kgm a function of space and time. This is particularly true for the vertical diffusivity kz (e.g., Canuto et al. 2004).
In previous studies, kgm and kredi were often set to O(1000) m2 s−1 while kz was O(10−5) m2 s−1, which are consistent with scaling arguments and observations. In contrast, only little is known about details of kgmskew. For example, Eden (2006) and Eden et al. (2007) diagnosed a wide range of values in the Southern and Atlantic Oceans, whereas from idealized model experiments Eden (2010) found that kgmskew should vanish or at least be very small. This is the first of the two motivations to estimate a global map of the mixing coefficients, the second follows.
The second motivation comes from the successful skill of the Estimating the Circulation and Climate of the Ocean (ECCO) synthesis in producing reasonable external forcing parameters and optimized ocean state (Stammer et al. 2002a,b, 2003; Wunsch and Heimbach 2006; Köhl et al. 2007; Köhl and Stammer 2008a,b; Romanova et al. 2010). In the ECCO approach, an ocean circulation model is being fit to observations (within error limits) in a dynamically consistent manner by adjusting uncertain model parameters. The approach is based on the Massachusetts Institute of Technology GCM (MITgcm; Marshall et al. 1997a,b) and its adjoint used to calculate the sensitivity of the ocean model output to the control parameters. The “adjoint model” is obtained from the tangent linear model of the MITgcm by replacing the resulting linear operator with its adjoint.
Most of the previous ocean state estimation work (Stammer et al. 2002a,b, 2003; Wunsch and Heimbach 2006; Köhl et al. 2007; Köhl and Stammer 2008a,b) made the assumption that the model uncertainties reside entirely in the initial conditions and the surface forcing fields. Although the ocean estimations from the ECCO efforts have moved considerably closer to the observations, certain misfits remain. For example, the simulated ocean state is not yet fully consistent with all data types and prior error assumptions (Köhl et al. 2007), nor is it consistent in detail with ocean transports from independent simulations (Stammer et al. 2003). Remaining differences in temperature and salinity are associated with errors in the vertical diffusion and the poor representation of the overflow physics, which points to the necessity of including an overflow parameterization and of including the mixed layer model parameters as controls into the optimization (Köhl et al. 2007). Although the optimization leads to a more realistic ocean state, unrealistic surface conditions become obvious at certain places. For example, the constrained model reproduces a better separation of the Gulf Stream by shifting its axis farther south at the expense of large forcing adjustments, which unphysically compensates for the model deficiencies (Stammer 2005).
In principle, errors in one model parameter can be compensated by adjusting other uncertain parameters. This inevitably can lead to a projection of errors on a large number of parameters. How much this is an issue during the optimization depends on the sensitivities of model biases to certain parameters and their prior error weights. Because of the prior information of the parameters in the optimization procedure, parameter changes are penalized and error projection is not completely avoidable, even in identical twin experiments where errors reside only in one estimated parameter. However, a more complete set of uncertain parameters should improve the optimization result. In an attempt to expand the ECCO approach to also estimate internal model parameters, Stammer (2005) estimated horizontal and vertical viscosity and diffusivity fields together with surface forcing fields; Ferreira et al. (2005) estimated an eddy stress field by formulating the velocity in terms of the residual mean.
The goal of the present work is to apply the German ECCO (GECCO) framework to obtain estimates of four eddy tracer mixing coefficients. As this first step, they are assumed to be spatially dependent but constant in time (in future efforts they might become temporally varying as well).
The structure of the remaining paper is as follows: In section 2, we will give a brief introduction of GECCO setups and the techniques required to guarantee the numerical stability of the models. The contributions from the adjusted mixing and external forcing parameters to the improvement of the ocean simulation are presented in section 3. Estimates of the four mixing parameters and their implications to the mixing theories are presented in section 4. Discussion and conclusions follow in section 5.
2. Estimation methodology
a. GECCO setups
GECCO is the German version of ECCO and was described in detail by Köhl et al. (2007). Its setup can be summarized here as the MITgcm with a horizontal resolution of 1°, a meridional coverage from 80°S to 80°N, and 23 vertical levels. The GM eddy parameterization and the K-profile parameterization (KPP) surface mixed layer scheme (Large et al. 1994) are used. Background values of kgm, kgmskew, kredi, and kz are 800, 0, 800, and 10−5 m2 s−1, respectively, which are consistent in magnitude with the previous work. In addition, a background horizontal Laplacian diffusivity of 102 m2 s−1 is imposed. The estimation window is 10 yr (1992–2001).
In the first forward run of the optimization procedure, which is called GECCO_Reference run (GR), the MITgcm is initialized with climatological potential temperature and salinity fields of Levitus (Levitus and Boyer 1994; Levitus et al. 1994) and is driven by the 6 hourly National Centers for Environmental Prediction (NCEP) atmospheric states (Kalnay et al. 1996). The surface fluxes of heat, salinity, and momentum are calculated based on the atmospheric state variables.
During every forward run over the 10-yr period, the misfit between the model simulation and ocean observations is measured by the cost function. The backward run of the adjoint model, which is driven by the misfits, follows when the forward run finishes and provides the sensitivities of the cost function to the control parameters (Marotzke et al. 1999). The control parameters are adjusted accordingly and applied in a new forward run. The loop as described above is performed over several iterations until the finally adjusted control parameters produce an optimized state that is consistent with the observations (the minimum of the cost function). The general methodologies are described by Wunsch (1996), Malanotte-Rizzoli (1996), and Fukumori (2001); details of ECCO approach are provided by Stammer et al. (2002b).
The cost function is constructed from data constraints (observations), background values of the control variables, and their corresponding prior errors. The suite of control parameters includes 12 variables, notably the initial potential temperature and salinity, air temperature, air humidity, net precipitation, longwave radiation, two components of wind speed, and four mixing parameters mentioned in the above section. We note that the contribution of the four mixing control variables to the total cost function is quite small compared to their time-varying counterparts. The data constrains are the same as in Köhl et al. (2007) and include a wide range of datasets, which include sea surface height (SSH), sea surface temperature (SST), sea surface salinity (SSS), surface drifter velocities, and in situ hydrographic temperatures and salinity (some of which could be seen from Table 2). The Levitus error profiles ranging from 0.5°C (0.13) near the surface to 0.05°C (0.01) at depth are used as uncertainties for the initial potential temperature (salinity). The standard deviations of the NCEP atmospheric state variables are used as the corresponding uncertainties. Patterns of those uncertainties can be found in Fig. 4 of Stammer et al. (2002b) and in Köhl et al. (2007). The uncertainties for the mixing coefficients kgm, kgmskew, kredi, and kz are set to 1000, 1000, 1000, and 10−6 m2 s−1, respectively.
To reduce the misfits and to bring the model into consistency with the observations, we assume that the model uncertainties reside in the initial conditions, the surface forcing fields, and in the mixing parameters. A justification for this assumption was provided by Stammer (2002b) and Köhl et al. (2007), based on the fact that these parameters are known to have large uncertainties and once they were adjusted (still within the data uncertainties) they can reproduce a solution that better fits the observations.
Hereafter, we will denote the two initial condition fields and six surface forcing fields as forcings and refer to the four mixing coefficients as mixings. The run with the finally estimated parameters is called GECCO_Forcing&Mixing (GFM). Using subsets of the 12 estimated parameters, we conducted five additional simulations to test the effects of forcings, mixings, and the four mixing coefficients individually, which will be discussed in section 3. The setups and the naming of the individual experiments are listed in Table 1.
b. Dealing with strong sensitivity in the adjoint model
The use of the adjoint method in applications with nonlinear general circulation models can be problematic, because the nonlinear dynamics of the model may lead to cost functions that do not converge to a global minimum because of multiple local minima (Pires et al. 1996; Li 1991). This is particularly true for the strong nonlinear GM and KPP schemes, so in the previous work the adjoint code to these schemes had to be excluded from the adjoint because these codes otherwise would cause exponentially growing gradients of the cost function (Köhl and Stammer 2008b).
Two main approaches can be employed to suppress the exponentially growing sensitivities. One is to reduce the model resolution or to increase the diffusivities in the adjoint model (Köhl and Willebrand 2002; Hoteit et al. 2005); the other is to modify the original adjoint model by excluding the unstable modules (Zhu et al. 2002; Weaver et al. 2003). Regardless of which approach is applied, the forward model remains unchanged and contains the full physics.
In the present study, following Zhu et al. (2002), we exclude from the adjoint to GM/Redi only the part that relates to the linearization of mixing coefficients with respect to the state variables; the remaining code enables the computation of gradients with respect to mixing parameters. To test the validity of the modified adjoint, we carried out two experiments with a 2-yr estimation period over which both the original and modified adjoints to GM are stable. Results (Fig. 1) suggest that the modified adjoint preserves the large-scale pattern of the gradient of cost function with respect to kgm in both the sign and the magnitude in regions that are dominated by strong density fronts, especially in the WBCs and the Antarctic Circumpolar Current (ACC) regions. A notable exception is the area centered around 20°S and 120°W, where the modified adjoint does not produce the positive pattern of the original adjoint; in addition, the magnitude of the gradient in the modified adjoint is a bit larger than that in the original adjoint in the ACC region.
3. Evaluation of the optimization
a. Cost functions
Results presented below are based on iteration 25; at that number of iterations, any further reduction in the cost function became insignificant. By then, the total cost function decreased by 54% relative to the first iteration; distributed over individual constraints, the misfits decreased between 6% and 90%, as summarized in Table 2. A major reduction in the cost function results from a better fit of the model to the time-varying temperature and salinity data from Argo (reduced by about 70%–80%) and to SST and SSS observations (reduced by about 65%–70%). The initial misfits in the climatological temperature and salinity were reduced by 43%–55%.
Overall, about 90% percent of the reduction in the total cost function originates from optimized forcings. The patterns of the estimated initial conditions and the atmospheric fluxes (not shown) are similar to the estimates in Stammer et al. (2004); however, the magnitudes of the adjustments are smaller, implying that the changes of forcings are still within the error bars. The adjustment of mixings leads to an additional 10% reduction in the total (global) cost function. However, the contribution of adjusted mixings to the reduction of the individual data misfits varies largely geographically and can be much higher locally, especially for time-mean quantities. For example, they contribute about 20% to the time-varying temperature in the CTD and Argo data and about 50% to the time-mean SSH field though only about 10% to the temperature and salinity climatologies. Among the four mixing parameters, the kgm parameter contributes most to the reduction of the total cost function and to most of its contributions as well, kgmskew and kredi follow, and kz contributes least. The comparison is provided with Table 2.
We note that all of the parameters kgm, kgmskew, and kredi are adjusted by changes in the order of the same error bounds after 25 iterations. The fact that kgm nevertheless contributes most to the improvement implies that a larger part of the initial model–data misfit points to the necessity to use spatially varying kgm. However, changes of kz, are limited the magnitude and region where they occur because of the relatively small a priori uncertainty and the short estimation time period (10 yr). Significant changes in kz occur only in the upper ocean, especially in low latitudes, where mixing time scales are shorter, which explains the small contribution of kz, in the midlatitudes and at depths here. Additional experiments (not shown) with longer estimation period and larger weights show that the structure of estimated kz and its effects are very different from the present results, especially below the upper ocean and in the extratropics.
b. The mean temperature fields
Differences of the 10-yr-mean potential temperature between the individual experiments and the Levitus temperature data are presented as root-mean-squares (RMS) in Fig. 2. Improvements are found over the whole water depth, especially between 100 and 1000 m, where the RMS decreases by up to 0.3°C. In the levels between 300 and 1000 m, which cover the main thermocline, both forcings and mixings contribute significantly to the decrease in RMS, whereas between 300 and 500 m mixings contributes about one-third of the reduction from forcings. The contribution from mixings in such levels certainly affects the structures of the thermocline and hence the ocean circulation. However, it is also obvious that in the top 100 m, where the tapering of the GM scheme renders the parameterization less efficient, and below 1000 m, where density slopes become flat and GM scheme does not take effect, that mixings parameters contribute little. Therefore, changes in the surface mixed layers and in the deep ocean are predominantly controlled by the other estimated parameters; namely, the surface fluxes fields and the initial conditions. The large remaining temperature RMS (>1°C) above 100 m, which is several factors larger than the prior uncertainties at these levels, may come from the incomplete adjustments of the mixing parameters due to tapering or incomplete adjustments of surface forcing fields, data errors, or prior errors, which we will not discuss here.
An example of differences in the 10-yr-mean temperature at a model level (310 m) of the main thermocline is shown in Fig. 3. The residuals between the nonoptimized experiment GR and Levitus (Fig. 3a) are large nearly over the entire ocean (the tropics are exceptions) with magnitudes around ±1.5°C. Very large residuals up to +5°/−4°C exist in the North Atlantic. The RMS in this depth level is 0.75°C. After the optimization, the residuals (Fig. 3b) are significantly decreased in the North Atlantic, including the Gulf Stream (by 3.5°C), because of adjustments in both forcings (Fig. 3c) and mixings (Fig. 3d). Improvements due to adjusted forcings are particularly also seen in the extension of Kuroshio, the global subtropical oceans, and some sections of ACC, whereas improvements caused by adjusted mixings are seen in the Kuroshio and its extension, over the Brazil Current, and in the entire ACC region. Nevertheless, large misfits remain in these areas. An evident feature is that improvements caused by adjusted mixings occur around the main currents and the maximum magnitude (4.75°C) of improvements is even larger than that by forcings (2.25°C) (Figs. 3c,d), which demonstrates regionally comparable contributions from mixings and forcings to the improvement of the model simulation.
Figure 4 shows that, among the four individual mixings parameters, kgm (Fig. 4a) contributes most to the improvement because both patterns and magnitudes are nearly the same as that from the joint mixings (Fig. 3d). Contributions from kgmskew and kredi (Figs. 4b,c) reside around the main currents, with the former a bit larger than the latter. Most of the improvements in the interior away from boundaries are contributed by the adjusted kz (Fig. 4d), especially in the tropical areas; changes in temperature caused by kgmskew are about half of the magnitude as that resulted from kgm.
We note that, in certain limited areas (e.g., around 40°N, 60°W over the Gulf Stream), impacts of the adjustments of kgm and kgmskew may counteract each other (Figs. 4a,b), leading to a smaller net effect than suggested by individual parameters. This is also true between the mixings’ and forcings’ effects, indicating that at the end it is a complex interaction of all parameters that leads to the final solutions and that all parameters need to be adjusted simultaneously.
For example, in the Kuroshio Extension region between 32° and 37°N and between 150° and 200°E, the adjusted forcings cause a slight increase of temperature while mixings lead to a bigger decrease (Figs. 3c,d). The most possible reason for the conflict is that one adjusted parameter reduces one constitute of the total cost function at the cost of increasing another constitute, which has to be reduced by other parameters, indicating the balance between different processes and also indicating that the final improvements should come from the joint effects of all the mixings and forcings. However, the conflict might be reduced to some extent in the future by the use of full error covariance matrices, which better reflect the actual degrees of freedom in the control variables and the space–time dependence of model–data differences.
c. Time-mean SSH
The difference between the mean Ocean Topography Experiment (TOPEX)/Poseidon (T/P) SSH observations (relative to the reference ellipsoid) and the University of Texas Gravity Recovery and Climate Experiment (GRACE) geoid model (GGM02; Tapley et al. 2003) is used here to constrain the mean dynamic SSH. Although the initial residuals between the nonoptimized experiment GR and the T/P observations (Fig. 5a) are smaller than ±10 cm over large parts of the ocean, residuals up to ±40 cm exist in the immediate vicinities of major boundary currents and in the high latitudes, indicating large departures of the unconstrained modeled currents from reality. These errors are improved during the optimization by adjusting both forcings and mixings (Fig. 5b). The adjusted forcings improve the SSH globally (Fig. 5c), with magnitude of improvements varying from 1 cm in the tropics through 5 cm in the subtropics to more than 20 cm in the main currents. Improvements caused by the adjusted mixings (Fig. 5d) reside generally over the main current regions, reaching there the same magnitude as improvements resulting from adjusted forcings, highlighting the importance of both parameters, as shown in Table 2. The improvements due to the adjusted kgm coefficient (not shown) are again quite similar to those of the net mixings, implying that kgm contributes most among the four mixing parameters.
Some overcorrections caused by adjusted forcings occur in the interior tropical Pacific, in some sections of ACC, and in some regions of the Gulf Stream (Fig. 5c), which are compensated by the adjusted mixings to some extent. Around the Gulf Stream and its extension, the adjusted forcings drives the Gulf Stream to separate from the slope at a more southward steering point indicated by the barotropic streamfunction (not shown) at the cost of a further deviation in the mean SSH between 35° and 45°N while improving the SSH in the Northwest Corner (NWC) (Fig. 5c); however, the adjusted mixings, while also producing a more consistent density field, draws the Gulf Stream slightly back to the coast indicated also by the barotropic streamfunction (not shown) and fits better with GGM02 without further modifying the NWC (Fig. 5d). The mixings change the density more than the forcings (which can be indicated by corresponding temperature changes in Figs. 3c,d) along the Gulf Stream path though increasing the misfit at the separation point. In the NWC, both mixings and forcings contribute and correct the density in the right direction; however, forcings are more important. In summary, the correct Gulf Stream path in the forcings experiment is more of a side effect of adjusted fields in the NWC, where the largest biases exist. However, it appears to be an overcorrection that is mitigated by the mixings, which further improves the SSH and density in the Gulf Stream.
Similarly, both mixings and forcings change the path of ACC in their own way. We note that in the Southern Ocean the changes of barotropic streamfunction due to both forcings and mixings show nearly the same pattern as those of SSH, except for the opposite signs of them, which enable us to check the ACC path change according to the SSH changes. For example, the adjusted forcings move the ACC path west of Drake Passage northward, corresponding with larger (negative) biases in the mean SSH (Fig. 5c). These biases, on the other hand, are mitigated to some extent by the adjusted mixings (mainly from kgm) (Fig. 5d), which meanwhile cause a southward movement of ACC path and a decrease of about 2 Sv (1 Sv ≡ 106 m3 s−1) in transport across Drake Passage. The (small) decrease in transport is accompanied with a spatially averaged increase of kgm in this region (Fig. 7). Note that both forcings and mixings improve the density here.
The above analyses indicate that by including mixing coefficients into the set of control parameters the optimization is able to more efficiently use the information from the sea level data. However, there might be conflicts between different data constraints, particularly that the mean SSH may have large errors in the WBCs (Rio et al. 2011).
The global (Euler) meridional overturning circulation (MOC) of the GFM experiment and the differences between the individual experiments are shown in Fig. 6. Near the surface, the forcings enhance the shallow tropical cells on both sides of the equator, decrease the overturning between 20°S and 60°N, and enhance the overturning north of 60°N and between 60° and 20°S (Fig. 6b). The maximum of the adjustments are slightly below 1000 m at 45°N and at 2000 m at 45°S, with amplitudes of about 8 and 5 Sv, respectively. MOC changes caused by mixings (Fig. 6c) are of the same pattern and reach up to 50% of forcings in the thermocline, with strengthened (weakened) overturning north (south) of 40°N above 1500 m. The Deacon cell above 2000 m is also enhanced. Changes of MOC caused by adjusted kgm (not shown) are again similar to the mixings’ effect, indicating that kgm contributes most to the joint mixings. However, the adjustments of mixings seem not to change the MOC much.
4. Estimated mixing parameters
Before interpreting the mixing parameters estimated in the experiment GFM, one has to answer the question of whether adjustments in control parameters reflect improvements of physical processes or merely reflect the correction of any other model deficiency; for example, lack of resolution or incomplete bottom topography. The former interpretation is supported if the estimated parameters are consistent with independent information (e.g., theories or observations) or at least if the adjustments match our understandings on errors of the concerned parameters. As an example, the estimated surface forcing fields in previous work have been well evaluated by comparing with independent estimates (Stammer et al. 2004; Romanova et al. 2010): buoyancy flux adjustments are found to be within the crude prior error bars and are consistent with known large-scale deficiencies in the NCEP products outside the boundary current regions, and wind stress adjustments are largely within the prior error bars and consistent well with satellite derived wind stresses. However, because of the limitation of observations and lack of knowledge about the spatial dependency of the mixing parameters, it is not easy to demonstrate whether the estimated mixing parameters are real or if they contain some compensation for deficiencies in representing other physical processes or model errors.
By adding here mixing parameters to the control vector in the present study (GFM), the estimated initial and surface forcing fields show even smaller adjustments compared to previous estimates without changing the mixing coefficients (Köhl et al. 2007), which means that the estimated forcings in the present study are still within the error bars. However, as an additional new feature, the model now computes the fluxes from the atmospheric state rather using surface fluxes from NCEP as was done previously, which also required the specification of different prior errors.
In the following, we describe our estimated mixing parameters while trying to evaluate the geographic patterns of them (especially of kgm) with those from theories and limited independent estimations from other models. We will also discuss whether the estimated parameters improve the physical processes in the model. Because the largest effects of adjusted mixings occur between 100 and 1000 m, the horizontal structures of mixing coefficients will be shown also at the model depth 310 m that represents the main thermocline [though in some limited regions (e.g., subpolar and mode water regions) this depth is seasonally still in the surface mixed layer].
a. Thickness diffusion kgm
1) Spatial structures of the kgm adjustments: δkgm
The adjustment δkgm to its background 800 m2 s−1 ranges roughly from −1500 to 1500 m2 s−1 (Fig. 7) and thus resides entirely within the limit (−5000 to 5000 m2 s−1) of previous estimates of eddy diffusivities from observations and models mentioned in the introduction. Positive (negative) adjustments imply that the background values (800 m2 s−1) were underestimated (overestimated) and should be increased (decreased). Generally, high absolute value of δkgm (more than 800 m2 s−1) occurs in the subpolar North Atlantic (except for the Labrador Current region), in the WBCs, and in the ACC regions; moderate absolute δkgm (between 100 and 500 m2 s−1) occur in the subtropical gyres and in the equatorial current systems, whereas small absolute δkgm (less than 100 m2 s−1) exit away from those regions. This general feature is consistent with other estimates (e.g., Eden 2006; Eden et al. 2007) and coincides with regions of high eddy–current interactions. The small adjustments in areas away from the strong and moderate currents result from the insensitivity of the model to kgm, which is discussed later, enabling us to consider that the real kgm in these regions may be quite small.
Generally, amplitudes of kgm are very small along the cores of the main currents or might even become negative (confined by the dashed lines in Fig. 7a where δkgm are less than −800 m2 s−1) in the cores of the Gulf Stream, Brazil Current, and Agulhas Current; in some regions of the ACC; and in the western coast of the North Atlantic. This kind of pattern extends from the bottom of the mixed layer down to about 1500-m depth along the ACC and the Gulf Stream (Fig. 8). Although this pattern is strengthening the previously too weak density fronts in the model to better fit the observations, it generally is plausible for geophysical flows that the mixing is weaker in the WBCs and stronger in the seaward flanks, as discussed by Shuckburgh and Haynes (2003). Weaker mixing along the main currents could be explained as the suppression effect of a strong current or a strong PV gradient on the eddy activities (Ferrari and Nikurashin 2010), which is discussed in more details below. Negative mixing diffusivities are also estimated in the main currents from a high-resolution model (Eden et al. 2007) and in the subtropical gyre in a low-resolution model (e.g., Ferreira et al. 2005). This indicates that the buoyancy (thickness) flux is up gradient and the available potential energy is fed by the eddy kinetic energy, which is demonstrated by evidences that modest transfer of energy to larger scales (an “inverse cascade”) exists (Venaille et al. 2011), and that eddies can generate an up-gradient PV flux to accelerate the mean current (Waterman and Jayne 2011). However, in the seaward flanks of the currents and in the region immediately north of the Kuroshio, southeast of the Gulf Stream and its extension, in the both flanks away from the ACC axis, and in other depth levels of the main currents (Fig. 8), kgm is positive. This guarantees a global sink of available potential energy.
From Fig. 7c, an apparent depth-decaying pattern of δkgm in the tropics and in the equatorward flanks of the subtropical gyres away from WBCs is consistent with previous findings (Ferreira et al. 2005; Olbers and Visbeck 2005; Phillips and Rintoul 2000; Ferrari and Polzin 2005). However, in the WBCs and in high latitudes, large δkgm occur in both the near surface layers and in the lower thermocline levels, which agrees well with Smith and Marshall (2009) and Abernathey et al. (2010), who found enhanced mixing in the mid depth in the ACC (between 1000 and 1500 m) in terms of effective diffusivity keff (Nakamura 1996). This pattern is often called the critical layer enhancement (e.g., Pierrehumbert 1990; Randel and Held 1991; Smith and Marshall 2009). Mixing is enhanced at the depth of the critical layer because the Rossby wave speed here nearly equals to the current velocity so that the eddy disturbances generated from the growing waves can efficiently displace fluid parcels to large meridional distances (Smith and Marshall 2009); this is also predicted by the linear instability theory (Killworth 1997) and confirmed by numerical experiments of varying complexity (Treguier 1999; Cerovecki et al. 2009).
As an example, the zonal average of δkgm between 100° and 60°W (where the ACC is roughly zonal and the zonal average is therefore along streamlines) is shown in Fig. 8a together with the zonal velocity and densities. This clearly shows that kgm is enhanced at depth levels where speeds of zonal velocities reside between 1 and 3 cm s−1, close to the observed phase speed of about 2 cm s−1 (Smith and Marshall 2009). The middepth enhancement of kgm is also estimated by Abernathey et al. (2010, their Fig. 4) in the form of an effective mixing keff, as defined by Nakamura (1996), although their critical layer depth is about 500 m shallower than found here in the present results. A pattern of enhanced kgm in regions represented by speeds of about 1–3 cm s−1 seems to be a general feature also in the Southern Ocean (not shown), regardless of some longitudinal variations there.
In areas with high velocities, on the other hand, the displacements caused by eddies cannot reach long distance (the length scale) before they are advected by the strong mean flow, resulting in weaker meridional mixing. This is referred to as strong current or PV gradient suppression (Ferrari and Nikurashin 2010) and usually occurs in the top layers of the main currents, where the velocities are high. This kind of suppression is quite evident in the thermocline levels in the ACC (Figs. 7a,c, 8a). The near surface suppression along the ACC is also estimated in the form of keff by Marshall et al. (2006) and Ferrari and Nikurashin (2010). A similar pattern is found in Gulf Stream (Fig. 8b), where a reduced kgm is in the thermocline levels where the PV gradient is large and enhanced kgm is around 1500 m. This is consistent with the observation of inhibited (intensified) exchange of floats in the upper ocean (at middepth) across the Gulf Stream jet (Bower et al. 1985; Bower and Rossby 1989; Bower and Lozier 1994; Lozier et al. 1997; Rogerson et al. 1999). We note that the close relation between the mean current velocity and the estimated kgm implies that the parameterization of an eddy mixing coefficient k is partly related to the difference between current velocity and wave phase speed (Marshall 1981; Killworth 1997; Abernathey et al. 2010). The sandwich structure of the estimated kgm in the ACC and the WBCs strengthened both the density fronts and the thermocline, which are too weak in the reference run compared to the observations. The changes in density field consequently improve the SSH, the velocity, and other fields. Additional to the improved fit of the model to the data, the estimated kgm contains reasonable physical structures, as discussed above.
2) Implications of the estimated kgm for possible parameterizations
A notable feature of δkgm is that it is quite small in the interior ocean away from the main currents, where the density slopes are weak (Fig. 7). Additional estimation experiments with background values 100 and 200 m2 s−1 for kgm and kredi (not shown) demonstrate that the regions with small δkgm remain the same, implying that the cost function is insensitive to kgm in these regions. We thus can assume that the background diffusivities in these regions are small; this could also be interpreted as that the processes due to baroclinic instability are very weak in low-density slopes. This interpretation agrees with Visbeck et al. (1997), who suggested a parameterization for kgm as
where the thermal wind relation has been applied to the second equation, α is a proportionality parameter denoting the mixing efficiency, N is the Brunt–Väisälä frequency, s is the density slope as denoted in the introduction, f is the Coriolis parameter, Ri is the Richardson number, and Lbc is the eddy mixing length scale. This parameterization shows that kgm is proportional to the density slope s.
Inspired by the above experiments and Visbeck et al.’s (1997) hypothesis, we assume arbitrarily that, in the insensitive regions where s is less than 2 × 10−4 and where |δkgm| is less than 50 m2 s−1, the background kgm is 50 m2 s−1, whereas in other regions the background kgm is kept as 800 m2 s−1 as before. We now obtain the estimated kgm by solving for δkgm according to the following specification:
The resulting kgm at 310 m is shown in Fig. 9b, which displays high similarity with the parameterized kgm by Eden et al. (2009, their Fig. 1d); additionally, the estimated kgm averaged between 50 and 1000 m (Fig. 9a) is also quite similar in both pattern and magnitude with the eddy transfer coefficient k estimated from observations by Stammer (1998, his Fig. 5c), both of whom applied the mixing length hypothesis, as done by Visbeck et al. (1997). Applying the mixing length hypothesis to our GFM data, we estimated the eddy diffusivity k according to
where Lbc and Tbc are estimates of eddy transfer length and time scales, respectively. Here, we follow Eden et al. (2007), who suggested that Lbc could be properly represented by the minimum of the estimated first baroclinic Rossby radius of deformation Lro and the Rhine’s scale Lrhi. The terms Lrhi and Lro can be estimated as and , where β is the meridional gradient of f, denotes the first baroclinic Rossby wave speed (Chelton et al. 1998), h is the local water depth, and Tbc is the inverse of the baroclinicity instability growth rate and can be estimated as but should be modified in weak stratification areas (Eden et al. 2009). This choice of the eddy length scale has shown better performance than others (Eden and Greatbatch 2008; Eden et al. 2009). Assuming α = 1, the estimated k with (3) at 310 m is shown in Fig. 9c, which resembles the large-scale pattern of kgm (Fig. 9b) well. Large differences exist in the high latitudes, where k is smaller than kgm, and in the north equatorial currents and the WBCs, where k is larger than kgm. The extremely large k in the tropical Pacific (Fig. 9c) might be due to the choice of the efficiency of eddy mixing parameter α, which may be too large here. Here, α is estimated to range from O(~0.01) (Visbeck et al. 1997; Stammer 1998) to 2 (Eden et al. 2009) but used as a constant. However, reasons for the differences in the high latitudes could be either from unrealistically estimated kgm due to other unresolved processes (e.g., ice–ocean interaction), from the incorrectly estimated eddy mixing length scale or a nonhomogeneous α in (3).
Using the mixing length hypothesis, we calculated k also directly from outputs of the eddy-resolving model STORM (https://verc.enes.org/community/projects-and-partnerships/projects/storm; J.-S. von Storch 2010, personal communication) as
where U is the characteristic eddy velocity, E is the eddy kinetic energy, and Lbc is the eddy mixing length scale defined as the same as in (3). Here, U and Lbc are averaged for 20 yr (1981–2000) and the estimated k according to (4) is shown in Fig. 9d. Again, similarities to kgm are apparent in pattern and differences rely on the magnitudes and the high latitudes. However, the similarity in pattern between the estimated kgm and the independent estimates of k demonstrates the skill for the estimation of kgm.
Given the good agreement between kgm and the estimates from the mixing length theory, the eddy mixing efficiency parameter α can be estimated according to
Its geographical structure (the dissimilarities from a uniform value) will provide further insight into modifications to the mixing length theory (Fig. 10). As an example, in 310 m and as a zonal average (Fig. 10a) the estimated α shows a primarily latitudinally dependent spatial pattern. High values in the high latitudes (the average in the Southern Ocean between 80° and 40°S is 1.9) are consistent with Eden et al. (2009), who used α = 2. Low values exist in the tropics (the average value between 20°S and 20°N is 0.5) (with regions of westward-flowing equatorial currents being exceptions because α here is higher than in the surroundings). Enhanced α occurs in the Southern Ocean between 1000 and 3000 m and in the Northern Hemisphere between 1000 and 2000 m: the critical layers (Fig. 10b). Larger α means more effective eddy transfer process subject to fixed time and length scales; on the other hand, because the estimates of Lbc and Tbc are also only approximations, a geographically varying α might suggest a demand of better estimation of the characteristic eddy mixing time and length scales.
b. The estimated kgmskew, kredi, and kz parameters
In this section, we will describe the spatial structure of the other three estimated mixing parameters. The physical interpretation of those changes is beyond the scope of this paper and will be left for future analysis.
Because the background value of kgmskew is 0, the adjustment of kgmskew is identical to its estimate. We note that, in contrast to kgm, which is anticipated to be positive to guarantee a down gradient mixing, there is no sign constraint on kgmskew (Eden et al. 2007) from physics. In the present work, positive (negative) kgmskew mean advection of mass properties along buoyancy contours with the higher density on the left (right) side.
A notable feature of kgmskew is that its magnitude O(±1000 m2 s−1) is of the same order as kgm (kgmskew should be 0 or at least very small for the down-gradient GM parameterization to be valid), although the largest values are mainly restricted to the Gulf Stream and Kuroshio regions (Fig. 11a). Large kgmskew also occur in some sections of the ACC, whereas moderate kgmskew occur in the equatorial currents systems. The distribution only slightly changes with depth, except that kgmskew is very small everywhere in the mixed layer (Fig. 11b). Two other physically interesting features are noteworthy: one is that large kgmskew changes occur in meandering regions of the currents and the spatial pattern of kgmskew is highly related to spatial variations of bottom topography; the other is that positive and negative kgmskew alternate along the current.
Taking the Southern Ocean as an example (Fig. 11b), positive and negative regions form regular cores with a radius of about 10°, correspondent with local bottom topographic scales. In the WBCs, the radiuses are different but still correspond to the local bottom topographic scales. We note that these structures are different from those obtained from high-resolution models by Eden (2006) and Eden et al. (2007), who found roughly latitudinal distribution of kgmskew, and also different from idealized models (Eden 2010), which suggested vanishing small kgmskew. The relation of the distribution of kgmskew and bottom topography suggests that a pure down-gradient hypothesis in GM may be modified by bottom topography at least in regions where PV is essentially influenced by the topography.
Figure 12 shows the distribution of δkredi. Generally, in the mixed layer and in the thermocline levels, regions of large δkredi resemble those of δkgm, except that δkredi is positive in the Gulf Stream and Kuroshio and is negative in the south flank of ACC. The kredi generally shows a (near) surface intensification and fast depth decay pattern in the Southern Ocean and in the WBCs. Negative kredi that exists in a few locations along the ACC is problematic from both physical and numerical perspectives, possibly compensating for too large numerical mixing or other model errors. However, in the high latitudes of the Northern Hemisphere, kredi is enhanced in the middle levels (between 1000 and 2000 m), which may exert significant influences on the deep water’s formation (Sijp et al. 2006; Sijp and England 2009).
In contrast to kgm, kgmskew, and kredi, which are adiabatic diffusivity coefficients, the vertical diffusivity coefficient mostly is a diapycnal diffusivity coefficient. The structure of δkz varies greatly in space. In the mixed layer (around 50 m; not shown), pronounced enhancements occur in the equatorial currents, in and around the South China, Philippines, and Caribbean Seas. The latter are all tropical regional seas and featured with a high complexity of topography at shallow depths, which enable strong mixing through the generation and breaking of internal tides. In the thermocline depths (Fig. 13a), positive adjustments generally occur in the western part of all tropical oceans and negative adjustments generally occur in the eastern Pacific tropics (especially true below 350 m). Below the thermocline, however, there are nearly no adjustments, which is mainly due to the short estimation period. This consequently exerts little effect on the change of MOC as expected (e.g., Bryan 1987; Mignot et al. 2006). The zonally averaged δkz (Fig. 13b) shows a representative north–south equatorial symmetry between 20°S and 20°N. It also shows two apparent cores of positive adjustments in depth in the tropical oceans, with one centered slightly above 100 m and the second centered slightly below 200 m, roughly at the location of the subsurface salinity maxima. In between the two positive cores is a thin layer of weak negative adjustments. The vertically varying structure in the tropics substantially influences the thermocline. Separating the last term of the tracer equation (1) into two parts
an evaluation for the GFM yields that the first part on the rhs, which is due to the vertical variation of kz, dominates over the second part. However, the adjustments of kz are small compared to its background, which are likely due to the small a priori error uncertainties.
5. Conclusions and discussion
Using oceanic state estimation technology, it is possible now to estimate not only the atmospheric forcing fields but also the internal model parameters so that the modeled oceanic state fields are consistent with global in situ and satellite observations. In the present study, we estimated four isopycnic/diapycnic mixing parameters jointly with initial conditions and six surface forcing fields over a period of 10 yr. After optimization, the model has moved closer to the observations that are imposed as constraints. The optimized tracer mixing coefficients reduce about 10% of the model–data misfit and changes occur mainly in the thermocline, where these coefficients take effect, and for the SSH, which depends on the distribution of density. However, misfits remain (Fig. 3b), especially in the main currents, in the subpolar North Atlantic, which cannot be further reduced by adjusting the present 12 control parameters because of insensitivity of the former to the latter. Instead, the misfits might be caused by the model resolution, the boundary conditions, or other deep ocean processes. To further improve the estimate, other factors, such as the inconsistency between different data, more accurate data, a priori errors, and error covariance, should be taken into consideration.
Moreover, the adjusted four eddy tracer mixing parameters reflect valuable information about the physical processes in the real ocean. Among the other diffusivities, we found that the thickness diffusion parameter kgm plays the most important role in redistributing water mass properties. The estimated kgm shows a vertical sandwich structure over the main currents; for example, it is reduced in the main thermocline layer while enhanced in the surface and in the middepth critical layers. This feature corresponds to the PV gradient suppression and critical layer enhancement hypotheses very well. In addition, the estimated kgm also agrees well in pattern with independent estimates of eddy transfer coefficients obtained from mixing length theory. The consistencies summarized above suggest that our estimated kgm shares many properties of the real mixing processes.
Another finding is the alternating patterns of kgmskew along the main currents. Acting as streamfunction of eddy-induced advection (Eden et al. 2007), the estimated kgmskew mimics cyclonic and anticyclonic modifications to the main currents. In addition, the structure of kgmskew is highly related to the bottom topography and may reflect some poorly resolved interactions with topography.
Finally, the estimation period is of concern if one wants to estimate parameters that are related to long time-scale processes in the oceans, which are usually associated with mixing or advection on climate time scale. For example, 10 yr is not enough for kz to fully take effect, except for the top 500 m, because it requires about O(102) yr for vertical mixing to diffuse the properties in the deep ocean. This underlines the fact that long dynamically consistent estimation efforts are required to improve our insight into ocean processes. We are especially awaiting longer runs to learn more about adjustments in the vertical mixing to lead to improved model simulations.
We are very grateful to the two anonymous reviewers for their very constructive comments, which greatly improved the manuscript. The integrations of the model simulations were performed at the German high-performance computing center DKRZ with support from University of Hamburg. The STORM results were kindly provided by Jin Song von Storch. Supported in part through the Exzellenzcluster, Integrated Climate System Analysis and Prediction (CliSAP) funded through the Deutsche Forschungsgemeinschaft.