1 Introduction

Central to the detection of climate change is the estimation of separate signal and noise patterns (e.g., Hegerl et al. 1997; Barnett et al. 2001; Drijfhout and Hazeleger 2007). The signal comprises the signature of climate change induced by increased anthropogenic greenhouse gas concentrations, while the noise is the climate variability induced by natural processes. In most detection studies it is implicitly assumed that the noise is stationary and can be estimated from a control simulation or historical database. This assumption does not rely on a fundamental property of the noise, but is born from the poor knowledge of how the noise (internal variability) reacts to a specific forcing. Recently, however, a few studies have appeared in which the change in natural climate variability was addressed that is due to anthropogenic forcing. In some studies tropical climate variability (Collins et al. 2005; Hoerling et al. 2007) features a more distinctive climate change signal than extratropical variability (Rauthe et al. 2004; Stephenson et al. 2006). A similar study of the anthropogenic change in ocean variability has not yet appeared.

One type of ocean variability that is likely to be sensitive to anthropogenic forcing is the internal variability of the Atlantic Meridional Overturning Circulation (MOC). The MOC itself displays a marked sensitivity to anthropogenic climate change (e.g., Schmittner et al. 2005). The forced change in the MOC is related to a different spatial structure of deepwater formation (Wood et al. 1999; Drijfhout and Hazeleger 2006). It seems very likely that such alterations in convective activity in the North Atlantic cause changes in the internal variability of the MOC. Also, a change in the atmospheric variability can be a reason for MOC-variability to alter. For instance, a possible change in the North Atlantic Oscillation (NAO) might cause the internal variability of the MOC to react (Dickson et al. 1996; Hurrell et al. 2006). The MOC is associated with a large northward heat flux (Ganachaud and Wunsch 2003). Variations in the MOC are related to ocean heat transport variations, which on their turn are coupled to atmospheric heat transport variations through the mechanism of Bjerknes compensation (Shaffrey and Sutton 2006; van der Swaluw et al. 2007). As a result, changes in MOC-variability can have a direct impact on climate variability over Northwestern Europe.

In this paper, we address the change in MOC-variability that is due to anthropogenic forcing by analyzing a large (62) member ensemble of climate model simulations, all forced with increasing levels of greenhouse gas concentrations in the atmosphere. With such a large ensemble, the forced signal and internal variability can be separated from each other and described more robustly than is possible with only a few runs. We simply assume that the forced signal is contained in the ensemble-mean average. Then, the internal variability in each individual ensemble member is defined as the anomaly with respect to the ensemble-mean signal. This definition allows for a complete description of the temporal and spatial character of both the forced signal and the natural internal variability. In the model used here, the MOC is too strong. However, both pattern and amplitude of the forced change, as well as the pattern and amplitude of the interannual variability compare well with other model estimates, see also Drijfhout and Hazeleger (2007).

2 The ensemble experiment

The model used for this study is version 1.4 of the Community Climate System Model (CCSM) of the National Center for Atmospheric Research. It simulates the evolution of the coupled atmosphere-ocean-sea-ice-land system under prescribed climate forcing. The atmosphere was run with a spectral resolution of T31 and 18 levels in the vertical. The land model distinguishes between specified vegetation types and contains a comprehensive treatment of surface processes. The sea-ice model includes ice thermodynamics and dynamics. The ocean model has 25 vertical levels and a 3.6° longitudinal resolution. The latitudinal resolution ranges from 0.9° in the tropics to 1.8° at higher latitudes. The coupled system does not require artificial flux corrections to simulate a realistic climate (Boville et al. 2001).

The simulations cover the period 1940–2080. Until 2000, the forcing includes specified estimates of temporally evolving solar radiation, temporally and geographically dependent aerosols due to manmade and natural emissions and time dependent greenhouse gases (GHGs). From 2000 onwards all types of external forcing are kept at the 2,000 values, except for the concentration of GHGs which increase according to a “business-as-usual” (BAU) scenario. This BAU scenario closely follows the Special Report on Emissions Scenarios (SRES) A1 scenario and CO2 concentrations increase to 640 ppmv by 2080, the end of the simulation period (Dai et al. 2001). For more details on the model and on the ensemble experiment we refer to Boville et al. (2001) and Selten et al. (2004), respectively.

An ensemble of 62 simulations was produced, each covering a 140-year period. The simulations differ only in a small random perturbation to the initial atmospheric temperature field, which leads, due to the chaotic nature of the atmospheric flow, to completely different weather patterns by the end of the first month. After a few years of spin up, the ocean circulation in each ensemble member has “forgotten” its initial state. Averaging over all ensemble members filters out most of the internally generated climate variations and allows for an estimate of the externally forced climate signal. Also, by analyzing the detrended time series of all ensemble members the internal variability can be much more adequately described than is possible for a single climate-scenario run.

In the following, we will disregard the first 10 years of each run, to allow the ocean circulation in each ensemble-member to get decoupled from its common initial state. For this reason, our analysis will start in 1950 instead of 1940. Figure 1 illustrates the ensemble experiment by showing the ensemble-mean, yearly-averaged overturning in the Atlantic at 35°N and 1,000 m depth (where the MOC-variability is maximal) as a function of time (thick black line), showing the response of the MOC to GHG forcing since 1950, together with the 62 ensemble members (purple dots) showing the spread due to the yearly-averaged internal variability. We discriminate two time periods; the past climate, running from 1955 to 1995; and the future climate, running from 2035 to 2075. The length of these periods (40 years) has been chosen to resolve decadal variations. In the following we will compare these two time periods and we refer to future change when we discuss the difference between the future climate (2035–2075) and the past climate (1955–1995).

Fig. 1
figure 1

The ensemble-averaged overturning at 35°N and 1,000 m depth as a function of time (thick black line) showing the response of the MOC to GHG forcing since 1950, together with the 62 ensemble members (purple dots). The green line shows the overturning for one arbitrary member. Units are in Sverdrups (=106 m3s−1). The two time periods discussed throughout this paper; past climate (1955–1995) and future climate (2035–2075) are denoted by a blue line and red line

3 The anthropogenic change in the MOC

3.1 The pattern of change

From Fig. 1 we see that in the early phase of the integrations the maximum strength of the MOC. The maximum occurs at 35°N and 1,000 m depth. Initially, it equilibrates at a peak value of about 33 Sverdrup (Sv) (1 Sv = 10 m6 s−1). This level of overturning is maintained until the year 2000. The following 80 years (2000–2080) show an almost linear decrease of about 5–6 Sv. In addition, the figure shows that the internal variability decreases: the probability density function for the MOC sharpens.

The Atlantic MOC is dominated by a single, basin-wide overturning cell. With a peak value of 33 Sv at 35°N, the MOC is definitely too strong compared to observational estimates (Ganachaud and Wunsch 2003). The too strong MOC results from model drift after coupling (Frank Bryan, personal communication). It is a particular feature of the coupled CCSM1.4 model. Although there appears to be a bias in the ensemble-mean overturning circulation, the response to enhanced greenhouse forcing compares well with estimates from other models. The climate sensitivity of the MOC is relatively low in the model, but the absolute value of the trend compares well to other model estimates (Schmittner et al. 2005), and also the internal variability compares well to other models, see Drijfhout and Hazeleger (2007) for more discussion.

The future change in the MOC is shown in Fig. 2, together with an average over the past years. The pattern of change shows an overall decrease of the MOC that extends over the full latitudinal range of the basin. The largest decrease is found at about 2,000 m depth, well below the level where the MOC itself attains its maximum. This implies that the MOC not only weakens, but also shoals. Apart from an almost uniform weakening of the MOC, Fig. 2 shows a conspicuous strengthening near 55°N, between 0 and 1,000 m depth. This indicates that there is an overall decrease in deepwater formation, but, as we will show later in more detail, at a few places deepwater formation increases or may establish itself for the first time. A discussion of the inherent change in North Atlantic DeepWater can be found in Drijfhout and Hazeleger (2006).

Fig. 2
figure 2

Future change in meridional overturning circulation (Sv in 80 years), overlaid by contours for mean overturning in the past climate (1955–1995) in Sverdrups. Future change is the average for the period 2035–2075 minus the average for the period 1955–1995

The MOC has been found to vary from subannual to interdecadal timescales (Bryden et al. 2005; Polyakov et al. 2005; Meinen et al. 2006). On shorter timescales MOC-variability is dominated by the wind-driven changes in Ekman transport and the barotropic return flow established beneath. On intermonthly to interannual timescales, changes in Ekman pumping may induce a Kelvin and Rossby-wave response (Hirschi et al. 2007). The associated baroclinic adjustment to Ekman transport changes transfers the overturning point from the base of the Ekman layer (10 m) associated with the fast barotropic response, downward to the depth where the first barolinic mode changes sign (1,000 m). On even longer timescales MOC-variability is linked to larger-scale density changes which may be associated with atmospheric variability through changes in heat flux and (deep) convection, or MOC-variations may be due to a self-sustained or damped oscillation which is excited by atmospheric noise (Delworth and Greatbatch 2000).

Figure 3a, b shows the MOC-variance for interannual and for interdecadal timescales, respectively. In this figure the variance that is due to an apparent numerical instability in the deep tropical regions (Gent, personal communication) has been removed. This variance increases with time and peaks at 3,500 m depth at the equator with a standard deviation of 4.5 Sv. To remove this variance we sequentially performed a regression of the 2D MOC time series on the time series of each suspicious point, which showed signs of numerical instability. The associated variability was removed by subtracting the time series of these regression patterns from the 2D MOC time series. The decrease of the variance due to this operation is further away from the tropics only less than 5%. Figure 3 shows the variance of the remaining variability. Interannual variability now attains maximum values at 35°N and 1,000 m depth, near the overturning maximum (see Fig. 2). We calculated the interannual variability from the variance of the annually averaged MOC time series. The amplitude of the variability as measured by the standard deviation is about 2 Sv. When subannual timescale are included the standard deviation exceeds 4 Sv at mid-latitudes and even 8 Sv in the tropics (Drijfhout and Hazeleger 2007). The future change in variance consists of an overall decrease of the variance that almost scales with the variance itself (Fig. 3a). The exception is the region at 55°N where the variance increases in the upper layers, associated with a local strengthening in overturning (Fig. 2).

Fig. 3
figure 3

The upper panel (a) shows the future change in variance at all interannual timescales overlaid by contours for the variance in the past climate. Units are in Sv2. The lower panel (b) shows the same for decadal variance only (11 year low-pass filtered)

At interdecadal timescales, obtained by smoothing the annual mean MOC time series with an 11-year Welch time filter, the standard deviation becomes less than 1 Sv (Fig. 3b). Maximum values initially occur at the same latitude as for the interannual variance (35°N), but now at larger depths: 1,500 m. Between 20°N and 40°N the variability reduces, but around 50°N there is an increase between 1,000 and 3,000 m depth. The upper layer-increase at 55°N is even more prominent at interdecadal timescales than at interannual timescales. The absolute amplitude is almost similar, indicating that subdecadal timescales hardly contribute to this feature. Eventually it starts to dominate the pattern of the interdecadal variance. The pattern of interdecadal change is more complex (multi-signed) than the pattern of interannual change, which shows a decrease almost everywhere. This suggests that it is (partly) related to changes in ocean dynamics. In the following we will focus on interdecadal changes. To further investigate this change in interdecadal variance we show the regression pattern of the time series of anomalous overturning on the time series of overturning at 35°N and 1,000 m depth, where the variance of the MOC is maximal. The regression pattern (i.e., the variability contours) we obtain (after applying an 11-year Welsh filter) and which is displayed in Fig. 4a, strongly resembles the first EOF. This confirms that both the regression analysis and the EOF decomposition recover the dominant mode of variability. Figure 4a shows a dipole pattern that is associated with the future change of this mode. Also the pattern of the change in variance (Fig. 3b) features a dipole with negative values between 20°N and 40°N and positive values between 40° and 60°N, but details are different, especially the depth where changes are maximal. Also the feature at 55°N, which is prominent in the pattern of the change in variance, is absent in the pattern of change for the dominant mode. The reason for these differences is that in Fig. 3b the displayed changes, in particular the alteration at 55°N, are partly associated with more local modes of variability which explain less variance than the dominant mode of variability shown in Fig. 4a. The dipole pattern of the future change in the dominant mode of variability can be explained as a poleward shift of this mode. In the remainder of this paper we attempt to elucidate the mechanisms that induce this poleward shift. We hypothesize that it is associated with increased overturning and variance at 55°N, although the changes at 55°N themselves do not show up in this pattern (they do show up with time lags of 5 year or more, but then the regression values become very low everywhere).

Fig. 4 °
figure 4

The upper panel (a) shows the future change in the pattern of the dominant interdecadal variability of the MOC, overlaid by contours of this pattern for the past climate. This pattern is obtained by regressing time series of the MOC on overturning time series at 35

N and 1,000 m depth. All time series were smoothed with an 11-year Welsh filter. Units are dimensionless. The lower panel (b) shows the change in pattern of MOC-variability that is associated with the NAO, obtained by a regression of the MOC on the NAO-index (index of Principal Component time series of the first EOF in pressure). Units are mSv hPa−1

3.2 The relation with the NAO

The MOC is thought to display an integrated response to the NAO (Eden and Willebrand 2001), implying that the regression between maximum overturning and NAO-index might peak at a certain lag. In the present study the regression peaks with the NAO leading with 1 year. This suggests that in the present model the response of the overturning to the NAO is dominated by the fast, wind-induced variations that are associated with the NAO; a blend of the almost instantaneous response to Ekman transport variations and the delayed response of internal Kelvin and Rossby waves to changes in Ekman pumping. The correlation between the MOC and the NAO, however, is not very strong. Time-filtering hardly increases the (lagged) correlation. The maximum value is 0.45. Despite the low correlation, a regression of the anomalous overturning on the NAO almost completely recovers the first EOF of MOC-variability, and the pattern displayed in Fig. 4a, b. In each case we obtain a monopole centered in the region where the overturning and its variance are maximal. But the NAO-related variability peaks 5° further north and we also find a smaller, separate monopole at 55°N. The pattern of change in the NAO-related variability, on the other hand, is different from the pattern of change shown in Fig. 4a. The response of the MOC to the NAO diminishes with time. In the present model the lagged response of the MOC to NAO-variations is quite weak for time lags of 5 years or more. Apparently, NAO variations project quite weakly on convective changes in the model. Even at decadal timescales the MOC response to the NAO is dominated by a baroclinic wave adjustment to changes in Ekman transport/pumping that occur preferably around 40°N in the storm track region. The expression of baroclinic adjustment in terms of meridional overturning circulation decreases due to the overall in crease in stratification of the ocean. This is arguably the reason why the MOC becomes less responsive to the wind-stress variability of the NAO. The poleward shift of the dominant mode of variability cannot be explained by changes in the relation between MOC and NAO. Also the NAO does not change in response to anthropogenic forcing in the present model.

Figure 5 shows the autocorrelation functions and spectra of the pattern shown in Fig. 4a for two time periods. The spectra (Fig. 5b) are red with two bands of increased variance; one near 4–8 years and one near 10–30 years. An inspection of the autocorrelation of this mode (Fig. 5a) reveals that it can be interpreted as a strongly damped mode with a pulse-like excitation phase. This mode cannot clearly be associated with an oscillation; there is no preferred timescale except for its damping timescale. In the past climate, the first crossing of zero autocorrelation occurs at 2 years which denotes one quarter of the cycle; the second crossing at 15 years. Especially the timescale associated with the second crossing of the zero autocorrelation line, which defines the overall timescale of the mode, increases with time, as is manifest from the change in spectrum shown in Fig. 5b. In the future climate the longer timescales become relatively more important; the mode becomes more persistent after the initial excitation phase. The stronger presence of longer timescales increases the time at which the zero crossing occurs, although the timescale associated with the excitation of the mode does not change; the width of the peak at zero lag does not widen (Fig. 5a). The shoulders between 2 and 8 years that arise in autocorrelation function in the future climate must be due to the stronger presence of longer timescales. On the basis of Figs. 4, 5 we now construct our central hypothesis. The poleward shift of the dominant mode of MOC-variability is due to a more important role for subpolar heat flux and convection changes in exciting this mode. In the following sections we will test this hypothesis.

Fig. 5
figure 5

The upper panel (a) shows the autocorrelation function of the dominant mode of the MOC for the past (blue) and future climate (red). The lower panel (b) shows the associated power spectrum in Sv2, together with the 95% confidence level (in black) for a red-noise spectrum

4 Changes in heat flux and convection

To test our central hypothesis we need to evaluate the GHG-forced changes in convection, which are possibly related to the poleward shift of the sea-ice margin that we expect to occur in a warming world. We will also look at net surface heat flux changes in relation to convection, as convection itself is not an easy observable, whereas the net surface heat flux is. First we investigate the changes in outgoing heat flux. Figure 6 shows a reduction of the outgoing heat flux in the Labrador Sea and central subpolar gyre. This suggests that convection diminishes there in future climate. At the same time there is an enhancement of the outgoing heat flux in the Baffin Bay, Irminger Sea and Greenland Sea. These enhanced outgoing heat fluxes indeed are related to the future change in sea-ice coverage, as depicted by Fig. 7a. In all regions where the outgoing heat flux increases, sea-ice coverage diminishes.

Fig. 6
figure 6

Future change in the outgoing heat flux, overlaid by contours for the heat flux in the past climate. Units are Wm−2

Fig. 7
figure 7

The upper panel (a) shows the future change in fraction of sea-ice coverage, overlaid by contours for sea-ice coverage in the past climate (dimensionless). The lower panel (b) shows the same for maximum (winter) mixed-layer depth in m, defined as the depth at which the surface to depth density difference exceeds 0.125 sigma units

Retreat of the sea-ice may lead to the opening of new convective areas. Because the calculation of mixed-layer depth is not a standard output diagnostic, the maximum mixed-layer depth was calculated for only 1 ensemble member. This member was selected for being the most representative of all 62 members with respect to heat flux changes, compared to the ensemble mean change in heat flux. Figure 7b confirms that indeed new convection sites appear southeast of Greenland and in the Baffin Bay as indicated by deeper maximum mixed-layer depths. In these areas the maximum surface density in late winter, when the mixed-layer depth is largest, increases (not shown), indicating that mixing to greater depth occurs, while everywhere else the maximum density decreases due to global warming. A particularly strong decrease of mixed-layer depth occurs in the Labrador Sea and in the Norwegian Sea where convection is reduced.

The changes in ocean mean state discussed above are caused by the GHG-forced changes in deepwater formation. But to study their impact on MOC-variability we also need to look at the changes in variability. Figure 8a shows the change in heat flux variance. The surface heat flux variability becomes stronger in the Labrador Sea and western Irminger Sea, where the mean outgoing heat flux diminishes. Vice-versa the variability weakens in the Baffin Bay and southeast of Greenland, where the mean outgoing heat flux increases. This can be understood by considering the associated changes in convective activity. In the Labrador Sea and western Irminger Sea the increased variance is due to an increased stability in these regions: convection becomes more irregular and less frequent. In general, when a transition occurs from convective to non-convective activity in a certain region, a regime is encountered where the flow is bi-stable with respect to convective instabilities. Both the on and off-states are possible solutions and hysteresis behavior can be found: periods with convection on and off tend to be prolonged (Lenderink and Haarsma 1994). This is what happens in the Labrador Sea and western Irminger Sea due to global warming. In the Baffin Bay and southeast of Greenland the variance decreases; convection becomes more frequent in these regions. But more important, heat flux changes associated with sea-ice variations lessen there.

Fig. 8
figure 8

The upper panel (a) shows the future change in (square root of) heat flux variance, overlaid by contours for the variance in the past climate (units are Wm−2). The lower panel (b) shows the same for the square root of the variance in maximum winter mixed-layer depth in m

This is substantiated in Fig. 8b. It shows that the increase in heat flux variance in the Labrador Sea and western Irminger Sea indeed is related to a growth of the variability in deepwater formation. Both the variance in maximum winter mixed-layer depth and maximum winter surface density (not shown) increase in these two regions in the future, while the average maximum winter mixed-layer depth and maximum winter surface density diminish. Deepwater formation is reduced and becomes more irregular, so the deepwater formation becomes more variable. Especially in the Labrador Sea the variability in deepwater formation on interannual timescales does not increase, while the variability on decadal timescales does. This is caused by the switches between off and on states, which occur on decadal timescales.

In the new convective regions; the Baffin Bay and southeast of Greenland, the variability in maximum winter mixed-layer depth is not enhanced but diminishes. On the other hand, the mean maximum winter mixed-layer depth increases in these regions. The variance in convection diminishes, but the amount of (deep) water formation grows. In the next section we will relate these changes to changes in the MOC.

5 The relation between heat flux changes and MOC-variations

Initially, the surface heat fluxes in the Irminger Sea, Labrador Sea and to a lesser extent in the GIN Sea are strongly related to the overturning at 35°N and 1,000 m depth. Heat flux variations southeast of Greenland and in the Baffin Bay initially have a negligible or even negative impact on MOC-variations.

Figure 9 shows the future change in the regression pattern of surface heat flux on the overturning time series at 35°N. The response time between the convection areas and 35°N consists of various timescales. Some of these timescales are longer than the decadal timescale that results from averaging with an 11-year Welch time filter. We have not considered lagged-regressions with time lags longer than 5 years, however, because lagged-regression values become quite low at these lags and they do not alter the picture that is implied by Fig. 9. The change in regression consists of increased values almost everywhere in the Northern Seas: larger regression values can be found in the Baffin Bay and Labrador Sea, the central subpolar gyre and western Irminger Sea and Greenland Sea. The interpretation of this almost uniform increase over the subpolar gyre region is twofold. First, in the new convective areas heat flux variance initially was related to sea-ice variations but heat flux variations become related to deepwater formation variability in the future. This interpretation is corroborated by the strong resemblance of regression and correlation patterns between heat flux and overturning time series (not shown): increased regression is always associated with increased correlation and vice-versa. Although the variability of surface heat fluxes reduces, the relation with the MOC at 35°N becomes stronger. Second, in the old convective regions deepwater formation becomes more irregular and heat flux variations grow there. As a result, the regression of these variations on the overturning at 35°N is also enhanced.

Fig. 9 °
figure 9

Future change in the relation between outgoing heat flux and overturning at 35°N, overlaid by contours of this pattern for the past climate. This pattern is obtained by regressing time series of the outgoing heat flux on overturning time series at 35

N and 1,000 m depth. Units are Wm−2 Sv−1. The pattern of hatched lines denotes the area in which the heat flux is integrated to calculate the regression pattern of Fig. 10

In particular in the Baffin Bay heat flux variations initially are weakly anti-correlated with MOC-variations. The variance in heat flux is large, but it is not associated with the MOC. Sea-ice variations dominate the heat flux variance. But as time proceeds, heat flux variance diminishes in this area, because sea-ice almost completely melts. The relation between heat flux variations and MOC-variations becomes stronger, but the impact of convection in this area on the MOC remains weak. We can disregard the role of this area in inducing the changes in the MOC that cause the poleward shift of its dominant mode of variability.

In the Labrador Sea and central subpolar gyre heat flux variations are initially positively correlated with MOC-variations. This is explained by the important role of these areas in forming deepwater and contributing to the deep return branch of the MOC. In a warmer climate these regions become less unstable as indicated by reduced winter mixed layer depths. Deepwater formation becomes more irregular and heat flux variations grow. Also the regression of heat flux variations on MOC-variability increases. In these regions a transition from a convective to a bi-stable state with convection switching between on and off modes causes the MOC-variability to be stronger related to heat flux variations caused by variations in deepwater formation.

In the western Irminger Sea, and the region southeast of Greenland, the heat flux variance is initially low and a relation between heat flux variations and MOC-variations is absent or weak. Most heat flux variations are related to sea-ice variability. Both heat flux variance and the relation between heat flux variance and MOC- variations grow with time. The increase of heat flux variability is due to the onset of deep convection. This region becomes a new convective area, but convection is still irregular.

In the Greenland Sea the heat flux variance initially peaks, especially northeast of Iceland. But these heat flux variations are only weakly connected to MOC-variations. A larger part of the variability is associated with sea-ice variability. We see a decrease of the variance here, the role of sea-ice variability dramatically lessens. At the same time, the relation between MOC-variations and heat flux variations becomes more important. Eventually, deepwater formation becomes established and it becomes less irregular than in the Irminger Sea. In the Norwegian Sea convection diminishes and heat flux variability and deepwater formation variability also decrease. There is no enhancement of the relation between heat flux variance and MOC-variability in this region.

So, can the hypothesis be confirmed that changes in deepwater formation drive a poleward shift of the dominant pattern of MOC-variability? To this end we compute the regression of the MOC time series on a heat flux time series averaged over a region in the Atlantic subpolar gyre. The region over which the heat flux is averaged is shown in Fig. 9. The change in the regression pattern is shown in Fig. 10. There are strong similarities with the change in the pattern of dominant MOC-variability shown in Fig. 4a. Figure 10 shows the change of all types of MOC-variability in response to a specific heat flux forcing, while Fig. 4a shows the change of one (dominant) mode of MOC-variability in response to all possible forcings. We cannot expect these figures to match completely, but the dominant dipole feature that describes the poleward shift of this variability is clearly present in both figures. We note that in Fig. 10 a strong change occurs near 55°N, 500 m depth which is absent in Fig. 4a. This feature describes the increased coupling between a local mode of variability at 55°N and the dominant mode, centered at 35°N. By and large, we may conclude that our hypothesis is indeed confirmed by Fig. 10. MOC-variability becomes stronger related to subpolar heat flux variations which result from variations in convection. This strengthened relation is associated with a poleward shift of the areas where convection occurs. But most important is that both the opening of new convective sites and decreasing deepwater formation in existing convective sites leads to an increase of the areas which are neither permanently (on interannual timescales) convective or non-convective, but which are subject to irregular switches between on and off states of convection. This makes the deepwater formation in the Nordic Seas more variable in time. As a result, the dominant pattern of internal variability shifts poleward.

Fig. 10
figure 10

Future change in the relation between the subpolar Atlantic heat flux and the MOC, overlaid by contours of this pattern for the past climate. This pattern is obtained by regressing time series of MOC on time series of the heat flux, averaged over the subpolar gyre region indicated in Fig. 9. Units are mSv W−1m2

6 Discussion and conclusions

The model results that we presented contain changes in MOC-variability that are driven by changes in convection due to anthropogenic warming. The record of detailed ocean observations is too short to validate/falsify these model results. Therefore, the main question we have to pose is whether the model results we have discussed are plausible and likely enough to be relevant to other models as well as the real world, or whether the results are model-dependent and non-representative.

The model suggests that in a gradually warming world convective areas may be subject to a particular life-cycle. Initially, in areas covered by sea-ice convection is inhibited. When the ice starts melting, the outgoing heat flux increases and periodically convection may occur. While more ice melts, convection increases and becomes well-established and less irregular. When heating further proceeds the water column starts to stratify again. First, in late winter the column will be subject to a transition from an unstable to a bi-stable regime, which is characterized by irregular switches between convective and non-convective periods. Later on, the column will become permanently stratified. In the transition between on and off-states of deep convection each area features periods in which convective events are irregular. It is in these periods that the heat flux variability on longer (decadal) timescales will attain maximum values.

The physical mechanism for such a life-cycle is general (Lenderink and Haarsma 1994) and not model-dependent. In the model that is used here the MOC is too strong. Too much deepwater formation occurs and some areas of deepwater formation (e.g., central subpolar gyre) are questionable. Also, the model seems to feature a numerical instability in the (deep) tropical ocean, which is a spurious source of mixing there. But these factors seem not to be relevant to the conclusion that future warming will be associated with increased irregularity in convective events in those areas where at present deepwater formation occurs on a more permanent basis, as the increased warming acts to stratify the water column in those areas. Another illustration of the mechanism of increased irregularity in deep convection due to anthropogenic warming is given by Schaeffer et al. (2002). These authors conclude that the inherent bi-stability of convective activity severely limits the predictability of the mostly abrupt transition from the on state to the off state of convection.

The non-linearity in changes in convection and its impact on the MOC is particularly well illustrated by the overturning strength at 55°N. The future change in overturning shows an increase for this region, as defined by the difference of the mean values for the period 2035–2075 and 1955–1995. However, if we inspect the overturning time series at this location (400 m depth; not shown), we find in particular that the overturning increases till 2045, after which the overturning decreases again. The variability peaks in the period 2020–2060 and reduces thereafter, see Fig. 11. In this period, most ensemble members feature an abrupt transition to larger overturning values, followed by a similar abrupt transition to lower values of overturning. The timing of these transitions, however, differs for each ensemble member.

Fig. 11
figure 11

The change in 19 year running mean standard deviation with respect to 1950 as a function of time and latitude at a depth of 364.2 m

In particular, the model suggests that the subpolar gyre as a whole may undergo a transition towards a larger irregularity in the occurrence of convection when the overall convection decreases. At present deep convection occurs at various places in the North Atlantic and it is documented to be subject to a certain degree of irregularity, in particular the episodic cessation of convection in, for instance, the Labrador Sea (Dickson et al. 1996). It is likely that the overall irregularity in deepwater formation will grow when the northern waters increase their stratification due to global warming.

It cannot be ruled out that in response to further warming convection might reestablish a more steady way of operation than in the near future, but occurring on a smaller scale than at present. This might happen when all “old” convective sites have been replaced by newer ones. When this indeed would occur, gradual warming, or cooling, would be associated with periodic transitions on a time-scale of hundred years, or more, from a more steady to a more irregular mode of convection. In that case, deepwater formation would be either characterized by an ensemble of regions where on or off-states dominate, or an ensemble of regions where in each region convection irregularly switches between on and off. The model suggests that if this would indeed occur, superimposed on a long-term trend of a poleward (equatorward) shift of the MOC and its variability that is associated with gradual warming (cooling), a centennial oscillation takes place in which the dominant MOC-variability is meridionally displaced. It would be characterized by an enhanced relation to northern heat flux variations in periods of increased convective variability. It would be interesting to test this hypothesis in a much longer climate model integration in which the climate system is imposed to a gradual anthropogenic warming. Moreover, such behavior might also be relevant to transition periods between ice ages and interstadials.

A second prediction from the model is that the timescale of MOC-variability increases when the relation between MOC and subpolar heat flux variations becomes stronger. This is due to the fact that in the model the relation between MOC-variations and the NAO becomes less, and, that the response of the MOC to the NAO is dominated by faster, wind-driven processes, compared to the response of the MOC to heat flux variations associated with changes in deep-water formation. For this reason the poleward shift is associated with an increased timescale of MOC-variability. Because the relation between the MOC and the NAO is not robust in climate models, we hesitate to attribute a model-independent generality to the increase of timescale that occurs in the present model when the relation between the MOC and subpolar heat flux variations intensifies.

It seems worthwhile to notice that in the model shown here that the overall MOC decrease resulting from global warming is partly counteracted by an increase found at higher latitudes. This result may provide some clue as to why the MOC is more insensitive to global warming in most current climate models compared to previous more idealized studies with convection sites consisting of very few grid cells that are more likely to completely switch off as a consequence of perturbations than in a case where the convection sites can shift to new locations.

For the present experiment we come to the following conclusions:

Future changes in convection drive a poleward shift of MOC-variability.

Old convective areas become more irregular, resulting in larger decadal variability in heat flux and deepwater formation.

Associated with sea-ice melt, new convective areas arise, that lie further poleward, and in which heat flux variability becomes associated with MOC-variability instead of sea-ice variations.