The German partner of the consortium for Estimating the Circulation and Climate of the Ocean (GECCO) provided a dynamically consistent estimate of the time-varying ocean circulation over the 50-yr period 1952–2001. The GECCO synthesis combines most of the data available during the entire estimation period with the ECCO–Massachusetts Institute of Technology (MIT) ocean circulation model using its adjoint. This GECCO estimate is analyzed here for the period 1962–2001 with respect to decadal and longer-term changes of the meridional overturning circulation (MOC) of the North Atlantic. A special focus is on the maximum MOC values at 25°N. Over this period, the dynamically self-consistent synthesis stays within the error bars of H. L. Bryden et al., but reveals a general increase of the MOC strength. The variability on decadal and longer time scales is decomposed into contributions from different processes. Changes in the model’s MOC strength are strongly influenced by the southward communication of density anomalies along the western boundary originating from the subpolar North Atlantic, which are related to changes in the Denmark Strait overflow but are only marginally influenced by water mass formation in the Labrador Sea. The influence of density anomalies propagating along the southern edge of the subtropical gyre associated with baroclinically unstable Rossby waves is found to be equally important. Wind-driven processes such as local Ekman transport explain a smaller fraction of the variability on those long time scales.
The meridional overturning circulation (MOC) of the ocean carries a large amount of heat northward in the Atlantic that subsequently is released to the atmosphere in mid- and high latitudes and thereby significantly influences the climate in northern Europe (Hall and Bryden 1982). The fundamental importance of MOC variations to Northern Hemisphere climate changes and possible impacts on society are of paramount interest, especially since recent model studies from Manabe and Stouffer (1993); Wood et al. (1999); Cubasch et al. (2001) hypothesized that the anthropogenic increase in the atmospheric concentration of CO2 and other greenhouse gases could lead to a measurable reduction in MOC strength in the Atlantic over the next century. The possibility that such large and rapid transitions of the MOC could actually occur over the next few decades in response to present-day climate shifts is still vigorously debated in the literature, without any obvious conclusion at hand. An analysis of the variability of the MOC in the North Atlantic over the last few decades and the underlying causes and mechanisms is therefore urgently needed before the likelihood for abrupt changes can be understood.
In the past, changes of the MOC strengths were repeatedly diagnosed from observations available at several latitudes in the North Atlantic. As an example, Koltermann et al. (1999) estimated changes of the meridional volume transport in different density classes from hydrographic observations taken repeatedly between 1957 and 1993 at latitudes 24°, 36°, and 48°N. While the authors could not find any significant MOC change at 24°N between 1957 and 1993, an intensification of the MOC appeared in their results north of this latitude in 1981. Estimates of MOC changes at 25°N in the North Atlantic were investigated again recently by Bryden et al. (2005) who extended the previous time series through 2004. The authors confirmed the absence of clear changes in the MOC strength before 1992. However, they report a significant reduction in MOC strength of about 30% in subsequent years, which, in their estimates, are compensated by a corresponding increase in the horizontal gyre circulation.
To determine long-term trends from noisy measurements and to separate them from low frequency variability has often proven to be difficult and error prone. Accordingly, the MOC analysis of Bryden et al. (2005) caused criticism because it is based on just a few realizations of single repeated cross-basin sections and could therefore be aliased by short-term baroclinic variability of the MOC (Thurnherr and Speer 2004). As a result, the estimated trend may not be statistically significant. The relevance of the later problem especially is now being addressed with major dedicated observational programs in Europe and the United States aiming at continuously monitoring and detecting changes in the MOC with a much higher time resolution (e.g., National Environment Research Council (NERC) RAPID; http://www.noc.soton.ac.uk/rapidmoc/). First results by Cunningham et al. (2007) demonstrate a large variability of the flow of 18.7 ± 5.6 Sv (Sv ≡ 106 m3 s−1), and other results suggest that several decades of measurements are needed before trends can be detected with any statistical significance (Baehr et al. 2007).
If used only in isolation, different sources of information may provide inconsistent answers to the question whether the MOC is changing. For a most reliable estimate of MOC changes, all available observations that potentially contribute information to the detection of MOC changes and information about their uncertainties should be used. The ideal approach is therefore to combine—in a dynamically consistent way—all available ocean data, including the existing satellite datasets with a numerical model and use the resulting model output as the basis for studies of the MOC. The consortium for Estimating the Circulation and Climate of the Ocean (ECCO) previously demonstrated the feasibility of such a dynamically consistent global ocean state estimate over more than a decade and discussed the skill of the estimated three-dimensional, time-varying oceanic state over the period 1992–2002 (Stammer et al. 2002, 2003, 2004; Köhl et al. 2007b; Wunsch and Heimbach 2006). This approach is, however, limited by the realism of the numerical model, and unresolved processes may bias the estimate. Köhl et al. (2007b), for instance, refrained from analyzing MOC trends in their 11-yr-long solution because trends in such a relatively short run are significantly influenced by numerical (spinup) artifacts of the model.
To address climate-relevant decadal and longer-term changes of the circulation, longer estimation efforts are required. It was therefore attempted by the German partner of the ECCO effort (GECCO) to estimate the ocean circulation over the 50-yr period 1952–2002 (hereafter the GECCO run). The aim of this GECCO effort is to bring an ocean circulation model into consistency with all in situ and satellite observations that were collected since the beginning of the 1950s. Based on this new 50-yr optimization, we investigate here decadal changes of the meridional overturning in the Atlantic Ocean and also provide details on the optimization setup. In analyzing MOC changes and their dynamical causes we essentially follow the previous study of Köhl (2005) who identified mechanisms that can potentially influence the MOC at 30°N in the Atlantic. Our aim here is to quantify the relative importance of each of these mechanisms in contributing to the observed MOC variability in the Atlantic Ocean. In a related study, Köhl and Stammer (2008) investigate regional and global sea level changes on decadal and longer time scales as they result from the same 50-yr GECCO estimate.
The structure of the remaining paper is as follows: In section 2 we present the method and approach taken here. The variability of the simulated MOC variability is investigated in section 3. Estimates of changes of the Atlantic MOC from the GECCO results are presented in section 4 for the last 50 years. Mechanisms for the observed and simulated SSH changes are discussed in section 5. A decomposition of the variability into the contributing processes is attempted in section 6, and concluding remarks follow in section 7.
2. Methodology and approach
a. Adjoint model framework
The assimilation approach used in this study is essentially identical to the 11-yr ECCO synthesis on a 1° horizontal global grid, described in detail by Köhl et al. (2007b). It is based on the ECCO–Massachusetts Institute of Technology (MIT) GCM and its adjoint and covers the global ocean from 80°S to 80°N with realistic topography based on the Smith and Sandwell (1997) dataset. Vertical mixing is represented by the K-Profile Paremeterization (KPP) scheme of Large et al. (1994). Background coefficients of vertical diffusion and viscosity are 10−5 m2 s−1 and 10−3 m2 s−1, respectively. Harmonic diffusion along neutral surfaces according to Redi (1982) in association with the parameterization of Gent and Mc Williams (1990, hereafter GM) for eddy-induced tracer advection represents unresolved eddy processes. The GM advection coefficient was set to 103 m2 s−1, and the isopycnal diffusivity and horizontal viscosity were chosen to be 102 m2 s−1 and 104 m2 s−1, respectively.
The model is operated in the hydrostatic configuration with an implicit free surface. It was forced with once-per-day surface heat flux and salinity flux as estimated from precipitation and latent heat flux, obtained from the National Centers for Environmental Prediction–National Center for Atmospheric Research (NCEP–NCAR) reanalysis project. Daily shortwave heat flux is treated separately. The vertical profile of absorption is modeled by the analytic formula of Paulson and Simpson (1977) for prescribed ocean water–Type Ib after Jerlov (1968). Wind stress forcing fields were provided twice per day.
Our solution is obtained by fitting the model simultaneously to all available data over the 50-yr period. As in the previous 11-yr optimization, this is done iteratively by first running the forward model to calculate the model data misfit formulated as a cost function, followed by an adjoint model run to calculate the gradients of this cost function. Similar to our experience with the previous 11-yr optimization, the adjoint code to the GM and KPP had to be excluded from the adjoint model because this code otherwise would cause exponentially growing adjoint variables. However, in contrast to the previous optimization, this reduced (approximate) adjoint still shows exponentially growing adjoint variables when integrated over the full 50-yr period.
Köhl and Willebrand (2002) related the appearance of exponentially growing adjoint variables to an exponentially growing number of secondary minima of the cost function, which prevent the optimization from converging toward a global minimum. A cost function defined in terms of statistical quantities will converge for long integration periods to one with a parabolic shape in which the secondary minima are transformed into stochastic perturbations. The authors demonstrated that adjoint gradients, calculated by using a modified coarse-resolution adjoint of the high-resolution forward model, provide a good approximation to gradients of this parabola and will enable a convergence of an optimization into the global minimum.
To overcome the problem of adjoint instabilities during the 50-yr optimization, we followed, to some extent, the strategy of Köhl and Willebrand (2002). However, rather than reducing the resolution of the adjoint model, we increased its diffusivity values (by a factor of 100), which likewise prevented the adjoint model from developing unstable modes without affecting the large-scale (size of typical atmospheric scales) patterns of the gradient. The gradients are thus smoother than in the previous shorter run, which basically means that only the large scale could be improved. At the same time, the forward model remained unchanged and also contains the full model physics. This strategy was also successfully adopted for an eddy-resolving assimilation in the tropical Pacific (Hoteit et al. 2005).
To bring the global 1° model into agreement with observations, initial temperature and salinity conditions as well as the time-dependent (10-day averages) surface fluxes of momentum, heat, and freshwater were adjusted by the adjoint method. The data constraints available after 1992 are the same as those used during the 11-yr optimization described by Köhl et al. (2007b), and before 1992 they essentially consist of an extensive database of upper ocean thermal (XBT and MBT) measurements. However, most available tide gauge and CTD data were used as well. A schematic of the optimization in terms of the temporal availability of all input datasets and in terms of control parameters is provided in Fig. 1.
Results presented below are based on iteration 23. Over those number of iterations the cost function was reduced to about 10% of its initial value, and further progress was slow. Because of its complexity (we estimate 251 947 317 parameters from 571 271 669 observations) the optimization did not converge completely after this number of iterations, but the model was brought into a level of agreement with the data that is comparable to the previous 11-yr optimization. The higher (by 30%) rms difference to the data can be rationalized through a reduced controllability of the state by the initial conditions. The effect of model drift due to unresolved physics (mesoscale eddies) is also larger over a longer integration period (Ferreira et al. 2005).
Although formally not yet fully converged, the solution is close enough (in particular, in better agreement with data than a simulation) to investigate MOC variability and underlying processes and causes. This is legitimate as long as our focus is primarily on potential processes and less on describing the actual changes that took place in the ocean. We believe that conclusions drawn about the analyzed variability, as well as about processes and mechanism, are more affected by the limitations of the numerical model than by lack of full convergence. This point will be further substantiated in the discussion below.
We note that results shown here (like for the earlier optimization) are based on two forward runs of the ECCO–MIT model. In the first case we use the estimated initial model temperature and salinity fields and the time-varying surface forcing, which were estimated in the optimization procedure. This run will be called GECCO run below. In an unconstrained reference run, we use Levitus initial conditions and NCEP surface forcing fields (including relaxation to Reynolds surface temperature and Levitus surface salinity), as it is traditional for process-oriented model runs. In both cases monthly mean averages of the time-varying model state were used as the basis of our investigation.
3. Atlantic MOC variability
In the following, we will use the GECCO results to infer changes in the MOC of the Atlantic on interannual to interdecadal time scales. To begin, Fig. 2a shows the time mean MOC of the Atlantic Ocean. A maximum of 18 Sv at about 48°N and values around 13 Sv farther south are in the range of previous estimates, although they are somewhat lower and higher than results published by Ganachaud and Wunsch (2000) for the lower latitudes and at 48°N, respectively.
In comparison to Bryden et al. (2005) or the much earlier estimate by Roemmich and Wunsch (1984), the GECCO circulation is shallower as indicated by the depth of the MOC maximum, which was already true for the shorter ECCO solution of Köhl et al. (2007b). Accordingly, the Upper and Lower North Atlantic Deep Water (UNADW and LNADW, respectively) form one single core, which is a common problem of level models. The reasons for this are described in detail by Willebrand et al. (2001), who explain that this is not an inherent defect of the model but a consequence of the excessive mixing associated with the overflows, which causes the associated part of the MOC branch to appear at shallower depth. Accordingly, Döscher et al. (1994) was able to reproduce the lower branch of the NADW with the same type of model by providing the correct density south of the overflows. A further option is the representation of the processes in the bottom boundary layer (BBL) by a parameterization (Döscher and Beckmann 2000). However, Eden (1999) showed in a model very similar to ours that the effect of this parameterization on the MOC is essentially confined to the layers below 1000 m, while the maximum MOC is not much affected. Moreover, Schweckendiek and Willebrand (2005) found that the MOC evolution in regional ocean models corresponds closely to that of different coupled simulations (one with an isopycnal and one with a level ocean model without BBL parameterization) and therefore does not depend on the coupling or detailed representation of the ocean processes. It then can be rationalized that. despite the deficiencies of the simulated MOC structure, the variability of the maximum MOC can be studied. However, if details of deeper layers are to be considered, one has to keep in mind that transports appear at lighter densities. Since the model adjusts during the first decade to the deficiencies in the production of sufficiently dense water, in all analyses we disregard the first decade of the simulation.
A decomposition of the MOC into the geostrophic, Ekman, and residual components is shown in Fig. 2. The Ekman component is calculated from the optimized zonal wind stress, assuming that a depth-independent compensation of the surface flow holds for which Kanzow et al. (2007) found observational indications. A further separation into the total barotropic and geostrophic vertical shear component [as performed by Lee and Marotzke (1998), not shown] leads to a barotropic component (created by the barotropic circulation flowing over zonally nonuniform topography), which is found in agreement with Sime et al. (2006) to dominate the MOC between 20° and 30°N. However, the thermal wind shear compensates for a large fraction of the external component in this latitude band. As explained in section 5, in the North Atlantic the two components cannot be regarded as independent. We therefore consider in the following only the sum of both components. To this end, the geostrophic velocity vg is determined by the hydrostatic pressure calculated from density ρ and the SSH elevation ζ, representing the ideal case in which the reference level is known:
where g is the gravity, ρ0 a reference density, and f the Coriolis parameter. A small spatially uniform velocity O(0.1 mm s−1) was removed afterward to insure zero meridional transport.
Outside the tropics (25°S–25°N) the flow appears primarily in geostrophic balance, while between both latitudes the Ekman part gains importance. The residual flow is collocated with the geostrophic component, but has a deeper maximum at about 1500 m.
The largest fraction is associated with an underestimation of the flow next to lateral or bottom boundaries. Friction (horizontal as well as vertical) and nonlinear terms in the momentum equation are one to two orders of magnitude too small to provide an explanation for the residual flow. However, we recall that the calculation of geostrophic velocities is also sensitive to numerical details (velocity and density grids are shifted on a C grid) and therefore is easily error prone. The residual flow thus represents mainly the error involved in the calculation of the geostrophic component and has no direct physical implication. The residual overturning of up to 5 Sv apparent north of 30°N is found to originate basically from the western boundary current region.
Shown in the left column of Fig. 3 is the standard deviation (STD) of the total MOC variability (based on monthly mean fields) and its contribution on subannual time scale (deviation from annual means) and on interannual and longer time scales (variability of annual mean fields). The total STD shows a maximum of more than 10 Sv at depths shallower than 500 m just north of the equator. In other regions the maximum STD is located around 1000 m with values close to 2.5 Sv (rms). The amplitude of the subannual component is only slightly smaller and explains most of the variability. In comparison, the annual component shows only values of around 1 Sv (rms) with a maximum at 1000 m. The separation of the full variability into Ekman, geostrophic, and residual components (right column) indicates that most of the shallow MOC variability is caused by the Ekman component, while the deeper (below 500 m) fraction seems to be mainly geostrophic in character. Finally, the STD of the residual MOC is small (less than 1 Sv outside the tropics). The residual MOC shows a vertically coherent structure as it would emerge from the transport through the bottom cells, indicating a large sensitivity of theses computations to details on the numerics, particularly near the boundary points.
4. MOC changes at 25°N
To facilitate a direct comparison with recent results provided by Bryden et al. (2005) for 25°N, we show in Fig. 4 a time series of the maximum MOC at 25°N (here and in all subsequent figures, the curves were smoothed with a 1-yr running-mean filter). The maximum is found at ∼1000-m depth. The figure suggests that since around 1960 the MOC slightly increased in strength throughout most of the remaining GECCO integration period. At the end of 2001 an overall increase of about 2 Sv is reached, which mainly takes place during the late 1970s. The corresponding change in the unconstrained run (also shown in the figure) is larger (3 Sv) since a further increase takes place in the late 1980s. The overall variability pattern on shorter time scales, including the initial drop in the MOC strength, is quite similar in both runs. However, apart from the additional increase after the 1980s in the reference run, we note an oscillation on a decadal time scale in the transport differences imposed through the ocean data. The corrections to the MOC in the optimized run have been increasing since the mid-1980s. The correlations between the reference and optimized runs before and after 1987 are thus r = 0.80 and r = 0.59, respectively. Clearly the impact of the larger data density after 1992 is showing an effect, which was also noted by Köhl and Stammer (2008) who reported that the amount of data before mid-1970 is not sufficient to correct biases in the NCEP surface fluxes in the same optimization.
Also shown in the figure is the time series of the MOC strength obtained from the shorter 11-yr optimization described by Köhl et al. (2007b). The corresponding changes resemble the curve of the longer run on shorter time scales (the part that is driven by the atmospheric variability) superimposed with a slower adjustment process, which basically resembles the first 10 years of the reference run: an increase during the first 2 years, followed by a steady state over 4 years, and a decrease thereafter. As already discussed by Köhl et al. (2007b), the slow adjustment part of the MOC variability is related to a transition to a state that the model is able to maintain, whereby a large fraction of the variability during the first decade has to be considered as unrealistic. The same adjustment is also visible in Wunsch and Heimbach (2006) from a short optimization from 1992 through 2004. We therefore excluded the first decade form the analysis. During the later 4 decades, our analysis remains within the 6 Sv error bar of the measurements of Bryden et al. (2005). Only their 1957 estimate is an exception but falls within the disregarded first decade. This estimate is also crucial for their conclusion that the MOC is decreasing. Both estimates, therefore, remain consistent, although our simulated variability and positive trend does not particularly agree with their measurements.
5. Mechanisms for North Atlantic MOC changes
We have seen in the last section that the MOC in the North Atlantic is changing essentially on all time scales and that the latitudinal pattern of those changes is complex. While up to now there has been speculation about how those variations appear, this is also true regarding the question of what drives them. In this section, we will try to address the second aspect as it concerns the GECCO results.
Based on an adjoint sensitivity study with essentially the same setup as used here, Köhl (2005, hereafter K05) investigated mechanisms that, in principle, can affect the strength of the MOC at 30°N in the North Atlantic. In summary, these mechanisms comprise Ekman transport, coastal upwelling, Rossby, and Kelvin wave-based processes. Two of the four processes are buoyancy driven and involve propagation of density anomalies across the basin and along the basin boundaries and, therefore, involve longer time lags relative to the creation time. Density anomalies in the Labrador Sea Water (LSW) are communicated by internal Kelvin boundary waves southward, and density anomalies off the coast of Africa are communicated by long Rossby waves westward along the baroclinically unstable region of the subtropical gyre. The latter changes the MOC by modulating the western boundary current along the coast of Brazil.
The two others are wind driven and are generally of shorter time scale. On monthly to interannual time scales, the local Ekman transport caused by zonal wind stress is the most prominent mechanism, which is also described by Jayne and Marotzke (2001) and Böning and Herrmann (1994) to be the most important process responsible for the subannual variability. Coastal jets, caused by up/downwelling due to changes in alongshore wind at the coast of West Africa are found to account for a larger fraction of the wind-induced MOC anomaly in the perturbation experiments performed by K05. In the following, we will use these four mechanisms as a basis to identify how much each contributes to the MOC changes estimated by GECCO on interannual to decadal and longer time scales. Additionally, Lee and Marotzke (1998) suggested that wind stress determines, via the time-dependent barotropic vorticity relation, the external component of the MOC in the Indian Ocean. In the Atlantic, the external component was found by Sime et al. (2006) to explain more than 40% of the low frequency MOC variability at 24°N in their simulation. However, they also suggest that the wind stress curl may only explain a very small fraction of the external component in the Atlantic Ocean, whereby dynamics must be different from the Indian Ocean. The shallow northward transport follows as a western boundary current through the Straits of Florida. Since the straits are relatively shallow with minimum depths of around 800 m, the deep western boundary current is topographically forced to stay east of the Antilles. This zonal offset between the shallow northward and deep southward transport leads to a large barotropic component, which is not wind driven but topographically steered.
As a starting point, Fig. 5 shows anomalies (relative to the 1962–2001 mean) of the maximum MOC as a function of latitude and time in the Atlantic Ocean north of 35°S as they appear from the GECCO optimization (left panel). Anomalies are defined here as the deviations of annual means from the time mean MOC. The figure shows negative anomalies of around 1–2 Sv in the subpolar North Atlantic during the first 15 yr that subsequently propagate southward, reaching the South Atlantic after about 3 yr. A southward propagation of positive anomalies follows after 1975. The propagation seems to reflect a principle adjustment process of the ocean to high-latitude forcing variations. Those frequent anomalies on interannual period seem to be superimposed to a long-term trend, enhancing MOC amplitudes over the entire period north of 20°N. The propagation of southward moving MOC anomalies is slow north of 40°N and accelerates in the subtropics. Although the equatorial dynamics seems to be somewhat different, the signal subsequently seems to appear in the South Atlantic, albeit with some indication of a negative phase speed. A fraction of the variability in the tropical band (20°S–20°N) appears to be independent of the anomalies farther north.
The GECCO results have to be compared with the right-hand side of the Fig. 5, showing the MOC response to perturbations of the density field in the Labrador Sea and off the coast of Africa from K05. The panel summarizes how buoyancy-driven mechanisms involving boundary waves along the western boundary and Rossby waves crossing the basin in a westward direction jointly create MOC perturbations at 30°N. Consistent with the left panel, quite similar propagation patterns are obvious in the sensitivity study, which have slower phase speeds north of 40°N but accelerate in the subtropics.
We note also the signature of a northward propagating anomaly south of the equator. Köhl (2005) described two designated pathways for the unstable Rossby waves, namely the two baroclinically unstable regions near the edges of the subtropical gyre in both hemispheres. The pathway in the Southern Hemisphere was not discussed in detail in his paper, since the mechanism seemed to be the same as for the Northern Hemisphere, nor was coherence across the equator, visible in our Fig. 5b. Rossby waves impinging on the western boundary modulate the MOC by modifying the density at the western boundary. The associated signal subsequently propagates equatorward as Kelvin waves; for the Southern Hemisphere this propagation is northward. In the Northern Hemisphere, this is a southward propagation at the western side triggered either by Rossby waves or by Kelvin waves originating from the subpolar region. On the eastern side this will be a poleward propagation triggered by equatorial Kelvin waves: after Rossby waves in the Southern Hemisphere reach the western boundary, Kelvin waves are triggered, travel northward, and subsequently cross the equator leading to a northward propagation across the equator. The modification of the density along the eastern boundary by the Kelvin wave also changes the MOC. The corresponding process also exists for the Rossby waves in the Northern Hemisphere, leading then to southward propagation across the equator. In the left panel of Fig. 5, northward propagation across the equator is visible in years 1966, 1971, and 1978 but southward propagation across the equator is visible in 1975 and during later years. The direction of propagation allows determination whether the cause of MOC anomalies at or around the equator is located in the Southern or Northern Hemisphere.
As a result, a signal origination form the north (south) will lead to a southward (northward) propagation across the equator. After crossing along the equator, the Kelvin wave actually triggers two waves: the one already discussed that continues in the original direction and one traveling in the opposite direction toward the opposite pole. Since the associated density signal of the second Kelvin wave is of the same sign, its effect on east–west density gradients and, thus, MOC anomalies opposes the original signal on the western boundary and will therefore reduce (or even erase) the MOC signal. There is an indication for such a process in Fig. 5 (e.g., for years 1975 and 1978) as the anomaly disappears in the source hemisphere once the signal has crossed the equator.
The Kelvin/Rossby wave mechanism described above leads to a coherence across the equator, apparent in Fig. 5a, concentrated on the band 20°S–20°N and suggests that unstable Rossby waves contribute an important fraction to the MOC variability in the tropics—probably the dominant mechanism for the MOC variability there. Westward traveling Rossby waves were described before by Hirschi et al. (2007) to cause a significant fraction of the MOC variability on even the subannual time scale.
Most prominent in both panels of Fig. 5 is the propagation of density anomalies southward from the subpolar ocean on decadal time scale. Consistent with K05, we find a boundary wave phase speed of 15 cm s−1 between 30° and 40°N for the same model as ours. Getzlaff et al. (2005) discuss a resolution dependence of the boundary wave speed and they report a speed of up to 3.5 cm s−1 for a 4/3° model and 20 cm s−1 for a 1/3° model in the subtropics. The resolution dependence of models formulated on an Arakaw B grid was already described by Hsieh et al. (1983), who report no such dependence for C-grid models. Since the ECCO–MIT model is formulated on a C grid, we find wave speeds close to what Getzlaff et al. (2005) reported for a 1/3° model.
A lagged-correlation analysis provides a statistical tool to identify mechanisms leading to actually simulated MOC variability. In Fig. 6 we show therefore the lagged correlation of the maximum MOC at 25.5°N with the horizontal streamfunction at 850 m for different lags. Note that at 25.5°N the horizontal streamfunction at the western boundary is identical to the MOC in 850-m depth. For zero time lag, the largest correlation can be found along the western boundary, demonstrating a large meridional (roughly from 10°S to 40°N) coherence of the MOC variability. While some fraction of this area is located poleward, most of it reaches toward low latitudes, consistent with the large-propagation speed of Kelvin waves there. Additionally, a large correlation exists along 25°N, which is found to be concentrated in the upper 200 m (maximum correlation of r = 0.72 at 70 m) and to be related to the Ekman drift. The correlation pattern also exists for a 12-month lead but with a deeper maximum at ∼600 m. The negative MOC response after switching of the zonal wind speed forcing, as is it visible at about 37°N for a 12-month lag, was also noted by K05. It has a deeper maximum at ∼2500 m. A small pattern with high correlation exists off the coast of Africa, which is surface intensified and represents the coastal jet mechanism described by K05. Because of the narrow spatial location of these latter patterns, a possible spatial coherence of wind stress anomalies could not be the explanation for the coherence of the MOC changes across the equator on shorter time scales.
With increasing time lead, the maximum correlation pattern propagates northward along the western boundary, occupying the Labrador Sea (LS) after 5 yr, in agreement with the Kelvin/boundary wave mechanism described by K05 and the propagation shown in Fig. 5. For longer lead times, the high-correlation region crosses the subpolar North Atlantic and can be found along the Greenland–Scotland Ridges for a 120-month time lag, indicating an enhanced transport through Denmark Strait. The correlation with lags of 42 and 60 months also shows a propagation signal, now crossing the Atlantic zonally between 10° and 25°N, which corresponds to the Rossby wave mechanism described by K05. Similar to the findings of K05, we find the maximum correlation at 300-m depth. The corresponding process in the Southern Hemisphere leads to a respective signature farther south (at 25°S). However, due to lower correlation values (r = 0.43) and shallower depth (z = 100 m), it is not resolved in the figure. For various lead and lag times a correlation along the North Equatorial Countercurrent (NECC) is visible. The model shows large Rossby wave activity in this region on seasonal time scales. Stramma et al. (2005) suggested that Rossby waves, which may be generated in the NECC, may modulate the strength of the North Brazil Current. This pathway was, however, not visible in the perturbation experiments performed by K05, and it is not clear if a causal relation to the MOC at 25°N exists. To address possible impacts of wind stress curl on the external MOC component at 25°N via the barotropic vorticity equation, a separate correlation with the wind stress curl was calculated. The pattern (not shown) reveals coefficients smaller than 0.2 along 25°N, although larger correlation exists farther south and north of this latitude.
6. Causes for MOC variations in the GECCO estimation
The discussion in the previous section suggested that all four processes, described by K05 as potentially involved in changing the MOC at 25°N, actually contribute to forming the MOC variability at that latitude in the GECCO synthesis. Here we will now determine, in a quantitative way, the fraction by which each process explains the simulated-MOC variability and consider the variability associated with each of them. For this purpose, Fig. 7 (top panel) shows (top) the time series of the maximum MOC at 25°N, but now together with time series of the Ekman transport at 25°N, the transport due to Rossby waves, and the transport due to the coastal jet off Africa. Also shown is the maximum MOC at 48°N.
The Ekman transport was calculated from the estimated wind stress along 25°N, which is highly correlated with the transport in the upper 100 m (on which the correlation in Fig. 6 is based). The transport time series due to the coastal jet along Africa is based on meridional velocities in the upper 500 m at 19°N, east of 24°W, that, according to K05, covers the region of the response to the alongshore winds; that is, it is also wind driven. The lagged correlation in this region, shown in Fig. 6, varies from 0.2 to above 0.6, suggesting that other processes become more important farther west. Correlation values larger than 0.6 are only found very close to the coast in a region that accounts for only very little transport variability, which thus does not allow a realistic estimation of the variability.
Quantifying the role of density-driven variability, communicated through Kelvin and Rossby waves, requires the use of proxy transports. For the Kelvin/boundary wave mechanism we select the MOC transport at 48°N. At this location only a small fraction of the MOC variability is due to Ekman transports (r = 0.08), and, as shown by K05, the Rossby wave mechanism as a contributor to the MOC variability also becomes less important farther north. To determine the contribution due to Rossby waves amplified in the baroclinically unstable subtropical gyre, the horizontal streamfunction at 15°N, 45°W in 300 m (i.e., the 0–300-m transport east of this point) close to the maximum correlation for the 42-month lag in Fig. 6 was chosen. We note, however, that the decomposition of the contributions to the MOC at 25°N is problematic for the latter two parts as these transports include components of the other three mechanisms.
From the time series shown Fig. 7, it appears that the two wind-driven processes can only explain a small fraction of the variability of the MOC at 25°N, as their variability is smaller than the other two processes and occurs on shorter time scales. Accordingly, the increase in MOC during the late 1970s is only visible in the time series of the two-wave-based processes but not in those related to the wind-driven mechanisms.
Table 1 summarizes correlations and lags for the transport components. On time scales longer than 1 yr, the local Ekman transport at 25°N could explain (with r = 0.34) only a small fraction of the MOC variability at the same latitude. Maximum correlation is close to zero time lag. The variability of the coastal jet off Africa is slightly better correlated (r = 0.38) with the MOC at 25.5°N. In contrast, the relation of the MOC changes at 25°N with the Rossby wave signal is much closer on the interannual time scales (r = 0.71) and also shows a lag for maximum correlation at 38 months. It was already suggested above that the southward propagation of density anomalies originating from the subpolar North Atlantic is an important mechanism on interannual to decadal time scales. The figure reveals now that the signal at 48°N leads the signal at 25°N by about 3 yr with a large correlation of r = 0.70.
In a first attempt to quantify how much variability of the MOC at 25°N can be explained by the variability of each of the four components, one might be tempted to use the individual pairwise correlations (shown in parentheses) to estimate the explained variance. Both wave-based mechanisms individually would then explain about half of the variability, while the two wind-driven mechanisms together would explain about one-third. However, the existence of cross-correlations between pairs of the transport time series, which are also shown in Table 1, generally reduces the variance explained by the components. Accordingly, a multilinear regression that integrates information from all four processes is able to explain only slightly more than 70% of the MOC variance at 25°N.
These cross-correlations exist because all processes are mainly driven by different aspects of the same external surface forcing. The measures for the different processes may also insufficiently isolate the mechanisms. Moreover, a direct influence of one process on another also exists. For instance, as already mentioned above and demonstrated by K05, the MOC at 48°N is also modified through Rossby waves farther south or coastal up/downwelling modulating the coastal jet will create a fraction of the Rossby waves that propagate westward. However, the existing time lags of about 3 yr between the subtropical Rossby waves and MOC at 48°N or up/downwelling and the Rossby waves should remove most of this interrelation.
Since the time series are not independent, the pairwise correlation overestimates the linear dependence because of the influence of the remaining variables on both of the time series. This influence can be eliminated by regarding partial multiple correlation coefficients as they evaluate the specific proportion of variance explained by one independent variable only. Here, rx1, x2, x3 . . . xn is the correlation between x1 and x2 after the influence of x3 . . . xn is eliminated: The reduction formula by Yule (1907),
allows the recursive calculation of the reduced correlations from the cross-correlations, shown in Table 1, and thus the determination of the associated explained variance for each of the variables xi.
The results of the analysis are also summarized in Table 1. Now the two wave-based components together explain 60% and the wind-driven processes together only 13% of the MOC variance at 25°N. We note that the westward Rossby wave contribution, with 35% explained variance, is slightly larger than the contribution from the Kelvin wave signal.
A total of 27% of the variance remains unexplained. Candidates for the associated processes are the Rossby waves in the South Atlantic, the barotropic mode driven by wind stress curl, and the Rossby wave along the NECC. However, including the Rossby waves in the South Atlantic revealed that only little more (≤2%) of the variance is explained and that many more such secondary processes need to be considered.
To further explore the sources of variability for the MOC, we show in Fig. 8 the time series of the maximum MOC at 48°N together with a regression of the density from a region south of Denmark Strait and of the Denmark Strait Overflow with the MOC curve from 48°N, respectively. Also shown is the regression of the North Atlantic Oscillation (NAO) index on the 48°N MOC changes as part of Fig. 7. After an initial decline, the MOC at 48°N shows an increase after 1960 with a second major increase in the mid-1970s. Superimposed on this long-term trend is interdecadal variability. For the last decade we added the estimate from the previous 11-yr optimization and results from Lumpkin et al. (2008). Both describe, in agreement with the present estimate, a maximum MOC around 1996–97 and a decline thereafter. Lumpkin et al. also provide an error bar for their estimations (about 1–2 Sv) and explore the sensitivity to different assumptions (about 1–2 Sv). Since the effect of the adjustment process propagates southward, it is completed in the subpolar area within a shorter time span. Only the early years of the 11-yr optimization might still be affected. Apart from the long-term trends, a fraction of this interdecadal variability is represented by the three regressions. Regression coefficients between the 48°N MOC time series and the above regressions are summarized in Table 2. In addition, similar parameters are shown for the mixed layer depth and density in the LS, the density of the Irminger Sea, and the local Ekman transport at that latitude.
In contrast to Böning et al. (2006), we find a low correlation between the 48°N MOC strength and the maximum mixed layer depth in the LS (r = 0.3). Instead, and in agreement with K05, an enhanced correlation to the density at 900 m in the LS (r = 0.89) and an even higher correlation to the density at 900-m depth from south of the overflows (r = 0.91) is found. The latter is supported by the high correlation between the MOC at 48°N and the Denmark Strait Overflow variability (r = 0.83, with lead 4.8 yr). Together with the lower correlation to the mixed layer depth variability in the LS, this indicates that the correlation to the density variability in the LS is substantially caused by density variations that are created by the overflows and advected into the LS rather than created locally. Although different from the results of Böning et al. (2006), our relation to the overflows is not surprising as, according to Dickson and Brown (1994), more than 13 Sv, and thus the dominant part of the North Atlantic Deep Water, originates from the overflows. We refrain from determining the relative importance of the overflows in comparison to the LS convections, as these processes are not very well resolved in our model. We note, however, that Pickard and Spall (2007) concluded from their calculation of the diapycnal mass fluxes derived from hydrographic data in the LS that during the years 1990–97 the Atlantic MOC was not affected much by LS convection. Both mechanisms are consistent with a correlation of the NAO index on the MOC changes at 48°N. Köhl and Stammer (2004) have discussed the relation of the NAO changes to variations of total transport through Denmark Strait (DS). Biastoch et al. (2003) describe a relation of the barotropic transport between the subpolar North Atlantic and the Nordic Seas to wind stress curl changes associated with variations in the NAO, consistent with the DS transport measurements by Macrander et al. (2005). Further observational evidence is provided by Kieke and Rhein (2006), who found from bottle and CTD data increasing overflows from the 1950s to the 1980s and a decrease during the 1990s, which corresponds to our findings. They were, however, not able to establish a clear relation to the NAO. A mechanism of a NAO-related production of dense water north of the sills proposed by Käse (2006) explains a relation of the Denmark Strait Overflow variability, monitored by interface height measurements north of the sills, to the NAO.
Previously, Eden and Willebrand (2001) also investigated the response of MOC changes at 48°N to atmospheric forcing anomalies over the LS and reported an increase in MOC strength after a lag of 2–3 yr relative to enhanced convection activity during high NAO states. In the absence of overflow measurements, the NAO index appears to be a good proxy for changes of the MOC at 48°N and three years later at 25°N. On the other hand, Köhl et al. (2007a) have discussed the possibility of using SSH data north of the DS to monitor the Denmark Strait throughflow. SSH observations north of the DS sill should therefore also be a good proxy to monitor changes of creation of dense water in the subpolar North Atlantic and to serve as an early warning system, even for changes of the MOC at 25°N. We tested this hypothesis by regressing SSH north of the sill [at 25.5°W, 66.5°N as in Köhl et al. (2007a)] on the Denmark Strait Overflow and the MOC at 48°N. The former (r = 0.70), is considerably smaller than in the earlier study (r = 0.84), and the correlation with the MOC at 48°N thus decreased from 0.84 (with the original Denmark Strait Overflow time series) to 0.6. Reasons for the reduced skill of SSH as a proxy in comparison to Köhl et al. (2007a) probably include the lack of resolution in this narrow passage. Other model results also support a MOC increase, which can be understood in that a large fraction of the MOC variations can be described as a lagged response to NAO variability (Eden and Willebrand, 2001). The general upward trend of the NAO since the 1960s suggests, according to Latif et al. (2006), an strengthening rather than a weakening of the MOC.
Additional independent information about MOC changes was obtained from measurements of the subpolar deep western boundary current (DWBC) or the intensity of the gyre circulation. In this context, Häkkinen and Rhines (2004) reported an increase of the SSH in the subpolar North Atlantic over the past decade, which they interpret as a proxy for the decline of the subpolar gyre circulation. The interpretation was supported by direct current meter and hydrographic observations, which additionally point to local buoyancy forcing as a possible cause for the decline.
Böning et al. (2006) verified in high-resolution simulations the close correspondence between SSH changes and boundary current transports in the Labrador Sea, and ultimately the MOC strength farther south. Although wind stress was found to alter the gyre circulation, leading to a decline in the 1990s, variations of the deep fraction of the boundary current and the MOC were found to be primarily of thermohaline origin. This explanation is consistent with changes in the LS convection concurrent with the NAO described by Lazier et al. (2002). Our estimation confirms the decline and its relation to SSH in the center of the LS (not shown) reported by Böning et al. (2006), which are in agreement with the downward trend of the MOC at 48° after the mid-90s, shown in Fig. 8. However, measurements of the DWBC east of Grand Banks during 1999–2005 and 1993–95 show no indication for a slowdown of the DWBC and thus the MOC (Schott et al., 2006). This view is also supported by observations of the Deep Labrador Current (Dengler et al. 2006), which actually indicate a moderate increase from the period 1996–99 to the period after 1999.
The global estimate of the ocean state over the period 1992–2002 by Köhl et al. (2007b) was extended into the past to cover the 50-yr period 1952–2001. The extension required the solution of additional problems, demonstrating that extensions of state estimation efforts in integration periods is not trivial. The obtained solution is in much better agreement with the data than the reference run, but has not reached the same level of agreement as the previous shorter solution during the period 1992–2001. The question remains why the mean MOC is lower than the reference run, as is consistently found to be true for the other two shorter optimizations. Long-term state estimation efforts are necessary to investigate changes in the circulation, especially the MOC, since the adjustment to the overflows dominates processes associated with the MOC variability during the first decade of the simulation. Over the last four decades we find a general increase of the MOC at 25°N by about 2 Sv. Although this disagrees with Bryden et al. (2005), who reported a reduction by 30% over the last 50 years, our estimate remains within error bars and thus consistent with their estimates. The reason for this seeming paradox is that the 1957 estimate is necessary for the significance of their reduction but is excluded from the comparison over the shorter period. At 48°N, the MOC shows, in agreement with Lumpkin et al. (2008), a maximum around 1996–97 and a decline thereafter.
Because of the complexity of our estimation problem, a formal error of our MOC estimates could not be achieved (although formally there is no obstacle to such an estimate), and a most promising evaluation of our results is likely a comparison with estimates available from ongoing RAPID measurements. For such an evaluation of the estimated MOC time series, the estimation has to be extended until the present day so as to establish an overlap with estimations from the ongoing RAPID estimates. Moreover, because the picture emerging from the published work regarding ongoing MOC changes is not entirely consistent, to date no final answer can be given, and our work is meant to provide an extra piece to the puzzle. Nevertheless, our analysis provides new insight about the processes that potentially influence the MOC variability at 25°N. A correlation analysis reveals, in accordance with K05, that four different processes contribute to the MOC variability there: namely, the two wind-driven processes that do not involve lags—the Ekman drift and the coastal jet off Africa—and two remote processes that are communicated by waves. In agreement with Sime et al. (2006), we find that the Ekman transport accounts for less than 20% of the interannual (pentadal) MOC variability at 24°N. The wave processes are based on Rossby waves in the baroclinically unstable region of the subtropical gyre and Kelvin waves that communicate density anomalies from the subpolar region toward the subtropical region. Both wave processes involve time lags. Accordingly, the variability of the MOC at 48°N leads the variability at 25°N by 2–3 yr and can be understood as responses to atmospheric forcing characterized by NAO. Among the mechanisms that create MOC variability at 25°N, the Rossby wave-based mechanism, with 35% explained variance, is the dominant component, followed by the southward propagation of the subpolar density variability, which is mainly caused by Denmark Strait Overflow anomalies and explains one-third of the MOC variability.
We thank all members of the ECCO consortium for their direct or indirect support and discussions that led to the results being presented in this paper as well as the anonymous reviewers for their comments. In particular, Dietmar Domenget, Charmaine King, Diana Spiegel, and Kyozo Ueyoshi helped with the processing of the input data. Reanalysis surface forcing fields from the National Center for Environmental Prediction–National Center for Atmospheric Research (NCEP–NCAR) are obtained through a computational grant at NCAR. Computational support from the National Partnership for Computational Infrastructure (NPACI) and the National Center for Atmospheric Research is acknowledged. Supported in part through ONR (NOPP) ECCO Grants N00014-99-1-1049 and N00014-99-1-1050, NASA Grants NAG5-11765, NAG5-12870, and NNG04GF30G, a contract with the Jet Propulsion Laboratory (1205624) and a grant from the Bundesministerium für Bildung und Forschung as part of the joint project NORDATLANTIK. This is a contribution of the Consortium for Estimating the Circulation and Climate of the Ocean (ECCO) funded by the National Oceanographic Partnership Program.
Corresponding author address: Armin Köhl, Institut für Meereskunde, Zentrum für Meeres- und Klimaforschung, Universität Hamburg, Bundesstrasse 53, 20146 Hamburg, Germany. Email: firstname.lastname@example.org