Abstract

An estimate of the time-varying global ocean circulation for the period 1992–2002 was obtained by combining most of the World Ocean Circulation Experiment (WOCE) ocean datasets with a general circulation model on a 1° horizontal grid. The estimate exactly satisfies the model equations without artificial sources or sinks of momentum, heat, and freshwater. To bring the model into agreement with observations, its initial temperature and salinity conditions were permitted to change, as were the time-dependent surface fluxes of momentum, heat, and freshwater. The estimation of these “control variables” is largely consistent with accepted uncertainties in the hydrographic climatology and meteorological analyses. The estimated time-mean horizontal transports of volume, heat, and freshwater, which were largely underestimated in the previous 2° optimization performed by Stammer et al., have converged with time-independent estimates from box inversions over most parts of the World Ocean. Trends in the model’s heat content are 7% larger than those reported by Levitus and correspond to a global net heat uptake of about 1.1 W m−2 over the model domain. The associated model trend in sea surface height over the estimation period resembles the observations from Ocean Topography Experiment (TOPEX)/Poseidon over most of the global ocean. Sea surface height changes in the model are primarily steric but show contributions from mass redistributions from the subpolar North Atlantic Ocean and the Southern Ocean to the subtropical Pacific Ocean gyres. Steric contributions are primarily temperature based but are partly compensated by salt variation. However, the North Atlantic and the Southern Ocean reveal a clear contribution of salt to large-scale sea level variations.

1. Introduction

In this paper we describe the results of a global synthesis obtained by the Estimating the Circulation and Climate of the Ocean (ECCO) Consortium (see Stammer et al. 2002a) on a 1° spatial grid (±80°) for the 11-yr period 1992–2002. The synthesis is obtained by forcing the ECCO ocean circulation model to near consistency (within specified error margins) with most of the World Ocean Circulation Experiment (WOCE) observations by using the model’s adjoint to modify the initial temperature and salinity conditions over the full water column and to adjust the time-varying meteorological forcing fields over the full estimation period on a daily basis. This work is an extension of a previous paper (Stammer et al. 2002b, hereinafter referred to as SEA02) that described a similar synthesis but one that was obtained on a coarser 2° grid and over a shorter period while using significantly less data and less accurate model physics. As before, our approach is to use a general assimilation method that produces dynamically self-consistent results. Stammer et al. (2004) discussed the skill of the surface forcing fields estimated as part of the same solution (see below).

We embark here on an evaluation and discussion of the resulting estimates of the 11-yr mean and interannual to decadal changes in the ocean state and property transports. Ocean variabilities on those time scales gained substantial attention during recent years and are a central focus of the World Climate Research Programme’s (WCRP) Climate Variability and Predictability Program (CLIVAR; information online at http://www.clivar.org). Based on observations, Levitus et al. (2001) reported changes in temperature and salinity of hitherto unknown amplitudes not only near the surface but also in the actual abyssal oceans. Cazenave and Nerem (2004) pointed out that the processes contributing regionally and globally to sea level changes are numerous and that our understanding of them is far from being complete; Antonov et al. (2002) and Munk (2003) discussed the importance of freshwater input, and Miller and Douglas (2004) suggest that the increase in mass (i.e., freshwater) is actually more important for sea level rise than the effect of warming.

Recently, Bryden et al. (2005) suggested from an analysis of a repeated hydrographic section along 24.5°N in the North Atlantic Ocean that the Atlantic meridional overturning circulation (MOC) has slowed down by about 30% between 1957 and 2004. Whereas in their analysis the northward transport of the Gulf Stream across 25°N remained nearly constant, the decrease of the MOC is evident both by a 50% larger southward-moving midocean recirculation of thermocline waters, and also by a 50% decrease in the southward transport of lower North Atlantic Deep Water between 3000 and 5000 m in depth. Discussions in the wake of this publication clearly demonstrate the controversy that can result from an isolated analysis of integral quantities. For a better understanding of the changes of the ocean’s transports, a wider analysis of changes in the ocean’s local and global mass, heat, and salt content, as well as of the related sea level changes, need to be considered.

In principle, an ocean synthesis provides the framework for such a combined estimate of the time-varying circulation of the ocean that is consistent with prior model and data error descriptions. The global ECCO ocean-state estimation effort supplies an estimate of the time-varying ocean circulation as the basis for quantitative studies of climate dynamics. The solution combined many ocean observations with a numerical model to produce a dynamical description of the changing ocean, which serves as a basis for understanding future changes in ocean circulation. It is expected to serve also as a basis for computations of not directly observable changes of integral ocean quantities that are consistent with the model physics, all available data, and the atmospheric forcing fields. Results shown below will be a first attempt to analyze regional and global changes in the ocean circulation as it results from the available ECCO 11-yr-long global ocean-state estimate.

In the following, we will provide in section 2 details of the model and the optimization. In section 3 the remaining model–data misfit is discussed and section 4 describes the basic features of the resulting model state. The overturning streamfunction and meridional heat and freshwater transports will be described in section 5. Section 6 then will deal with changes in heat content, while section 7 will discuss variations in sea surface height (SSH) on decadal time scales and their relation to steric and mass variations. A concluding discussion will follow in section 8.

2. Method

First, we briefly summarize the model setup and the optimization, which builds on the ECCO forward and adjoint model. SEA02 provide further details of the general assimilation approach and the underlying infrastructure.

a. The model

The ECCO model is based on the Massachusetts Institute of Technology (MIT) general circulation model (Marshall et al. 1997a, b), which consists of conservation equations for horizontal and vertical momentum, volume, heat, and salt, as well as an equation of state, which are solved on a staggered C grid (Arakawa and Lamb 1977). Spatial coordinates are longitude, latitude, and height. A detailed description of the model is provided by Adcroft et al. (2002). For our application, the ECCO model was configured on a 1° horizontal grid over ±80° in latitude and on a vertical grid with 23 levels (as in SEA02 with one additional 500-m-thick level).

For the present purposes, we use a hydrostatic version with an implicit free surface. A Gent–McWilliams (GM) eddy parameterization scheme (Gent and McWilliams 1990) was implemented and a full K profile parameterization (KPP) surface mixed layer model (Large et al. 1994) is used. Background coefficients of vertical diffusion and viscosity are 10−5 and 10−3 m2 s−1, respectively. The GM advection coefficient is 103 m2 s−1. Isopycnal diffusivity and horizontal viscosity were chosen as 102 and 104 m2 s−1, respectively.

The optimization is initialized with the Levitus et al. (1994a, b) January potential temperature and salinity fields after the baroclinic flow field was spun up during a 1-month forward integration. Daily surface heat and freshwater forcing fields, obtained from the first National Centers for Environmental Prediction–National Center for Atmospheric Research (NCEP–NCAR) reanalysis project (Kalnay et al. 1996), and twice-daily wind stress fields are used to force the model at its surface. The meridional boundaries are closed; that is, no volume flux occurs normal to them. In addition, no relaxation of the model fields toward the climatological temperature T and salinity S fields was applied at the poleward lateral boundaries to avoid artificial sinks and sources of heat and freshwater. Moreover, no relaxation was applied in the Strait of Gibraltar since the Mediterranean Sea is part of the model domain. As a result, in our runs all sources for heat and freshwater reside in the surface boundary only.

The model topography was derived from the Smith and Sandwell (1997) topography provided on a 1/30° horizontal resolution and 5′ Gridded Earth Topography (ETOPO5) poleward of 72.006° latitude. In each 1° box, the median of the ETOPO5 data was calculated and the resulting topography was subsequently edited manually to improve the representation of coastline, islands, straits, and ridges. Special attention was paid to the topography in critical areas such as (a) the Caribbean Sea and the Straits of Florida, (b) the Denmark–Scotland overflow, (c) the Indonesian Throughflow, and (d) Drake Passage. The tuning of the topography resulted in mean transports through these key regions as summarized in Table 1.

Table 1.

Volume and enthalpy [with reference to 0°C to enable comparison with Ganachaud and Wunsch (2000)] transports through channels and passages.

Volume and enthalpy [with reference to 0°C to enable comparison with Ganachaud and Wunsch (2000)] transports through channels and passages.
Volume and enthalpy [with reference to 0°C to enable comparison with Ganachaud and Wunsch (2000)] transports through channels and passages.

b. The optimization

A schematic of the resulting optimization setup is given in Fig. 1. Further details are described in the appendix. The top part of Fig. 1 shows the data constraints and the periods for which the data are available. For a detailed summary of all datasets that are being used in the estimation and their prior error statistics, see Lu et al. (2002). They include several satellite datasets, surface drifter velocities, and in situ hydrographic temperature and salinity profiles, as well as WOCE hydrographic sections (see caption for details). We also constrain the model’s monthly mean climatology of temperature and salinity to the Levitus et al. (1994a, b) climatology. Surface forcing fields are constrained through daily NCEP surface fluxes and scatterometer wind stress. In contrast to our earlier estimation experiments (e.g., SEA02; Köhl et al. 2002) we use the difference of the time-mean Ocean Topography Experiment (TOPEX)/Poseidon (T/P) SSH field minus the modern Gravity Recovery and Climate Experiment (GRACE) geoid model (GGM01S; Tapley et al. 2003) to constrain the time-mean dynamic surface topography. In setting up the cost function, J, the mean and time-dependent components of surface elevation were separated, thus isolating errors owing to the geoid from the distinctly different ones in the time-evolving components.

Fig. 1.

Schematic of the optimization. (top) The data constraints imposed upon the model. The lines indicate times when data are available; mean or climatological data are shown as available throughout the whole period. (bottom) The control parameters that were changed during the optimization. These include the adjustments to the initial potential temperature (θ) and salinity (S) fields, as well as 2-day averages of the surface forcing fields over the full 11 yr.

Fig. 1.

Schematic of the optimization. (top) The data constraints imposed upon the model. The lines indicate times when data are available; mean or climatological data are shown as available throughout the whole period. (bottom) The control parameters that were changed during the optimization. These include the adjustments to the initial potential temperature (θ) and salinity (S) fields, as well as 2-day averages of the surface forcing fields over the full 11 yr.

The lower part of Fig. 1 shows the control variables that are being adjusted during the optimization so as to bring the model into consistency with the data. In the present calculations those control parameters include the adjustments to the initial potential temperature (θ) and S fields, as well as 2-day averages of surface forcing fields over the full 11 yr. A combination of variations in these controls is determined that leads to the best fit in a least squares sense of the model state with respect to observations and their uncertainties over the full data period. Surface forcing fields were split into time-mean and time-varying components. We assume here that the model uncertainties reside primarily in the initial conditions and the surface forcing fields; accordingly, the model simulation was brought into consistency with the observations through changes, within error bars, in both the initial conditions and the NCEP–NCAR reanalysis surface forcing. A justification for this assumption was provided by SEA02 and Stammer et al. (2007), based on the fact that these parameters are known to have large uncertainties and that under similar conditions previous estimation attempts obtained a solution that was consistent with prior model–data uncertainties. A major challenge for obtaining the best possible estimates of interannual and longer-term variability of the ocean in ocean-state estimates separating the numerical model drift from the time-varying processes in the ocean and minimizing the former without affecting the latter. Moreover, model drift is related to uncertainties in the surface forcing that we intend to estimate along with the estimate of the ocean circulation and is therefore expected to decrease toward the minimum of the cost function.

To help in this process, estimated changes in the initial temperature and salinity fields were required to stay within the uncertainty limits of the Levitus hydrography; that is, the model initial T, S conditions were required to stay—within uncertainty limits—close to the Levitus T, S January climatology. In addition, the monthly T and S climatology from the 11-yr period (as opposed to individual monthly mean fields) were constrained by the monthly mean Levitus et al. (1994a, b) climatology, thereby constraining water masses while at the same time allowing interannual variability to occur. Last, we constrained the pointwise model drift in temperature and salinity over the 11-yr period by minimizing a term that measures the rms difference between the last-year mean minus the first-year mean. The function of this latter term is to suppress artificial (numerical) model drift, arising, for example, from global imbalances in the NCEP–NCAR forcing, in regions with insufficient in situ data.

Without any extra constraint, changes in T and S would be independent and not bound by any dynamical constraint. Such changes, when they are not in geostrophic balance, are known to lead to large inertial oscillations [the initialization problem; see Daley (1991)]. To control the geostrophic imbalance, which will locally cause large horizontal divergences and convergences and thus vertical velocities, an explicit constraint on the difference between the first- and the last-year time-mean amplitudes of the vertical velocity was introduced. The rms difference was required to stay within the rms limits of 2 × 10−6 m s−1 (the corresponding value of “natural” fluctuations of an unperturbed simulation). The effect of this constraint is that initial vertical velocities normal to the topography are characterized by an enhanced amplitude of the bottom torque term,

 
formula

and thus an anomalous large initial vertical circulation is minimized.

We note that observed drifts of the sea surface height, and in temperature and salinity fields over the full model depth, are not eliminated by the constraints on the model drift. The effect of the constraints will be discussed in more detail below.

3. Model–data misfit

The optimization started from the earlier solution provided by (SEA02) on a 2° spatial grid by taking into account its estimated time-mean adjustments of the NCEP–NCAR surface forcing fields. A total of 48 iterations were required subsequently for the solution to converge to the point that further reductions in the cost function between successive iterations became insignificant. In this context it is important to recall that one cannot expect a solution like this to fully converge (within error bars) given this number of iterations. Without costly preconditioning, such a convergence can only be expected after as many iterations as the ratio of the rank of the system to the rank of the noise. With about 108 formally independent control parameters, there is no reason to believe that such a system is only of the order of rank 100 and would be converged after the given number of iterations. However, we use a simple diagonal preconditioning and we do expect diminishing returns as the eigenvalues of the remaining errors become whiter in the spectrum, and therefore one would also start to fit the model to the observational noise. Given our limited computer resources, we decided (somewhat subjectively and based on the slow additional convergence and reduced adjoint gradients) to stop the optimization at this point and start scrutinizing it for later improvements.

As described in the appendix, the total cost function consists of contributions obtained from model–data differences and from the control terms. During the course of the optimization, the data-related part of the total cost function decreased to about 50% of its initial value. Before we turn to an analysis of interannual and longer-period fluctuations in the resulting solution, we will first provide a discussion of the remaining model–data residuals, the estimated ocean state itself, and its time-mean transports.

To begin with, Fig. 2 shows the contributions of individual datasets to the cost function. The top panel in Fig. 2 shows the logarithm of all terms that contribute to the total model–data misfit of the optimized state. Color coding of the results, discriminating between the dominant, important, and less influential datasets, demonstrates that the remaining total misfit is mostly determined by only a few variables, notably the daily time-varying SSH fields, the monthly mean Reynolds SST fields, and the monthly mean Levitus temperature and salinity climatology. A reduced contribution results from the Tropical Rainfall Measuring Mission (TRMM) Microwave Imager (TMI) SST dataset, the time-varying surface forcing, and the XBT temperature profiles. Argo temperature profiles and the model’s drift in temperature, salinity, and vertical velocity produce an even smaller contribution. All other components are orders of magnitude smaller. However, this does not necessarily imply insignificance for the final state. As an example, changing only the mean SSH data in the optimization [from a T/P–Earth Gravity Model (EGM96) to a T/P–GRACE field] reveals that the mean SSH field projects strongly on the corrections to the initial conditions even though its relative contribution to the total cost function is small (see Stammer et al. 2007 for details).

Fig. 2.

Cost function contributions from individual datasets for the optimized state. (top) The bars represent the logarithm of the total model–data misfit and each individual contribution. (bottom) The square root of those values after normalizing each column by the number of the respective observations. The labels T and S describe the Levitus terms followed by the initial conditions TO and SO. Here, τx and τ y and Hq and Hs correspond to the zonal and meridional wind stress and the heat flux and freshwater flux corrections (the subscript m indicates the mean). In addition, drift(T), drift(S), and drift(w) are the drift terms, scatx and scaty correspond to the zonal and meridional scatterometer wind stress data, and “velocity” corresponds to the mean velocity field from the drifter data. Other labels stand for the data type.

Fig. 2.

Cost function contributions from individual datasets for the optimized state. (top) The bars represent the logarithm of the total model–data misfit and each individual contribution. (bottom) The square root of those values after normalizing each column by the number of the respective observations. The labels T and S describe the Levitus terms followed by the initial conditions TO and SO. Here, τx and τ y and Hq and Hs correspond to the zonal and meridional wind stress and the heat flux and freshwater flux corrections (the subscript m indicates the mean). In addition, drift(T), drift(S), and drift(w) are the drift terms, scatx and scaty correspond to the zonal and meridional scatterometer wind stress data, and “velocity” corresponds to the mean velocity field from the drifter data. Other labels stand for the data type.

The remaining model–data misfits can provide information about model errors, about remaining data errors, or inconsistencies in the prior data error statistics. This will be analyzed in the following to assess the success of the optimization. To that end, we show in the bottom panel of Fig. 2 the square root of the individual cost function terms after normalizing them by the number of the corresponding observations. Although formally the approximate character of the errors will always prevent the result from passing the χ2 test (Bennet 2002), a value much different from one indicates that either the model or the estimation of the error is largely incorrect and cannot be adjusted without changing the result significantly. Any value larger than unity corresponds to the factor by which the rms model–data difference exceeds the a priori error information. Relative values smaller than 1 are obtained for all forcing fields and the time-varying SSH anomalies. A value of about 1 is obtained only for the mean velocities from drifters. Contributions between 1 and 2 can be considered to be approximately consistent with the prior error. The rms values of all other data reside above 2 and indicate significant inconsistencies between the existing solution and those data or their prior error information. We recall, however, that the latter terms had a small contribution to the total misfit and possibly to the solution itself. In other words, the solution did not converge to the respective constraint and the residuals are larger than the a priori error bars.

The largest columns in Fig. 2b result from the model’s drift in temperature, salinity, and w, which indicates either an underestimation of the associated errors or too large a model drift. Our prior error values for the temperature drift, applied in the cost function, span 0.1°–0.33°C (over the 11-yr period) in the upper 700 m. For a comparison, Levitus et al. (2005) recently provided amplitudes of the observed horizontally averaged rms differences of the temperature anomalies from 1992 to 2002 that were between 0.25° and 0.85°C for the same depth range. This indicates that our prior statistics for temperature and salinity drifts in the ocean were possibly a factor of 2 too small. After adjusting for those errors, the final normalized temperature and salinity drift terms would be less than 2 and thus closer to a statistically consistent solution.

The question remains as to how critical the underestimation of the drift error is for our solution; that is, to what extent do we overconstrain the long-term changes in our estimate that actually take place in the ocean? In this context we recall that the estimate is only sensitive to the details of the drift term if this term is as large as the part of all the observational data that also controls directly or indirectly the trends in the solution. Figure 2a shows that the contribution from all surface and subsurface temperature observations [i.e., contributions SST + TMI + CTD(T) + XBT + ARGO(T)] exceeds by far the term originating from the model temperature drift [drift(T)] despite the overestimation of the drift weight. With a ratio of 1, of the relative contributions from the drift term to the temperature data, we would anticipate a significant impact from the overconstrained drift term on the model’s temperature tendency since most of the temperature data also include information about the temperature trend. Instead, the ratio is only 0.14 in our solution (it would be close to 0.034 with a more correct, smaller weight for the drift constraint, which would reduce the present value to only ¼ of it). Moreover, there are many data sources that indirectly impose information about long-term changes, such as the SSH data. The ratio between the drift constraint term and all terms resulting from time-dependent data, which might indirectly constrain the trend, decreases to 0.02 (this also holds for salinity). Our conclusion is therefore that the temperature trend in our solution is primarily imposed by the temperature and probably also other time-dependent data.

Analyzing the salinity drift [drift(S)], the condition reverses in that many fewer salinity observations are available that can directly impact the estimated salinity tendency. The ratio of the salt drift term to the sum of all time-dependent direct salinity cost function contributions [the sum over SSS + CTD(S) + ARGO(S)] is 4.7, suggesting that the solution is strongly determined in its salinity drift by the salt drift constraint and not by the direct observations (indirect observations such as SSH still contribute). However, a predominance of the estimated salt trend through the drift constraint is still unlikely because in this case we would expect a significantly reduced drift(S) term in comparison with the drift(T) term. Instead, both are of similar magnitude (cf. Fig. 2b), leading us to the hypothesis that the estimation of salinity is mostly determined by indirect observations such as SSH plus temperature data.

A closer inspection of the spatial distribution of the misfits between the model and XBT profiles points to local model problems or alternatively to local inconsistencies with the a priori representation error. As an example, Fig. 3a shows the horizontally averaged rms difference of our optimized solution with the XBT data as a function of depth and time. Prior to the averaging, values were normalized by the inverse of the local a priori error so that, again, a value around 1 would indicate the consistency of the solution. As can be seen from the figure, relative rms values of 2 or larger are mostly restricted to regions below 500 m and are mainly associated with the subpolar North Atlantic (Fig. 3b). A comparison of the model–XBT misfits obtained from the unconstrained model (Figs. 3c and 3d) shows that the assimilation is able to correct for a substantial fraction of the initial model–XBT discrepancies in most regions by changing the initial conditions or surface forcing fields. However, the enhanced (relative to other geographical locations) total model–XBT misfit remains in the North Atlantic, which presumably results from erroneous northern boundary conditions: as will be discussed below, problems in the vertical mixing/convection, and especially those associated with the Denmark Straight overflow, cause the model to deviate from the observed temperature distributions in those regions. A similar overflow problem is also present west of Gibraltar and is associated with Mediterranean outflow. In both cases, the poor representation of the overflow physics results in too shallow a depth of the signature of the overflowing water.

Fig. 3.

Horizontal and vertical distribution of (top) the rms difference between the optimized and XBT data and (bottom) the reference run and XBT data. The rms values are shown after multiplication with the local weight.

Fig. 3.

Horizontal and vertical distribution of (top) the rms difference between the optimized and XBT data and (bottom) the reference run and XBT data. The rms values are shown after multiplication with the local weight.

Figure 4 provides information about the model bias in the form of time series of the horizontally averaged model–data differences and maps of the vertically averaged time-mean differences. These results are shown separately for the optimization run and the first-guess solution. Figure 4 reveals too cold a temperature field in the optimized run above 100-m depth and too warm an intermediate range between 100 and 300 m that shifts to 300 to 500 m after 1997. In the deep ocean the bias appears to be much smaller than the rms differences. The horizontal distribution of the vertically averaged bias shows a positive bias only in the subpolar regions and the Indian Ocean. The drift in the unconstrained run is much larger near the surface but is similar below 200 m although in the optimized run the rms differences are also significantly smaller below 200 m (see Fig. 3). The largest temperature bias can be found in the subpolar North Atlantic and near boundary currents.

Fig. 4.

Horizontal and vertical distribution of the vertically and horizontally averaged difference between (top) the optimized run and XBT data and (bottom) the reference run and XBT data. The mean values are shown after multiplication with the local weight.

Fig. 4.

Horizontal and vertical distribution of the vertically and horizontally averaged difference between (top) the optimized run and XBT data and (bottom) the reference run and XBT data. The mean values are shown after multiplication with the local weight.

It appears that one of the main reasons for the remaining model–data inconsistencies with many datasets is caused by incorrect water mass formation processes in the model’s subpolar North Atlantic. This points to the necessity of including overflow parameterizations and mixed layer parameters as controls into the optimization. Despite those remaining inconsistencies, Fig. 3 demonstrates the improved skill of the optimized state in simulating in situ observations as compared with the control run and suggests that our results are better suited for an analysis of the changing ocean circulation over large parts of the global ocean.

Some remaining isolated, larger misfit values exist, which upon inspection usually emerge as the result of bad observations that should have been removed. Testing the data quality through model–data comparisons may be a valuable feedback for CLIVAR and its ongoing reanalysis efforts. Since we assume uncorrelated errors and since the forcings are mapped fields, the associated degrees of freedom are biased toward those of the forcings that limit the size of the corrections. This problem will be remedied in the future with the use of error covariance matrices, which would better reflect the actual degrees of freedom in the forcing and the model–data differences and would give more weight to the model–data misfit terms rather than the mean forcing/control terms.

4. The estimated ocean state

Figure 5 (top) shows the estimated mean sea surface height field from the 11-yr period. In comparison with the earlier solution on a 2° grid (Fig. 7 in SEA02), the higher resolution of the present model improved the structures in the simulated circulation. Nevertheless, the model still lacks mesoscale variability as well as the observed sharpness of the boundary currents and, therefore, cannot be in full agreement with the data. As compared with the 1° unconstrained model simulation (first guess), the parameter adjustments lead in general to stronger gyres by [O(10 Sv); Sv ≡ 106 m3 s−1]. Equally important, the optimization leads to a more realistic separation of the Gulf Stream by shifting its axis farther south and by allowing cold fresh slope water to flow southward along the shelf slope between Flemish Cap and Cape Hatteras. It does this at the expense of large forcing adjustments. The pattern and size (maximum values are about a factor of 5 larger than the prescribed error) of those forcing corrections indicate that the wind stress corrections, in particular, compensate for model deficiencies. Caution is therefore advised if these corrections are to be used in model configurations that do not show these deficiencies. However, Menemenlis et al. (2005) show that applying the corrections to a ¼° resolution model leads also to a better simulation of the observed data. Moreover, Stammer et al. (2004) demonstrated by comparing with satellite-based wind stress estimates the skill of the wind stress corrections in regions away from the energetic boundary currents, which allows us to distinguish between the compensation for model deficiencies and realistic corrections.

Fig. 5.

(top) Time-mean SSH [contour interval (CI) = 10 cm] and (bottom) rms SSH variability (CI = 2 cm). A linear trend (shown in Fig. 13) was removed from the daily SSH fields before calculating the rms variability.

Fig. 5.

(top) Time-mean SSH [contour interval (CI) = 10 cm] and (bottom) rms SSH variability (CI = 2 cm). A linear trend (shown in Fig. 13) was removed from the daily SSH fields before calculating the rms variability.

Figure 5 (bottom) shows the rms model SSH variability. Much of the variability in the boundary currents and their extensions is associated with the seasonal warming–cooling cycle. In contrast, the variability in the Tropics originates from the changes of the flow field on seasonal time scales and from Rossby wave activity there. However, most of the enhanced variability on the western and eastern sides of the Pacific Ocean and in the Indian Ocean are associated with the 1997/98 ENSO event. In the Southern Ocean, enhanced variability southwest of Australia and upstream of Drake Passage is associated with high-latitude barotropic activity.

Figure 6 shows the mean difference of the estimated SSH minus T/P–GRACE observations. Residuals are of the order of ±10 cm and reach ±40 cm on relatively short spatial scales near steep topography and over the Southern Ocean. An inconsistency with the prescribed error of the mean SSH of 4.5 cm is obvious (we recall that the normalized rms SSH misfit is around 1.5 in Fig. 2b; i.e., slightly larger than the expected value of 1), especially in those latter regions. It is unclear whether the estimation from the model, or from the TP data, or from the geoid contains the largest error. Significant errors most likely reside in all three components, as discussed in detail by Stammer et al. (2007).

Fig. 6.

Difference between (top) the estimated mean SSH and the T/P–GRACE mean SSH observations (CI = 10 cm) and (bottom) rms difference of the daily SSH fields with respect to T/P observations (CI = 2 cm).

Fig. 6.

Difference between (top) the estimated mean SSH and the T/P–GRACE mean SSH observations (CI = 10 cm) and (bottom) rms difference of the daily SSH fields with respect to T/P observations (CI = 2 cm).

Figure 6 (bottom) shows the rms model–T/P difference. Areas with the largest deviation from the observed SSH variations coincide with the most energetic regions of the ocean because of the missing mesoscale dynamics of the model. This is in accord with the prior assumption that 50% of the T/P SSH variance cannot be resolved and has to be attributed to the model representation error (note that the respective misfit term is less than unity in Fig. 2b). The low model variability apparent in the eastern north equatorial region between 5° and 10°N is due to the coarse resolution used, which leads to an unrealistically large damping of annual Rossby waves originating from the eastern side of the basin and of tropical instability waves. It is also noteworthy that upstream of Drake Passage the pattern of the rms misfit shows a striking resemblance to that of the high-energy barotropic region there, which points at a possibly large aliasing error in satellite altimetric observations (cf. Stammer et al. 2000).

Examples of the time-mean flow field are provided in Fig. 7, which shows vertical sections of the zonal and meridional velocity from the Atlantic, Pacific, and Indian Oceans. Those sections were chosen to represent each basin by crossing major current systems. Because of the seasonal reversal of the monsoon winds, we show only the August climatology, which represents the summer monsoon state. As is to be expected, the circulation reveals more details than were available from the previous 2° solution. In particular, the stack of several zonal currents and countercurrent bands of smaller meridional and vertical extent in the tropical Atlantic and Pacific (Figs. 7a and 7c) is now better resolved than in the 2° solution. We note that in Fig. 7a the zonal section along 30°W shows a pronounced Azores Current approximately at the observed position at 35°N (Käse and Krauss 1996) whereas the control run (not shown) reveals only a weak current farther south.

Fig. 7.

(left) Meridional and (right) zonal sections of the time-mean zonal and meridional velocities, respectively. Sections are shown for the (top) Atlantic, (middle) Pacific, and (bottom) Indian Oceans. Note the nonlinear color scale.

Fig. 7.

(left) Meridional and (right) zonal sections of the time-mean zonal and meridional velocities, respectively. Sections are shown for the (top) Atlantic, (middle) Pacific, and (bottom) Indian Oceans. Note the nonlinear color scale.

The Deep Western Boundary Current (DWBC) in the North Atlantic is only represented as a single core at about 2500 m in Fig. 7b, and no separation in either of the two separate cores exists as they formed by the convection and the overflow water type branches of the North Atlantic Deep Water (NADW), respectively. In the Pacific the flow field across 20°N appears to be surface intensified (Fig. 7d) with only weak velocity amplitudes present at intermediate-depth levels. However, near steep bottom topography, we encounter enhanced bottom-intensified amplitudes. Most noticeable for the Indian Ocean (Fig. 7f) are the strong subsurface currents associated with the seasonal overturning. The Pacific shows a deep western boundary current (Fig. 7d); however, it is weak and at about 1200-m depth [similar to the simulation by Maltrud et al. (1998)] shallower than the observed (below 2000 m).

5. Overturning, heat, and freshwater transports

Because the meridional mass transport streamfunction and the meridional transports of heat and freshwater are among the climatically important elements of the ocean circulation, much effort was spent in the past on determining their time-mean amplitudes, especially during the World Ocean Circulation Experiment (WOCE). The mean meridional mass transport streamfunction,

 
formula

is plotted in Fig. 8 (top) for the global ocean. The maximum downward transport is close to 18 Sv in the Northern Hemisphere and is associated mostly with the Atlantic Ocean. The deep inflow of Antarctic bottom water amounts to around 10 Sv and can be found mostly in the Pacific and Indian Ocean sector of the Antarctic Circumpolar Current (ACC).

Fig. 8.

Time-mean meridional overturning streamfunctions (Sv) for the global ocean calculated on (top) isobaric and (bottom) isopycnic surfaces (CI = 5 Sv).

Fig. 8.

Time-mean meridional overturning streamfunctions (Sv) for the global ocean calculated on (top) isobaric and (bottom) isopycnic surfaces (CI = 5 Sv).

Relative to the previous 2° solution (SEA02), the largest differences in the MOC estimate are found in middepths levels where the 1° solution reaches higher values and has more meridionally coherent structures. Increased model resolution and the inclusion of the Gent and McWilliams (1990) parameterization are the prime candidates for the improvement. Figure 8 (bottom) shows the mean overturning streamfunction, but computed along density (σ) surfaces according to

 
formula

It follows from the figure that most of the cross-isopycnic transport occurs as part of the near-surface Ekman cell. Across the ACC we can find an additional loop transporting surface water northward across isopycnals from 28 to 26 σt. Basin-wide exchanges of up to 15 Sv are limited to densities above 26 σt. This loop reaches across the entire Atlantic basin from the ACC to 60°N, transforming water from about 26.5 σt to 27.5 σt; the southward returning water leaves the Atlantic at a density of close to 27.8 σt.

Figure 9 shows transport time series of the Gulf Stream and the Deep Western Boundary Current, both evaluated at 43°N. Note that the model’s Gulf Stream transport (the sum of all northward transports in the core cell at 43°N) declines after 1994. Also shown in Fig. 9 is the strength of the meridional overturning circulation at 900-m depth along the same latitude. Both transport time series agree fairly well with the temporal changes observed in the MOC; we note that the transports are actually leading the MOC variations by a few months. An explanation for the decline of the western boundary current transports is provided below.

Fig. 9.

Anomalies of meridional cell transports. Only velocities in the direction of the mean current of the cell are summed up. (top) The transport anomalies are shown as they result at 43°N, by integrating the deep southward branch (black, mean value of 22.5 Sv) and northward Gulf Stream branch (red, mean value of 28.5 Sv) between 50° and 35°W. Also shown is the amplitude of the MOC at the same latitude at 900-m depth (green, mean value 17.5 Sv). (bottom) The southward transports are shown from the same latitude after integrating the southward velocities between 35° and 10°W above 1000-m depth (black, mean value of 9.5 Sv) and below 1000-m depth (red, mean value of 10.0 Sv). Also shown are the anomalies of the barotropic streamfunction computed at 35°W (green, mean value of 7.0 Sv) and the NAO index (blue).

Fig. 9.

Anomalies of meridional cell transports. Only velocities in the direction of the mean current of the cell are summed up. (top) The transport anomalies are shown as they result at 43°N, by integrating the deep southward branch (black, mean value of 22.5 Sv) and northward Gulf Stream branch (red, mean value of 28.5 Sv) between 50° and 35°W. Also shown is the amplitude of the MOC at the same latitude at 900-m depth (green, mean value 17.5 Sv). (bottom) The southward transports are shown from the same latitude after integrating the southward velocities between 35° and 10°W above 1000-m depth (black, mean value of 9.5 Sv) and below 1000-m depth (red, mean value of 10.0 Sv). Also shown are the anomalies of the barotropic streamfunction computed at 35°W (green, mean value of 7.0 Sv) and the NAO index (blue).

Figure 9 (bottom) shows southward transports of the eastern Atlantic roughly above and below 1000-m depth (see the caption for details). Those variations are very different from what can be seen at the western side. Rather than showing a relation to the varying strengths of the MOC, the eastern transport variations agree more with the evolution of the barotropic streamfunction. The time changes in the eastern transports can be understood as a response of the circulation to the North Atlantic Oscillation (NAO) (shown as a blue curve); this response showed a high index until 1995, which is associated with a strong barotropic circulation (Curry and McCartney 2001). Likewise, a strengthening of the barotropic circulation near the subpolar front (43°N, 35°W in our case) is described by Eden and Willebrand (2001) as the fast response to a high NAO. They showed that a strengthening of the subtropical gyre is described as the delayed response to a high NAO. While the in-phase relation due to the fast response is clearly established, the delayed mechanism is not apparent in our solution.

Table 1 shows the time-mean volume transports through key channels and passages as they result from the reference and the optimized runs. The associated transport of enthalpy is also shown for comparison with the values provided by Ganachaud and Wunsch (2000). Changes between the constrained and unconstrained solutions are in general small and within ±15% of the first-guess values, for example, for the transport through the Indonesian Archipelago and the Straits of Florida. The Straits of Florida transport flow is in the range of the observations [the classical value from the work of Schmitz and Richardson (1968) is 30 Sv] and is reduced by 1.4 Sv relative to the unconstrained run through the adjustments of the wind stress south of 24°N that lead to a lower Caribbean Current transport. The associated variabilities of both time series are very similar, showing an increase from 29 Sv in 1992 to 35 Sv in 1998 and a reduction back to 29 Sv after 2000. The increase from 1992 to 1998 is in accord with the time series of Baringer and Larsen (2001) and transport estimates from the calibration cruises (information available online at http:www.aoml.noaa.gov/phod/floridacurrent/). The former shows a slightly smaller maximum value of 34 Sv and the latter also shows the decline after 1998 but with values that are larger by about 1 Sv. The anticorrelation of the transport with the NAO index explains the closeness of the transport estimates to those of the unconstrained solution and points to very small formal estimation errors and model errors that may be dominant.

Figure 10 shows the time-mean global meridional heat and freshwater transports and their interannual variations computed as the standard deviation from individual annual mean values. A selection of values from recent independent transport estimates obtained from box inversions (Ganachaud and Wunsch 2003; Macdonald and Wunsch 1996) and from those from Wijffels et al. (1992) shows similar shapes and values, indicating the degree of reliability of recent ocean transport estimates as a legacy of WOCE. Although no formal error estimates are attached to our results (because of the computational burden associated with it), the interannual variability in our estimates (also shown in Fig. 10) indicates how variable those transports can be in the ocean. The associated variability in the heat transport (relative to the typical mean values) is higher (up to 50% of the mean near 10°S) than the corresponding values of the freshwater transport, which reach a maximum of only 20% around 5°S. The highest variability occurs in the southern equatorial region.

Fig. 10.

Global time-mean (top) meridional heat (PW) and (bottom) freshwater (Sv) transports. The envelopes represent one std dev of those quantities computed from 11 individual annual mean values. The symbols with error bars are estimates of Ganachaud and Wunsch (2000) in the top panel and of MacDonald and Wunsch (1996) (squares) and Wijffels et al. (1992) (circles) in the bottom panel.

Fig. 10.

Global time-mean (top) meridional heat (PW) and (bottom) freshwater (Sv) transports. The envelopes represent one std dev of those quantities computed from 11 individual annual mean values. The symbols with error bars are estimates of Ganachaud and Wunsch (2000) in the top panel and of MacDonald and Wunsch (1996) (squares) and Wijffels et al. (1992) (circles) in the bottom panel.

While in the past much emphasis was put on the estimation of time-mean quantities, the temporal variation of the MOC and meridional transports of heat and freshwater has recently gained much attention, especially through programs like CLIVAR and the Rapid Climate Change Programme (RAPID; information online at http://www.soc.soton.ac.uk/rapid/rapid.php). However, for the present estimate the analysis of the overturning trend is problematic because in z-coordinate models the flow over the Iceland–Scotland Ridge is associated with substantial mixing that dilutes the water masses (Beckmann and Döscher 1997). As a result, the dense overflow water is too light and enters the deep ocean south of the ridge at too shallow a depth (see also Willebrand et al. 2001). Therefore, the linear trend of the model’s meridional streamfunction reveals a reduction of the deep overturning, which is mainly associated with a shoaling of the Deep Western Boundary Current in the Atlantic. The realistic part of the trend cannot be isolated and the interpretation of the MOC trends is not straightforward. The respective result is therefore not shown from our solution.

The trend of the barotropic streamfunction, shown in Fig. 11, is wind driven and thus much less affected by the model details. It reveals a strengthening of the gyres in the Pacific Ocean whereas the north equatorial gyre slightly weakens but its zonal extension to the east increases. In the Indian Ocean the south equatorial gyre slightly strengthens and the north equatorial gyre slightly weakens. The subtropical gyre shows inconsistent trends. In agreement with Häkkinen and Rhines (2004), a general weakening is observed in the Atlantic Ocean, which has the strongest signal in the subpolar Atlantic. Although the total transport of the ACC is slightly reduced, the frontal structures within the currents intensify. According to the Sverdrup balance, changes in barotropic streamfunction can be related to changes in the wind stress curl. The drift of the wind stress curl reveals (not shown) roughly similar patterns although they include considerably smaller scales. This drift is consistent with conclusions drawn by Roemmich et al. (2007), who describe an increase of the South Pacific circulation as a result of the intensification of the wind stress curl east of New Zealand.

Fig. 11.

Linear trends of the barotropic streamfunction (Sv yr−1) estimated by least squares fitting a line to the 11-yr-long time series of monthly mean fields.

Fig. 11.

Linear trends of the barotropic streamfunction (Sv yr−1) estimated by least squares fitting a line to the 11-yr-long time series of monthly mean fields.

6. Secular heat content changes

For climate applications, our knowledge of the connection between the changing atmospheric forcing and changes in the deep ocean over hundreds of years is rudimentary. In this discussion, changes in the ocean’s inventory of heat, salt, mass, and sea level are all relevant, as is ice melting. We will use in the following the model’s results to compute changes in the global model heat and salt content. Shown in Fig. 12 are time series of the global top-to-bottom heat content of our solution together with those computed separately over the top 510 m, those from 510 to 2200 m, and those from 2200 m to the bottom. A temporal mean was removed from each individual curve.

Fig. 12.

Time series of the global heat content of our solution (blue) plotted together with those computed separately over the top 510 m (green), from 510 to 2200 m (red), and from 2200 m to the bottom (cyan). Also shown is an estimate of global heat content increase as it would if following from Levitus et al. (2005) for the 11-yr time series (thick solid line). A temporal mean was removed from all curves.

Fig. 12.

Time series of the global heat content of our solution (blue) plotted together with those computed separately over the top 510 m (green), from 510 to 2200 m (red), and from 2200 m to the bottom (cyan). Also shown is an estimate of global heat content increase as it would if following from Levitus et al. (2005) for the 11-yr time series (thick solid line). A temporal mean was removed from all curves.

The global top-to-bottom increase in heat content is about 0.75 × 1022 J yr−1 over the 11-yr period of 1992–2002. This number is close to the estimate provided by Levitus et al. (2005), as their Fig. 1 shows about 7.5 × 1022 J yr−1 from 1992 to 2003 in the upper 700 m, which corresponds to an increase of about 0.7 × 1022 J yr−1. It is interesting to note that the change in heat content of the unconstrained run is 1 × 1021, nearly an order of magnitude smaller, which indicates that the assimilation brought the heat balance of the model into consistency with the independent findings. The model results suggest that the top layer contributes about half of the total (0–700 m: 0.39 × 1022 J yr−1) to this increase, but that a significant component of the net heat storage change actually occurs below 2000-m depth. We note that in our solution the total increase in heat content occurs after the last ENSO event. In the past, temporal changes in individual depth ranges essentially compensated each other by a simultaneous cooling of the upper layer and a warming at depth. In the middle layer this warming is at first gradual and then accelerates toward the later years. In contrast, the bottom layer shows a steady increase in heat content with almost no fluctuations on shorter time scales visible as one might expect if individual atmospheric events or the NAO cycle would imprint itself in this heating of the deep ocean. It is important to note that the long-period changes discussed here, although not identical in detail, did not change totally in character between the solution here and the earlier solution on a 2° spatial grid described by SEA02. Moreover, the estimates vary only slightly between the last iterations, and we consider them therefore a somewhat robust solution of our system.

7. Sea surface height changes

Figure 13 shows the model trends in sea surface height, estimated by least squares fitting a straight line to the sea surface height fields over the period 1991–2002. The largest changes are of the order of ±2 cm yr−1 and can be found in the subpolar regions of both hemispheres. Smaller changes are visible in the tropical and subtropical oceans. All of these changes show intriguing gyre or circulation structures. Sea surface height changes are positive along a ridge across the Pacific in both hemispheres. A somewhat similar structure also seems to be present in the North Atlantic. In many locations, the altimeter-alone changes are similar to those estimated by the data–model combination. The biggest differences appear in the amplitude of the changes, rather than their spatial structures: a comparison of the two panels in Fig. 13 reveals that the model in some locations produces only about one-half (yet in others 2 times) the observed amplitudes in the regional changes in sea level (a linear trend of the globally averaged sea level increase was removed from the T/P data in the figure to be consistent with the model results that do not change on a global average because of the Boussinesq approximation of the model). Although the pattern of changes in sea level of the unconstrained run show similarities with those of the optimized run, unrealistically large changes in the higher latitudes appear to be greatly reduced and much more realistic in the optimized solution. As compared with the constrained run, the unconstrained first-guess run of the model [which included a relaxation of the surface temperature and salinity fields toward Reynolds and Smith (1994) monthly mean temperature and climatological monthly mean Levitus salinity fields] shows substantially different (both in amplitude and pattern) trends in the model SSH field. In particular, the unconstrained solution shows large trends near Antarctica that were largely (but not entirely) removed by the assimilation.

Fig. 13.

Linear least squares trend of the SSH changes (cm yr−1) model estimated from (top) mapped Aviso SSH fields, (middle) the optimization, and (bottom) the unconstrained simulation fields.

Fig. 13.

Linear least squares trend of the SSH changes (cm yr−1) model estimated from (top) mapped Aviso SSH fields, (middle) the optimization, and (bottom) the unconstrained simulation fields.

Because SSH changes on global scale are much smaller than most observed regional SSH trends, estimates of regional changes shown in Fig. 13 are not significantly affected by errors in global imbalances. To discuss what causes regional changes, we decompose in Fig. 14 the estimated model SSH trend over the 11-yr period 1992–2002 into a steric sea level change of the model potential density field and a nonsteric change in SSH, computed here as the residual between the total SSH trend and its steric component. It is obvious from the figure that most regional changes in the model SSH can be explained by steric shifts. But in several locations a significant contribution to local sea level change results from mass redistribution within the quasi-global model domain that to some degree compensates for the steric change. As an example, the Labrador Sea loses 10 cm of SSH because of mass flux out of that region. This compensates for half of the steric SSH increase so that the net effect on SSH changes in the North Atlantic is significantly smaller than a simple analysis of the temperature data would have suggested. Other regions that are significantly affected by changes in mass include the subpolar North Pacific and several places along the ACC. Moreover, the subtropical gyres of the Pacific are gaining mass in both hemispheres and are associated with a simultaneous increase in gyre strength, amplifying the steric SSH increase there. In the North Pacific the mass increase enhances the temperature-driven SSH increase; in the South Pacific the salt-driven SSH increase is enhanced through the influx of mass. Except for the shelf regions, we find, in accord with Dickey et al. (2002), an equatorward transport of mass in the 1990s whereas this is only a partially independent result since.

Fig. 14.

Model (a) steric sea level changes and (b) nonsteric SSH change estimated over the 11-yr period of 1992–2002 from the model potential density field. Corresponding (c) thermosteric and (d) halosteric sea level change estimated over the same 11-yr period. All fields are in centimeters per year, but note the different color scales.

Fig. 14.

Model (a) steric sea level changes and (b) nonsteric SSH change estimated over the 11-yr period of 1992–2002 from the model potential density field. Corresponding (c) thermosteric and (d) halosteric sea level change estimated over the same 11-yr period. All fields are in centimeters per year, but note the different color scales.

In the bottom row in Fig. 14 the linear steric sea level change is decomposed into thermosteric and halosteric contributions during 1992–2002. A comparison of the changes in SSH with those in temperature and salinity reveals that over the Pacific and North Atlantic the steric change is dominated by temperature but reduced by salinity effects whereas the opposite is true for the South Atlantic and the southern part of the Indian Ocean. Overall, partial cancellation between the temperature and salinity effects is observed in most areas, which might be explained by the conclusion that much of the exchange takes place on isopycnals. Because mixing is in general small, density can only be changed by convergences in one layer that are compensated by divergences in a different layer as well as heat and salt flux through the surface.

To understand what causes the changes in SSH that were observed by T/P, we can use the model results to provide possible mechanisms. The thermosteric effect is caused by an imbalance of the local heat flux through the surface and an associated divergence of the lateral heat transport in the underlying water column. Observed changes are consistent with a surface imbalance of less than 1 W m−2. Since variability exists on all time scales, a trend may also be a result of short extreme events, such as the 1997/98 El Niño, which explains most of the pattern between 20°S and 20°N. Outside the Tropics, the correspondence of the SSH trend patterns with those of the trend of the barotropic circulation in Fig. 11 reveals that most of the observed changes are wind driven. For a more complete answer, we have to refer to a later study involving an analysis of model and adjoint model sensitivities.

8. Discussion

The purpose of this global time-dependent data synthesis, made as a contribution of the ECCO Consortium, is to obtain a description of the time-dependent ocean circulation that is, within error bars, consistent with the global datasets and the dynamics embedded in the numerical model. We presented results from a global 11-yr-long synthesis on a 1° global grid that incorporated all major global datasets available from the period 1992 through 2002 and that as such is one of the most complete WOCE syntheses available. The resulting cost function terms demonstrate that we were able to achieve an improvement of the model simulation, although the resulting state is not yet fully consistent with all data types and prior error assumptions.

We emphasize that our optimization is performed over an 11-yr period only. Over this relatively short period, internal model errors, such as the background mixing coefficients, are less critical to the final solution than are the initial conditions and external forcing fields, although they may influence the equatorial dynamics. A detailed discussion of the remaining model–data differences supports the presumption that our results are of sufficient quality to start a first analysis of long-period variations in the estimate, discuss their relation to the real ocean, and analyze potential mechanisms that are responsible for those changes in the real ocean.

Of interest for us here are quantities such as the flow field, sea level changes, transports, mixed layer depths, overturning, and especially the surface heat, freshwater, and momentum fluxes. We note that Wunsch and Heimbach (2006) offer a somewhat similar discussion from a slightly extended version (through 2004) of our ECCO ocean-state estimation system. In summary, the mean meridional heat transport, largely underestimated in the previous 2° optimization (Stammer et al. 2003), lies now in the lower range of other recent estimates. Time changes in the model’s heat content are slightly larger than those reported by Levitus and correspond to a global net heat uptake of about 1.1 W m−2 over the model domain.

Estimates of global sea level rise are of significant public concern in the context of global change. Unfortunately, we were not able to provide such estimates since the global sea level is not a direct diagnostic of the model and has to be deduced from steric SSH contributions, which are affected by a relatively large imbalance in the global salinity budget in our estimates. It is important to recall that observed sea level variations occur not just on the global scale but also show significant regional variations. As discussed by White et al. (2005), it is the regional, especially coastal, SSH changes that have potential implications for coastal planners, not just global averages. While not very useful at the current resolution for studies of coastal SSH change, our results can be used to study changes in sea level on regional and basin scales. Because they are an order of magnitude larger than global values, they are less affected by imperfect global imbalances in surface heat or freshwater fluxes. The estimated trends in sea surface height over the period 1993–2002 is consistent overall with the observations of the TOPEX/Poseidon altimeter over most of the global ocean, although some regions show different amplitudes. Most noticeable, our estimate overestimates the SSH increase in the western subtropical Pacific.

Sea surface height changes in the model are primarily steric but show a contribution of mass redistribution in high latitudes. Steric contributions are primarily temperature based but are partly compensated by salt variations. However, the South Atlantic and the area south of the Indian Ocean reveal the dominance of salt in large-scale sea level variations. In principle, the model results can be used to provide an explanation for the changes in SSH that were observed by T/P. In practice it is not a simple matter to identify the mechanisms leading to observed changes but will require ensembles of sensitivity studies or extra optimization runs. Such experiments are beyond the scope of this paper.

In principle any estimate such as that provided here should be given with error bars. However, estimating the a posteriori error covariance matrix is very important but almost a parallel problem to the estimation problem itself. Before it is solved, the community needs to start using the existing ocean-state estimation results in scientific studies, demonstrate their usefulness and skill, or help to identify deficiencies so that they can be improved in future runs.

One difficulty in estimating climate changes in the ocean is in distinguishing between interannual-to-decadal oscillations such as ENSO fluctuations, and longer-term trends. Equally important, if we are to estimate the deep-ocean state in the presence of those climate fluctuations, a simple numerical model drift has to be separated from changes in the deep ocean that take place as a result of past climate changes. Separating real physical shifts in sea level and hydrographic properties from regional numerical drifts will be a challenge for longer-term climate-oriented state estimates. During the short estimation period of 11 yr, model adjustments are still important since some of the physics are not sufficiently resolved, especially the representation of the Greenland–Scotland Ridge overflow. To improve the ocean estimates, subsequent estimation approaches should include overflow parameterization schemes as well as bottom boundary layers and add northern boundary conditions into the control vector. In addition to those changes, several further improvements will be required in future estimates. These include an extension of the estimation period to 50 yr, paralleling the NCEP–NCAR reanalysis. Moreover, the model needs to be truly global and needs to include an ice model, and the limitations of the Boussinesq approximation need to be corrected. Then we can expect quantitative answers to more of the issues related to sea level rise and decadal ocean variability.

Table A1. Prior error description, used as the model–data misfit weighting factor in Eq. (1).

Table A1. Prior error description, used as the model–data misfit weighting factor in Eq. (1).
Table A1. Prior error description, used as the model–data misfit weighting factor in Eq. (1).

Acknowledgments

We thank all members of the ECCO Consortium for their direct or indirect support and discussions that lead to the results that are presented in this paper. In particular, Elisabeth Remy, Christian Eckert, and Patrick Heimbach were instrumental for setting up the ECCO 1° assimilation infrastructure. Youyu Lu, Charmaine King, Diana Spiegel, and Kyozo Ueyoshi helped with the processing of the input data. We also thank Lee Fu and two anonymous reviewers who suggested many improvements to the manuscript. Don Chambers from the University of Texas provided the T/P-GRACE dynamic topography (based on the GRACE model GGM01S). Reanalysis surface forcing fields from NCEP–NCAR were obtained through a computational grant at NCAR. SSH drifts were estimated from the Ssalto/Duacs maps, produced as part of the Environment and Climate EU Enact project (EVK2-CT2001-00117) and distributed by Aviso, with support from CNES. Computational support from the National Partnership for Computational Infrastructure (NPACI) and NCAR is acknowledged. Support was provided in part through ONR (NOPP) ECCO Grants N00014-99-1-1049, NAG5-11765, NAG5-12870, and NNG04GF30G, and through a contract with the Jet Propulsion Laboratory (1205624). This is a contribution of the Consortium for Estimating the Circulation and Climate of the Ocean (ECCO) funded by the National Oceanographic Partnership Program.

This paper is a tribute to Carl Wunsch and his long-held vision of the way the ocean state and its changes should be diagnosed through a modern observing system and through a quantitative model–data combination. One of his outstanding contributions to physical oceanography is his pioneering development of inverse methods in oceanography that today enables the community to perform routine ocean-state estimation. This is equivalent to weather prediction for the atmosphere and is intended to support operational oceanography as well as climate research. Carl had not only this vision but also the energy to move the field forward from “business as usual” into a new way of thinking about the ocean and how it should be observed and modeled—a vision that finally led to the establishment of a consortium among MIT, the Jet Propulsion Laboratory, and the Scripps Institution of Oceanography: Estimating the Circulation and Climate of the Ocean. It has also led to the ongoing international Global Ocean Data Assimilation Experiment (GODAE).

Thinking about Carl’s work, many other pioneering contributions to physical oceanography and to our understanding of how the ocean works come to mind. His role in the conception of the World Ocean Circulation Experiment (WOCE) is well known. An important observational component of WOCE was the U.S.–French satellite altimeter mission, called TOPEX/Poseidon. Perhaps some do not know that this satellite mission goes back to the vision of a few, including Carl, who understood the importance of having continuous coverage of satellite altimeter observations of the ocean: a satellite is the only tool that can observe ocean variability on all relevant time and space scales. History now has proven the concept to be correct by making the TOPEX/Poseidon mission the longest and most successful NASA mission to date and establishing altimetry as an important component of ongoing and future climate observing systems.

Beyond his own important scientific contributions and his leadership, Carl’s impact can also be measured by the number of students and young scientists he has advised since 1966 or had influenced significantly during their careers. From our experience we can say that Carl has the unique ability to stimulate and entrain young scientists into cutting-edge research discussions, which in many ways is as important as his outstanding expertise in the field. For many peers, Carl’s stimulating enthusiasm about their work has helped to make significant progress and evolution in the field. It is this sparkle, as much as his scientific accomplishments, that demonstrates Carl’s truly remarkable leadership.

REFERENCES

REFERENCES
Adcroft
,
A.
,
J-M.
Campin
,
P.
Heimbach
,
C.
Hill
, and
J.
Marshall
,
cited
.
2002
:
Mitgcm release 1. [Available online at http://mitgcm.org/sealion/.]
.
Antonov
,
J. I.
,
S.
Levitus
, and
T. P.
Boyer
,
2002
:
Steric sea level variations during 1957–1994: Importance of salinity.
J. Geophys. Res.
,
107
.
8013, doi:10.1029/2001JC000964
.
Arakawa
,
A.
, and
V. R.
Lamb
,
1977
:
Computational design of the basic dynamical processes of the UCLA general circulation model. General Circulation Models of the Atmosphere, J. Chang, Ed., Methods of Computational Physics, Vol. 17, Academic Press, 173–265
.
Baringer
,
M. O.
, and
J. C.
Larsen
,
2001
:
Sixteen years of Florida current transport at 27°N.
Geophys. Res. Lett.
,
28
,
3179
3182
.
Beckmann
,
A.
, and
R.
Döscher
,
1997
:
A method for improved representation of dense water spreading over topography in geopotential-coordinate models.
J. Phys. Oceanogr.
,
27
,
581
591
.
Bennett
,
A. F.
,
2002
:
Inverse Modeling of the Ocean and Atmosphere.
Cambridge University Press, 352 pp
.
Bryden
,
H.
,
H. R.
Longworth
, and
S. A.
Cunningham
,
2005
:
Slowing of the Atlantic meridional overturning circulation at 25°N.
Nature
,
438
,
655
657
.
Cazenave
,
A.
, and
R. S.
Nerem
,
2004
:
Present-day sea level change: Observations and causes.
Rev. Geophys.
,
42
.
RG3001, doi:10.1029/2003RG000139
.
Cooper
,
M.
, and
K.
Haines
,
1996
:
Altimetric assimilation with water property conservation.
J. Geophys. Res.
,
101
,
1059
1077
.
Curry
,
R. G.
, and
M. S.
McCartney
,
2001
:
Ocean gyre circulation changes associated with the North Atlantic Oscillation.
J. Phys. Oceanogr.
,
31
,
3374
3400
.
Daley
,
R.
,
1991
:
Atmospheric Data Analysis.
Cambridge University Press, 457 pp
.
Dickey
,
J. O.
,
S. L.
Marcus
,
O.
de Viron
, and
I.
Fukumori
,
2002
:
Recent earth oblateness variations: Unraveling climate and postglacial rebound effects.
Science
,
298
,
1975
1977
.
Eden
,
C.
, and
J.
Willebrand
,
2001
:
Mechanism of interannual to decadal variability of the North Atlantic circulation.
J. Climate
,
14
,
2266
2280
.
Ganachaud
,
A.
, and
C.
Wunsch
,
2000
:
Improved estimates of global ocean circulation, heat transport and mixing from hydrographic data.
Nature
,
408
,
453
457
.
Ganachaud
,
A.
, and
C.
Wunsch
,
2003
:
Large-scale ocean heat and freshwater transports during the World Ocean Circulation Experiment.
J. Climate
,
16
,
696
705
.
Gent
,
P. R.
, and
J. C.
McWilliams
,
1990
:
Isopycnal mixing in ocean models.
J. Phys. Oceanogr.
,
20
,
150
155
.
Giering
,
R.
, and
T.
Kaminski
,
1998
:
Recipes for adjoint code construction.
Assoc. Comput. Mach. Trans. Math. Software
,
24
,
437
474
.
Gilbert
,
J. C.
, and
C.
Lemaréchal
,
1989
:
Some numerical experiments with variable-storage quasi-Newton algorithms.
Math. Programm.
,
45
,
407
435
.
Häkkinen
,
S.
, and
P. B.
Rhines
,
2004
:
Decline of subpolar North Atlantic circulation during the 1990s.
Science
,
304
,
555
559
.
Kalnay
,
E.
, and
Coauthors
,
1996
:
The NCEP/NCAR 40-Year Reanalysis Project.
Bull. Amer. Meteor. Soc.
,
77
,
437
471
.
Käse
,
R. H.
, and
W.
Krauss
,
1996
:
The Gulf Stream, the North Atlantic Current, and the origin of the Azores Current. The Warmwatersphere of the North Atlantic Ocean, W. Krauss, Ed., Gebrüder Borntraeger, 291–337
.
Köhl
,
A.
,
Y.
Lu
,
P.
Heimbach
,
B.
Cornuelle
,
D.
Stammer
, and
C.
Wunsch
, and
the ECCO Consortium
,
2002
:
The ECCO 1° global WOCE synthesis. ECCO Rep. 20, 33 pp. [Available online at http://www.ecco-group.org/reports.html]
.
Large
,
W. G.
,
J. C.
McWilliams
, and
S. C.
Doney
,
1994
:
Oceanic vertical mixing: A review and a model with nonlocal boundary layer parameterization.
Rev. Geophys.
,
32
,
363
403
.
Levitus
,
S.
,
R.
Burgett
, and
T.
Boyer
,
1994a
:
Salinity. Vol. 3, World Ocean Atlas 1994, NOAA Atlas NESDIS 3, 99 pp
.
Levitus
,
S.
,
R.
Burgett
, and
T.
Boyer
,
1994b
:
Temperature. Vol. 4, World Ocean Atlas 1994, NOAA Atlas NESDIS 4, 117 pp
.
Levitus
,
S.
,
J. I.
Antonov
,
J. L.
Wang
,
T. L.
Delworth
,
K. W.
Dixon
, and
A. J.
Broccoli
,
2001
:
Anthropogenic warming of Earth’s climate system.
Science
,
292
,
267
270
.
Levitus
,
S.
,
J. J.
Antonov
, and
T.
Boyer
,
2005
:
Warming of the World Ocean, 1955–2003.
Geophys. Res. Lett.
,
32
.
L02604, doi:10.1029/2004GL021592
.
Lu
,
Y.
,
K.
Ueyoshi
,
A.
Köhl
,
E.
Remy
,
K.
Lorbacher
, and
D.
Stammer
,
2002
:
Input data sets for the ECCO Global 1° WOCE synthesis. ECCO Rep. 18, 37 pp
.
Macdonald
,
A. M.
, and
C.
Wunsch
,
1996
:
An estimate of global ocean circulation and heat fluxes.
Nature
,
382
,
436
439
.
Maltrud
,
M. E.
,
R. D.
Smith
,
A. J.
Semtner
, and
R. C.
Malone
,
1998
:
Global eddy-resolving ocean simulations driven by 1985–1995 atmospheric winds.
J. Geophys. Res.
,
103
,
30825
30853
.
Marotzke
,
J.
,
1992
:
The role of integration time in determining a steady-state through data assimilation.
J. Phys. Oceanogr.
,
22
,
1556
1567
.
Marotzke
,
J.
, and
C.
Wunsch
,
1993
:
Finding the steady-state of a general-circulation model through data assimilation–Application to the North Atlantic Ocean.
J. Geophys. Res.
,
98
,
20149
20167
.
Marotzke
,
J.
,
R.
Giering
,
Q. K.
Zhang
,
D.
Stammer
,
C. N.
Hill
, and
T.
Lee
,
1999
:
Construction of the adjoint MIT ocean general circulation model and application to Atlantic heat transport sensitivity.
J. Geophys. Res.
,
104
,
29529
29548
.
Marshall
,
J.
,
C.
Hill
,
L.
Perelman
, and
A.
Adcroft
,
1997a
:
Hydrostatic, quasi-hydrostatic, and nonhydrostatic ocean modeling.
J. Geophys. Res.
,
102C
,
5733
5752
.
Marshall
,
J.
,
A.
Adcroft
,
C.
Hill
,
L.
Perelman
, and
C.
Heisey
,
1997b
:
A finite-volume, incompressible Navier–Stokes model for studies of the ocean on parallel computers.
J. Geophys. Res.
,
102C
,
5753
5766
.
Menemenlis
,
D.
, and
Coauthors
,
2005
:
NASA supercomputer improves prospects for ocean climate research.
Eos, Trans. Amer. Meteor. Soc.
,
86
,
95
96
.
Miller
,
L.
, and
B. C.
Douglas
,
2004
:
The rate of twentieth-century global sea level rise and its causes are the subjects of intense controversy. Most direct estimates from tide gauges give.
Nature
,
428
,
406
409
.
Munk
,
W.
,
2003
:
Ocean freshening, sea level rising.
Science
,
300
,
2041
2043
.
Niiler
,
P.
,
2001
:
The World Ocean surface circulation.
Ocean Circulation and Climate, G. Siedler, J. Church, and J. Gould, Eds., International Geophysics Series, Vol. 77, Academic Press, 193–204
.
Reynolds
,
R. W.
, and
T. M.
Smith
,
1994
:
Improved global sea surface temperature analysis using optimum interpolation.
J. Climate
,
7
,
929
948
.
Roemmich
,
D.
,
J.
Gilson
,
R.
Davis
,
P.
Sutton
,
S.
Wijffels
, and
S.
Riser
,
2007
:
Decadal spinup of the South Pacific subtropical gyre.
J. Phys. Oceanogr.
,
37
,
162
173
.
Schmitz
,
W. J.
, and
W. S.
Richardson
,
1968
:
On the transport of the Florida Current.
Deep-Sea Res.
,
15
,
679
693
.
Smith
,
W. H. F.
, and
D. T.
Sandwell
,
1997
:
Global seafloor topography from satellite altimetry and ship depth soundings.
Science
,
277
,
195
196
.
Stammer
,
D.
,
1998
:
On eddy characteristics, eddy transports, and mean flow properties.
J. Phys. Oceanogr.
,
28
,
727
739
.
Stammer
,
D.
,
C.
Wunsch
, and
R. M.
Ponte
,
2000
:
De-aliasing of global high frequency barotropic motions in altimeter observations.
Geophys. Res. Lett.
,
27
,
1175
1178
.
Stammer
,
D.
,
C.
Wunsch
,
I.
Fukumori
, and
J.
Marshall
,
2002a
:
State estimation improves prospects for ocean research.
Eos, Trans. Amer. Geophys. Union
,
83
.
289, 294, 295
.
Stammer
,
D.
, and
Coauthors
,
2002b
:
The global ocean circulation during 1992–1997, estimated from ocean observations and a general circulation model.
J. Geophys. Res.
,
107
.
3118, doi:10.1029/2001JC000888
.
Stammer
,
D.
,
F.
Wentz
, and
C.
Gentemann
,
2003
:
Validation of new microwave SST measurement for climate purposes.
J. Climate
,
16
,
73
87
.
Stammer
,
D.
,
K.
Ueyoshi
,
A.
Köhl
,
W. B.
Large
,
S.
Josey
, and
C.
Wunsch
,
2004
:
Estimating air–sea fluxes of heat, freshwater, and momentum through global ocean data assimilation.
J. Geophys. Res.
,
109
.
C05023, doi:10.1029/2003JC002082
.
Stammer
,
D.
,
A.
Köhl
, and
C.
Wunsch
,
2007
:
Impact of accurate geoid fields on estimates of the ocean circulation.
J. Atmos. Oceanic Technol., in press
.
Tapley
,
B. D.
,
D. P.
Chambers
,
S.
Bettadpur
, and
J. C.
Ries
,
2003
:
Large scale ocean circulation from the GRACE GGM01 geoid.
Geophys. Res. Lett.
,
30
.
2163, doi:10.1029/2003GL018622
.
White
,
N. J.
,
J. A.
Church
, and
J. M.
Gregory
,
2005
:
Coastal and global averaged sea level rise for 1950 to 2000.
Geophys. Res. Let.
,
32
.
L01601, doi:10.1029/2004GL021391
.
Wijffels
,
S. E.
,
R. W.
Schmitt
,
H. L.
Bryden
, and
A.
Stigebrandt
,
1992
:
Transport of freshwater by the oceans.
J. Phys. Oceanogr.
,
22
,
155
162
.
Willebrand
,
J.
, and
Coauthors
,
2001
:
Circulation characteristics in three eddy-permitting models of the North Atlantic.
Progress in Oceanography
,
48
,
Pergamon,
.
123
161
.
Wunsch
,
C.
,
1996
:
The Ocean Circulation Inverse Problem.
Cambridge University Press, 442 pp
.
Wunsch
,
C.
, and
P.
Heimbach
,
2006
:
Estimated decadal changes in the North Atlantic meridional overturning circulation and heat flux 1993–2004.
J. Phys. Oceanogr.
,
36
,
2012
2024
.

APPENDIX

Description of the Cost Function

The ECCO adjoint model is obtained from the forward code by using the Transformation of Algorithms in FORTRAN (TAF) Model Compiler (Giering and Kaminski 1998; see also www.fastopt.de online) as described by Marotzke et al. (1999). The adjoint model is used to provide the gradient of the model–data misfit, J, with respect to the control variables; this gradient is used to modify the control variables, u(t), so as to reduce the value of J. Details are described in this appendix. Because of the large number of control variables, we solved the large set of simultaneous coupled nonlinear equations [“normal equations”; see Wunsch (1996) for details] iteratively, using a quasi-Newton descent algorithm (Gilbert and Lemaréchal 1989).

The nonlinearity of the KPP and GM schemes causes the adjoint related to these parameterizations to develop unstable modes that grow rapidly in time. For short time scales (shorter than a few months) the amplitudes of these modes remain very small. Within this period, the inspection of changes in the gradients after the removal of the adjoint to the parameterizations revealed only minor changes due to the code omission, which enables the exclusion of the respective adjoint codes. Over longer times the instability of the adjoint model to the KPP scheme becomes dominant but describes only the growth of the unstable modes that are damped in the full nonlinear forward model. A certain capability of the adjoint model is, however, lost. For instance, the adjoint to the KPP model describes forcing changes that would affect the vertical mixing coefficients, which is only a second-order effect. We note that only the performance of the optimization is affect but not the forward model.

The explicit form of the cost function J reads as follows:

 
formula

In (A1), T and S stand for monthly averages, T and S are climatological monthly T and S fields (i.e., averages of all model Januarys, Februarys, etc.), δT0 and δS0 are changes in the initial conditions, and the remaining δ terms represent 2-day averages of the surface flux fields. The terms T, S, and with the subscripts “first” and “last” are annual mean values over the respective years. The T/P and European Remote Sensing Satellite (ERS) altimeter anomalies are evaluated on a daily basis with along-track data averaged over 1° grid cells, and the time-mean T/P SSH field minus the GRACE GGM01S geoid model (Tapley et al. 2003), ηTP, is imposed over the entire 11-yr period.

The weighted matrices in each term of J determine the solution to the optimization problem. In principle, one should specify the full a priori error covariance matrix for each data type. The mean Levitus hydrography terms were weighted by error profiles [σT(z) and σS(z)] provided by Levitus et al. (1994a, b) and ranging in the vertical direction from 0.5° near the surface to 0.05°C at depth for potential temperature, and 0.13 to 0.01 PSU for salinity. To constrain in situ data along hydrographic sections or from XBT and float programs, the weight reflects the measurement errors and a representation error of the model that describes unresolved mesoscale features in the profile data. To construct a representation error, we used the Cooper and Haines (1996) relation based on the rms SSH variability and Levitus climatological temperature data. However, an evaluation of XBT data from the Pacific revealed that the Cooper and Haines eddy amplitudes would underestimate the actual temperature and salinity eddy variability by roughly a factor of 5 and we enhanced those values accordingly (in reality this is regionally varying). The final error weight for the climatological and profile hydrographic data is 0.25 times the inverse of the sum of the squares of the representation error and the background error. We thereby take into account that not every layer is statistically independent and that only a few vertical modes are necessary to describe the ocean. See Marotzke (1992) and Marotzke and Wunsch (1993) for a detailed discussion. To obtain a simple approximation to the weight for the temperature and salinity drift, the globally averaged Reynolds SST drift over the 11-yr ΔT was extrapolated into depth using the error profiles provided by Levitus et al. according to [ΔT/σT(1)]σT(z) and [ΔT/σT(1)]σS(z), respectively.

Imposed weighting matrices and their sources are summarized in Table A1. For the time-mean SSH field, a geographically uniform error of 4.5 cm was imposed, which accounts for uncertainties in the geoid and in the time-mean altimetric measurements. This error was provided in agreement with the GRACE project. The SST error is uniformly 0.5°C, in agreement with Reynolds and Smith (1994) and Stammer et al. (2003). Errors of mean drifter velocities are specified as provided by Niiler (2001). We note that some of the error weights were adjusted during the optimization; for example, the SST data were first excluded from shallow regions. Also, because of a coding error, the weight on the in situ data was slightly stronger than the prescribed error would imply.

We expect the model to have the least skill in simulating observed observations over shallow regions, partially because the model resolution is too low to properly simulate the coastal dynamics. All data were therefore excluded from the cost function in regions where the bathymetry is shallower than 1000 m. However, this data omission resulted in poorly constrained SST values in those regions. SST observations from those shallow regions were therefore included in the cost function during the last eight iterations. During those few iterations the solution did not fully correct the coastal SST fields, which resulted in enhanced SST misfits over some shelf regions.

Footnotes

Corresponding author address: Detlef Stammer, Institut für Meereskunde, Zentrum für Meeres- und Klimaforschung, Universität Hamburg, D-20146 Hamburg, Germany. Email: stammer@ifm.uni-hamburg.de

This article included in the In Honor of Carl Wunsch special collection.