1 Introduction

The subpolar gyre (SPG) influences the North Atlantic climate by modulating the transport of heat and salt between the North Atlantic and the Arctic oceans and, particularly, into the Nordic and Labrador seas, where deep water formation takes place (e.g., Grossmann and Klotzbach 2009). Variations in the strength and shape of the SPG can lead to major changes in, for example, the intensity of the Atlantic meridional overturning circulation (AMOC) (e.g., Hátún et al. 2005; Moreno-Chamarro et al. 2015) or the distribution of sea ice in the Arctic (e.g., Yoshimori et al. 2010; Jungclaus et al. 2014). A recent paleoceanographic reconstruction of the Atlantic mid-depth (600–1300 m) gyre circulation over the past 1500 years suggested that the SPG weakened between the Medieval Climate Anomaly (MCA; roughly the 950–1250 CE period; years hereafter always refer to the Common Era) and the Little Ice Age (LIA; roughly 1450–1850) (Copard et al. 2012). Although the abruptness of the SPG weakening could not be well defined due to the relatively poor temporal resolution of the proxy used, the finding is consistent with other climate reconstructions from the North Atlantic that do indeed present an abrupt—in some decades—change in the thirteenth and fourteenth centuries: climate reconstructions, for example, of upper-ocean temperature (e.g., Sicre et al. 2008, 2014; Ran et al. 2011; Miettinen et al. 2012, 2015; Moffa-Sánchez et al. 2014a) and salinity in the North Atlantic (e.g., Jensen et al. 2004; Moffa-Sánchez et al. 2014a), summer temperatures of northeastern North America (Gennaretti et al. 2014), ice cap growth of the Arctic Canada or off Iceland and East Greenland (Massé et al. 2008; Miller et al. 2012; Miettinen et al. 2015), and of the inflow of Atlantic waters into the Nordic Seas (Dylmer et al. 2013). Yet attribution and the physical mechanism responsible for such an abrupt SPG weakening remain unclear, as well as the connection between this event and the other climate changes reconstructed in the North Atlantic during the LIA. This study thus aims at better understanding the SPG variability on decadal and longer time scales, its related dynamics, and how it contributed to shaping the North Atlantic climate during the last millennium.

A reduction in the meridional heat transported by the AMOC was invoked in earlier studies to explain anomalous cold conditions in the North Atlantic reconstructed during the LIA (e.g., Lund et al. 2006; Miller et al. 2012; Lehner et al. 2013). This hypothesis, however, found no support in a recent reconstruction of the AMOC strength (Rahmstorf et al. 2015). Other studies suggested that the cold North Atlantic could have been associated with a shift in the phase of the North Atlantic Oscillation (NAO) (Trouet et al. 2009; Copard et al. 2012), as it is well known that the NAO influences the SPG variability through wind forcing and surface buoyancy fluxes (e.g., Eden and Willebrand 2001; Deshayes and Frankignoul 2008). The NAO reconstruction of Trouet et al. 2009 was nevertheless contested in a later study Lehner et al. (2012), and Ortega et al. (2015) presented an alternative NAO reconstruction in which no phase shift was found between the MCA and LIA. Further uncertainty arises from the fact that the previous studies attributed these North Atlantic climate anomalies during the LIA to either internal climate variability (Drijfhout et al. 2013) or external forcing, including a series of strong, decadally-paced volcanic eruptions (e.g., Miller et al. 2012; Schleussner and Feulner 2013) and changes in the incoming solar radiation during the thirteenth century (Lehner et al. 2013).

Here, we use an ensemble of three transient coupled climate simulations of the last millennium to investigate the mechanisms of decadal to centennial SPG variability under pre-industrial conditions and their implications for the multicentennial climate variability in the North Atlantic. This study concentrates on the dynamics of an abrupt weakening shift in the SPG strength that is simulated in one ensemble member; such a shift separates two relatively stable, multicentennial North Atlantic climate regimes, so that post-shift features resemble those reconstructed during the LIA. We thus aim to answer the following questions: How does this abrupt SPG transition compare with other dominant features of simulated SPG strength variability on decadal and longer time scales? (Sect. 3). What are the main differences between the two climate regimes before and after the shift? What mechanism allows the later regime to be stable for at least two centuries? (Sect. 3.1). And, what are the potential triggering factors of such an abrupt shift? (Sect. 3.2.1). We also use sensitivity experiments to assess the relative role of the initial climate conditions and the imposed volcanic forcing during the shift (Sect. 3.2.2). Our results are discussed and compared against some of the available reconstructions for the LIA in Sect. 4. Conclusions are given in Sect. 5.

2 Model description and experimental setup

This study employs the Max Plank Institute Earth System Model for paleo-applications (MPI-ESM-P). The atmosphere general circulation model ECHAM6 (Stevens et al. 2013) is run at a horizontal resolution of spectral truncation T63 (1.875°), with 47 vertical levels that resolve the stratosphere up to 0.01 hPa. The ocean/sea-ice model MPIOM (Marsland et al. 2003; Jungclaus et al. 2013) has a conformal mapping horizontal grid of nominal 1.5° resolution and grid poles over Antarctica and southern Greenland. The convergence of the mesh size toward the poles translates into a grid spacing of 15 to 100 km in the North Atlantic and, thereby, provides a relatively high resolution in the regions of interest for this study, namely the northern North Atlantic. Vertically, MPIOM applies 40 unevenly spaced z levels, with the first 20 covering the upper 700 m of the water column.

We use an ensemble of three transient simulations of the preindustrial last millennium, between 850 and 1849 (hereafter Past1000-R1, R2, and R3). These follow the protocol of the Paleoclimate Modeling Intercomparison Project, Phase 3 (PMIP3) (Schmidt et al. 2011; Braconnot et al. 2012). Prescribed external forcing factors are: variations in volcanic aerosols (Crowley and Unterman 2012), in total solar irradiance (TSI) (Vieira et al. 2011) and in the atmospheric concentration of the most important well-mixed greenhouse gases and anthropogenic aerosols (Schmidt et al. 2011); changes in global land-cover and agricultural areas (Pongratz et al. 2008); and annual values of the orbital parameters, that is, eccentricity, obliquity, and perihelion. In addition, we use an 1156-year-long control simulation under constant preindustrial (1850) boundary conditions (hereafter PiControl) as a reference for the forced experiments. The integration of the Past1000s is initiated after a 400-year-long spin-up with constant 850 boundary conditions, itself initiated from the last integration year of PiControl. Further information on the model and the experimental setup is provided by Jungclaus et al. 2014.

Three ensembles of sensitivity experiments are, in addition, created to assess the robustness of the simulated SPG shift under different initial conditions and applied external forcing. The first ensemble focuses on long-term preconditioning and consists of 5 all-forcing simulations initiated in 1593 (hereafter AllFor-1 to 5). This ensemble includes the original Past1000-R3 as AllFor-1. The second ensemble focuses on short-term preconditioning and consists of 5 all-forcing simulations initiated in 1600 (hereafter AllFor-6 to 10). The third ensemble focuses on the role of the volcanic forcing and consists of 10 simulations initiated in 1593 in which there is no volcanic forcing for the period 1593–1650 (hereafter NoVol-1 to 10). All simulations are initiated in January 1st, and are run until, at least, 1699. Each ensemble member is generated from the original Past1000-R3 (AllFor-1) by inducing an infinitesimal perturbation in the atmospheric vertical diffusivity during the first integration year, except for the first member of the NoVol ensemble. This method for ensemble initialization is a common procedure with MPI-ESM (e.g., Zanchettin et al. 2013), as it produces a rapid divergence of climate trajectories and, therefore, an adequate ensemble spread.

3 Results

We use an index for the SPG strength to describe its variability above decadal time scales during the last millennium. This index is herein defined as the absolute value of the barotropic streamfunction spatially averaged between 50°N–65°N and 10°W–60°W, further smoothed with an 11-year running mean filter (Fig. 1a). Overall, the SPG variations in all transient simulations fall within the range of variability in PiControl (Fig. 1a; dot-dashed and dotted lines). Two remarkably different behaviors can nevertheless be identified in the Past1000s: whereas the SPG weakens progressively at an almost constant rate throughout the entire last millennium in R1 and R2 (with trends of the order of −0.10 Sv/century), the SPG in R3 shows a rapid shift from a strong to a weak state around 1600. An algorithm for discrete level identification (Little et al. 2011) allows an objective characterization of the shift: it occurs in about two decades, accounts for a weakening of 0.72 Sv, and separates two stationary states, each covering at least 250 years (Fig. 1a, horizontal solid lines for Past1000-R3, and Fig. 6). Discrete levels in Past1000-R1 and R2, in contrast to R3, succeed each other more rapidly in time and are separated by smaller steps, a behavior that corresponds well with a smooth, long-lasting weakening. Levels in the three Past1000s take similar values during the last two centuries, ranging between the 5th and the 25th percentiles of the SPG strength in PiControl. A weak SPG state thus robustly develops throughout the last millennium in the three simulations, either smoothly or abruptly.

Fig. 1
figure 1

a Index for the SPG strength (in Sv) during the last millennium: light-gray lines are annual averages, and solid black lines are the corresponding 11-year running mean. Horizontal dash-dotted and dashed lines respectively account for the 25th–75th and 5th–95th percentiles of the SPG strength in PiControl. Horizontal solid lines are the discrete levels extracted following the algorithm described in Little et al. (2011). b Annual globally averaged values of the reconstructed TSI (black line, in W/m2) (Vieira et al. 2011) and the simulated volcanic forcing, shown as net adiative flux anomalies, long plus short waves (LW + SW, in W/m2), at the top of the atmosphere with respect to the climatological mean. The volcanic forcing is for Past1000-R3. Timing of volcanic eruptions is the same in all Past1000s as they are all based on the reconstructed volcanic aerosol properties by Crowley and Unterman (2012), although the implementation of the forcing slightly differs among the simulations (see Jungclaus et al. 2014)

3.1 North Atlantic climate regimes before and after the shift

This section compares the two climate regimes in the North Atlantic before and after the shift in the SPG strength in Past1000-R3, and explores the mechanism that renders the later regime stable for at least 250 years. The analysis is performed assuming the pre-shift period as reference regime and the post-shift period as the anomalous one. Figure 2 illustrates the differences between the mean states after (1650–1749) and before (1450–1549) the shift. These two periods are chosen to exclude years immediately around the SPG shift, when climate anomalies are exceptionally large (see, for instance, Fig. 3). Similar results are obtained if longer period are considered (e.g., for the period 1350–1549 and 1650–1849).

Fig. 2
figure 2

Differences between the 100-year climate means after (1650–1749) and before (1450–1549) the shift in Past1000-R3. Only significant anomalies above the 95 % confidence level are shown, based on the likelihood of occurrence of the signal in PiControl. a Barotropic streamfunction (in Sv, shading) with contours representing the climatological mean in PiControl at 10 Sv intervals (dashed/solid lines correspond to cyclonic/anticyclonic flow). b SSS contribution to the ocean surface density (in kg/m3), assuming a linear relationship between both variables. c Late-winter (March) mixed layer depth (in meter; shading), with contours representing the climatological mean in PiControl at 750 m intervals. d SST contribution to the ocean surface density (in kg/m3), assuming a linear relationship between both variables. Note that color scale is the same as in (b). e AMOC (in Sv), with contours representing the climatological mean in PiControl at 5 Sv intervals. f Sea ice concentration in March (in percentage of area). g Wintertime (DJF) SLP (in hPa; shading), with contours representing the climatological mean in PiControl at 5 hPa intervals. h Wind stress curl (in 10−7 N/m3; shading), with lines accounting for the zero line before (red) and after (blue) the shift

Fig. 3
figure 3

ag Standardized anomalies of different variables in Past1000-R3. Standardization is done after applying an 11-year running mean filter. Values in parentheses correspond to the mean ± the standard deviation in the last millennium. a SPG strength (12.77 ± 0.80 Sv), defined as in Fig. 1. West–east gradient of upper-ocean densities (see text, Sect. 3.1; 0.285 ± 0.055 kg/m3). b Late-winter (March) mixed layer depth (MLD) and SSS in the Labrador Sea, both averaged between 45°N–63°N and 40°W–60°W (749 ± 165 m and 33.46 ± 0.14 psu, respectively). c Meridional overturning index (MOI), defined as the meridional overturning streamfunction averaged between 35°N–45°N at 1000 m depth (20.89 ± 0.79 Sv). Meridional heat transport (HTTR) in the Atlantic Ocean at 26°N (0.975 ± 0.030 PW). d Atlantic Ocean HTTR at 55°N (0.485 ± 0.22 PW). Nordic Seas SST, averaged between 65°N–80°N and 20°W–20°E (3.89 ± 0.40 °C). e Late-winter (March) sea ice concentration in the North Atlantic sector, averaged between 50°N–85°N and 80°W–40°E (25.4 ± 1.3 % of area). Total (liquid + solid) freshwater transport (FWTR) through the Denmark Strait (4309 ± 370 km3/year). f Wind stress curl averaged over the subpolar region (see text, Sect. 3.1; (1.559 ± 0.075)·10−7 N/m3/s). g NAO index, defined as the first mode of the wintertime (DJF) SLP over the North Atlantic, between 20°N–90°N and 80°W–40°E and after applying an 11-year running mean filter (61 % of the explained total variance; anomaly pattern not shown). h Northern Hemisphere (NH) near-surface temperature anomalies (in  °C) with respect to the period 1500–1849: overlap of the reconstructions [in  %; shading, from Masson-Delmotte et al. (2013)] and simulated in Past1000-R1, R2, and R3 (lines, smoothed with 11-year running mean)

The reduced SPG strength after the shift is the result of a broad slowdown of the subpolar barotropic circulation, with the exception of the Irminger gyre (Fig. 2a). The center of the subtropical gyre is also significantly weaker during the later regime, although most of the changes in the gyre circulation lie within the subpolar region. Despite post-shift changes in the barotropic streamfunction, a meridional shift in the location of the boundary between the subtropical and the subpolar gyres does not occur after the shift, as evidenced by the fact that the location of the zero wind-stress-curl line stays unchanged between both regimes (Fig. 2h). Moreover, the abrupt shift in the SPG strength can be detected in the second variability mode of the North Atlantic barotropic streamfunction above decadal time scales; this mode resembles the mean circulation in the region (not shown), and is linked with oceanic convection activity in the Labrador Sea, the SPG strength itself, and the northward heat transport at subpolar latitudes (not shown).

Broadly speaking, the strength of the SPG is largely controlled by the horizontal gradient in upper-ocean densities between the gyre’s center and its boundaries, and by the wind stress over the subpolar area (e.g., Häkkinen et al. 2011; Born and Stocker 2014). The former is herein represented as the difference between upper-ocean (500 m) densities in the western subpolar basin, averaged between 55°N–60°N and 50°W–60°W, and the eastern basin, averaged between 55°N–60°N and 10°W–20°W, further smoothed with an 11-year running mean (Fig. 3a). The zonal gradient in upper-ocean densities thus seeks to represent the difference between the gyre’s center (or maximum) within the Labrador Sea, where upper-ocean density reaches maximum values, and the gyre’s eastern boundary (contours in Fig. 2a). The zonal density gradient shows relatively high correlation with the SPG strength (maximum correlation coefficient of 0.76, p < 0.05, at lag 0 years), which indicates that it largely controls SPG variations above decadal time scales. On the other hand, the contribution of the wind stress curl is herein accounted by spatially averaging the wind stress curl over the subpolar region, between 50°N–65°N and 10°W–60°W, also smoothed with an 11-year running mean (Fig. 3f). The wind stress curl shows a significant but relatively weak connection with the SPG strength on decadal and longer time scales for the whole millennium (maximum correlation coefficient of 0.34, p < 0.05, with the wind stress curl leading by 3 years), although it can occasionally be linked to notable variations in the SPG, for example for the period 1350–1400 (see Fig. 3a, f). The reorganization of the climate after 1600 leads to changes in both the upper-ocean density gradient and the wind stress curl, as discussed in the following.

The zonal gradient of upper-ocean densities appears strongly weakened as a result of the anomalous fresh, hence light, surface conditions in the western subpolar North Atlantic after 1600 (Fig. 2b), whereas post-shift changes in upper-ocean densities in the eastern subpolar North Atlantic are much smaller. In the Labrador Sea, in particular, the sea-surface salinity (SSS) drops by about 0.20 psu between regimes (Fig. 3b) and dominates the density increase due to surface cooling (compare Fig. 2b and d). This surface freshening results, first, from an increase in the export of sea ice from the Arctic Ocean and Nordic Seas into the Labrador Sea (Fig. 3e, discussed below), where it melts and freshens the surface layers, and second, from a reduced westward salt transport along the Irminger Current (not shown) due to weaker gyre circulation. The SPG therefore self-maintains its weaker regime after 1600 by keeping the western subpolar basin anomalously fresh, hence light.

Changes in the deep water formation in the western subpolar North Atlantic also contribute to sustaining the weak SPG after the shift. In its climatological mean, late-winter cooling of the upper Labrador Sea induces strong vertical mixing down to a depth of several hundred meters (Fig. 2c, contours) that forms a core of very dense waters within the SPG’s center, which in turn strengthens the zonal gradient of upper-ocean densities and, thereby, the gyre’s circulation (Born and Stocker 2014). After the SPG shift, however, ocean deep mixing is locally hampered as fresher, hence lighter, surface conditions develop in the Labrador Sea (Figs. 2c, 3b). This contributes to the reduction of the upper-ocean densities in the western subpolar region and, hence, of the strength of the SPG.

Late-winter mixed layer depth in the Labrador Sea, although weaker than pre-shift conditions, still reaches values above 2000 m after the shift. Deep waters are thus still formed within the Labrador Sea. Likewise, no remarkable changes occur in oceanic deep convection in the Nordic Seas after 1600 (Fig. 2c) or in the strength of the Greenland Scotland Ridge overflow (not shown). Hence, the SPG shift is associated with non-significant changes in the circulation pattern of the AMOC (Fig. 2e). The AMOC strength (Fig. 3c) does not exhibit any long-lasting shift after 1600, in contrast to the SPG, even though it temporally weakens between 1600 and 1640 during the shift.

Changes in the large-scale circulation of the Atlantic Ocean result in changes in the poleward oceanic heat transport. The meridional heat transport at 26°N is only slightly reduced after the SPG shift, because it is largely driven by the AMOC (Fig. 3c). At subpolar latitudes, where the meridional heat transport is principally driven by the SPG, the meridional heat transport shows a conspicuous reduction from 1600 onwards (Fig. 3d); a smaller amount of upper-ocean heat reaching the intergyre region is consequently transported poleward along the eastern SPG. On the contrary, heat tends to accumulate at mid-latitudes in the upper North Atlantic after the shift, thus resulting in an upper-ocean (upper 1000 m approximately) warming along the path of the North Atlantic Current up to 50°N–55°N (Fig. 2d). Farther downstream, a deficit in upper-ocean heat leads to significantly colder surface conditions in both the western subpolar region and Nordic Seas, where, in particular, sea-surface temperatures (SSTs) drop by about 0.4 °C on average after the shift (Fig. 3d). Near-surface atmospheric temperature anomalies between the two climate regimes mirror this pattern (similar to the one shown in Fig. 8d). The SPG shift alone is thus able to induce large-scale, long-lasting colder surface conditions in the North Atlantic sector.

Cooling of the upper Nordic and Labrador seas leads to a post-shift overall expansion and thickening of the sea ice in the North Atlantic and Arctic oceans (Fig. 3e); for the sea ice extent this is particularly prominent in the Labrador, Nordic, and Barents seas (Fig. 2f). As a result, a larger amount of sea ice is exported toward the Labrador Sea after the shift. The total (liquid plus solid) freshwater transport (FWTR) through the Denmark Strait hence increases by about 200 km3/year on average after the shift (Fig. 3e), mainly due to an increase in the solid component (not shown). Larger FWTR into the Labrador Sea contributes to the surface freshening (Fig. 2b) and, thereby, helps sustain the weaker SPG regime after the shift. Besides, the insulating effect due to expanded sea ice in the northern Labrador Sea weakens the underlying surface circulation along the West Greenland and Labrador currents (Fig. 2a).

Colder surface conditions and expanded and thicker sea ice reduce surface heat losses in the Labrador Sea (not shown). This increases stability in the overlying atmosphere after the shift (Hoskins and Karoly 1981), inducing significant positive anomalies in winter sea-level pressure (SLP) over the western subpolar North Atlantic that describe a blocking-like structure (Fig. 2g). The Icelandic Low appears displaced eastward from its climatological position after the shift, as illustrated by the significant negative anomalies in the winter SLP to the east of its center. This anomalous SLP pattern does not resemble any particular NAO phase, and the NAO index does not exhibit any regime change associated with the SPG shift (Fig. 3g). The development of an anomalous blocking-like SLP structure south of Greenland after the shift causes a significant reduction of the wind stress curl over the subpolar basin, particularly to the west (Figs. 2h, 3f), therefore contributing to the SPG weakening as well. This illustrates the role of the atmosphere in the weaker SPG regime from 1600 onwards.

In summary, two contributions mainly sustain the weak SPG regime after the shift: first, a long-lasting upper-ocean freshening in the Labrador Sea, which reduces the zonal gradient in upper-ocean densities in the subpolar North Atlantic and is maintained by the own weak SPG and by larger export of sea ice from the Nordic Seas and the Arctic; and second, a long-lasting reduction in wind stress curl due to an anomalous blocking-like structure over the subpolar North Atlantic.

3.2 Mechanism triggering the shift in the subpolar gyre strength

This section first describes the chain of events that leads to the development of the SPG shift in Past1000-R3, and then highlights the possible triggering mechanisms. Our dynamical interpretation is confirmed by means of sensitivity experiments.

3.2.1 The SPG shift in Past1000-R3

The volcanic cluster that ends in 1600 with the eruption of the Huaynaputina volcano in Peru (Verosub and Lippman 2008) marks the onset of the shift between North Atlantic climate regimes described in Sect. 3.1. This cluster is composed of volcanic eruptions of relatively small magnitude, compared to other major and better-known events for the last millennium (Fig. 1b). Nonetheless, the cluster accumulates a negative anomaly in the globally-averaged net radiative flux at the top of the atmosphere (TOA) down to −5 W/m2 right after 1600, and of approximately −2 W/m2 during the following 20 years (Fig. 4a). Over the 10 years after these eruptions, the Northern Hemisphere total sea ice volume increases by 10 % from 28,000 km3 to 30,800 km3 (Fig. 4a). The increase is remarkable north of Greenland, where the sea ice thickens, and in the Barents and Nordic seas, where the sea ice edge advances extensively southward (not shown). An anomalous, albeit of smaller amplitude, growth in the Arctic sea ice is found in this same simulation also after the strong volcanic eruptions in 1250 or 1815 (Fig. 3e, for the North Atlantic sector). Arctic sea ice growth typically occurs in all three Past1000s after some, though not all, volcanic eruptions (Supplementary Material in Moreno-Chamarro et al. 2015), in agreement with earlier studies (e.g., Zhong et al. 2011; Zanchettin et al. 2012). Other relatively abrupt events of sea ice growth are, in contrast, attributed to internal cold events in the subpolar North Atlantic driven by a SPG weakening alone, as for example around 1150 or 1350, and should not be mistaken for a response to strong volcanic eruptions (Moreno-Chamarro et al. 2015).

Fig. 4
figure 4

Standardized annual anomalies of variables in Past1000-R3 throughout the SPG shift, between 1550 and 1680 (thin lines), with thick lines as their corresponding 11-year running means. Values in parentheses correspond to the mean ± the standard deviation in the last millennium. a Northern Hemisphere (NH) late-winter (March) sea ice volume (28,000 ± 1700 km3). Also, simulated volcanic forcing (thin line), as in Fig. 1, and accumulated volcanic forcing (thick line), both in W/m2. b Liquid and solid components of the FWTR through the Denmark Strait (3192 ± 436 m3/year and 1115 ± 402 km3/year, respectively). c Mixed layer depth (MLD) and SSS in the Labrador Sea, for the same region as in Fig. 3 (748 ± 273 m and 33.46 ± 0.19 psu, respectively). d SPG strength and west–east gradient of upper-ocean densities, for the same regions as in Fig. 3 (12.7 ± 1.1 Sv, and 0.286 ± 0.076 kg/m3, respectively)

The abrupt (in just a decade) rise in the sea ice translates into an abrupt increase in the sea ice export from the Arctic Ocean and the Nordic Seas toward the subpolar North Atlantic through the Denmark Strait (Fig. 4b). The solid component of the FWTR, that is sea ice, reaches a maximum positive anomaly of about 900 km3/year in the 1620s, and remains predominantly positive thereafter until the 1660s. As also shown in Born et al. 2010, freshwater bound in ice can very effectively impact both oceanic deep convection and the SPG strength because ice remains concentrated at the ocean surface during the advection, and because its maximum export is reached in late winter, during the main convection season in the subpolar North Atlantic. As the sea ice flows southward, warmer surface waters melt it progressively, decreasing the surface salinity of the upper subpolar North Atlantic and, especially, of the Labrador Sea (Fig. 4c). The largest salinity reduction in Past1000-R3, of about 0.5 psu, thus occurs around 1620, synchronously with the largest freshwater export from the Arctic. Fresher, hence lighter, conditions in the upper Labrador Sea considerably reduce the west-east gradient of upper-ocean densities in the subpolar region (Fig. 4d). The SPG thereby weakens by about 5 Sv until a minimum strength of nearly 9 Sv (Fig. 4d). Increase in water column stratification due to surface freshening additionally induces a strong reduction in late-winter oceanic convective mixing: minimum values of only 100 meters for the mean mixed layer depth in the Labrador Sea are found between 1620 and 1630 (Fig. 4c). Deep water formation is thus shut down, leading to substantial weakening of the AMOC for about 20–30 years (Fig. 3c).

The development of the initial anomalies in the Labrador Sea triggers two positive feedbacks within the SPG that lead to further weakening of the gyre. Firstly, a debilitated SPG transports less heat poleward between 1600 and 1620 (Fig. 3d). As a result, the upper Nordic Seas cool down, inducing anomalously strong sea ice growth in the area. Consequently, the sea ice export through the Denmark Strait reinforces (Fig. 4b), leading to further salinity drop in the upper Labrador Sea (Fig. 4c). Secondly, a weaker SPG also transports less salt toward the western subpolar basin (not shown), thus sustaining the anomalously fresh conditions in the upper Labrador Sea. As a weaker SPG also transports less heat westward, temperatures in the upper Labrador Sea are eventually cold enough to increase upper-ocean densities, despite the haline-driven lightening (from 1620 onwards, Fig. 4d). As upper-ocean density increases, deep water formation progressively resumes, mixing waters from the upper and deeper levels. This process reduces the magnitude of the negative salinity anomaly in the upper Labrador Sea, hence strengthening the gradient of upper ocean densities and the SPG. Yet, the recovery of the SPG strength is not complete, as sufficient salinity anomalies remain in the western subpolar basin to keep a weaker SPG regime stable afterward, as discussed in Sect. 3.1.

3.2.2 Testing the mechanism triggering the shift

The sensitivity experiments described in Sect. 2 are designed to assess the relevance of the long-term (AllFor-1 to 5) and short-term preconditioning (AllFor-6 to 10) for the development of the SPG shift, as well as to determine whether the volcanic cluster peaking in 1600 is necessary to trigger this shift (NoVol-1 to 10). Note that AllFor-1 corresponds to the simulation that originally presents the shift (Past1000-R3), hence the target to reproduce. The SPG strength for all simulations is plotted in Fig. 5, overlying their associated discrete levels. The change in the SPG strength between the two discrete levels in years 1620 and 1590 is additionally shown in Fig. 6 (along the horizontal axis). The choice of these particular years is arbitrary as levels extend for more than one. Although we are mainly interested in shifts occurring around 1600, as in the target, we also include the delayed shift in AllFor-8 around 1640 (Fig. 6, asterisk; using discrete levels in years 1640 and 1670) for it occurs after a large volcanic eruption (Fig. 1b) and exhibits features similar to the target (see below).

Fig. 5
figure 5

Index for the SPG strength (in Sv) in the sensitivity experiments, as in Fig. 1: thin and solid lines are respectively annual averages and the 11-year running mean. Horizontal solid lines are the discrete levels. Vertical dotted lines indicate the starting year of each integration

Fig. 6
figure 6

Scatter plots for the sensitivity experiments: a and b Change between the two discrete levels in years 1620 and 1590 in the Nordic Seas SST (in  °C) and Labrador Sea SSS (in psu) against the analogous change in the SPG strength (in Sv). Variables are defined for the same regions as in Fig. 3. c 30-year trend, between 1600 and 1629, in the liquid (crosses) and solid (diamonds) FWTR through the Denmark Strait (in km3/year per decade) against the change between discrete levels in years 1620 and 1590 in the SPG strength (in Sv). The FWTR trend of the later shift in AllFor-8 is calculated between 1640 and 1669 (asterisk in a and b), and between 1590 and 1619 in AllFor-7, since this simulation undergoes a faster increase than the others do (first cross starting from the left in a and b)

The magnitude of the change in the SPG strength after 1600 varies widely across the AllFor simulations, ranging between 0.05 and 0.78 Sv (Figs. 5, 6). Values are generally larger in the short-term preconditioning ensemble (AllFor-6 to 10), ranging from 0.21 to 0.78 Sv, than in the long-term preconditioning ensemble (AllFor-2 to 5), between 0.05 and 0.31 Sv, except for the 0.72 Sv shift in the target (AllFor-1). Only AllFor-7, AllFor-8, and AllFor-9, which all belong to the short-term preconditioning ensemble, show a SPG shift larger than 50 % of the target. The chance of reproducing a shift as strong as the target increases in experiments initiated in 1600, presumably because they are initiated closer to the onset of the shift, and their climate trajectories therefore cannot sufficiently diverge. On the other hand, four AllFor simulations show a relatively small reduction in the SPG strength, below 0.25 Sv. The substantial differences found in the SPG evolution across the experiments suggest that the volcanic forcing alone does not suffice to robustly trigger a SPG shift. The volcanic cluster must then occur together with a favorable combination of background climate conditions.

Within the NoVol ensemble, a SPG shift occurs only in NoVol-3 (0.62 Sv) and in NoVol-6 (0.36 Sv). The volcanic forcing is hence not a necessary factor to trigger a shift in the SPG. However, since most of the NoVol experiments show no shift after 1600 at all, we can also conclude that, within the limits of our small ensemble size, the volcanic cluster seems to increase the likelihood of the occurrence of a SPG shift.

Discrete levels are also obtained for the Labrador Sea SSS and Nordic Seas SST in the sensitivity experiments, for the same regions as in Fig. 3 (time series not shown). Their associated changes around the shift plotted against the change in the SPG strength (Fig. 6a, b, respectively for SSS and SST) allow evaluating the magnitude of changes in the Nordic Seas with respect to those within the subpolar North Atlantic, and assessing whether the shift simulated in each experiments involves the entire North Atlantic/Arctic region, as shown to be the case of the target (see Sect. 3.1). According to Fig. 6, the larger the change in the SPG strength is, the colder and fresher conditions develop, respectively due to reduced heat and salt transports by the SPG, in the upper Nordic and Labrador seas after the shift. Changes in the SPG strength and Nordic Seas SST appear to follow a non-linear relationship, with an even colder Nordic Seas surface for the strongest SPG shifts. This result might be related to the presence of positive feedbacks between the subpolar North Atlantic and the Nordic Seas involving the upper-ocean temperatures in the Nordic Seas, the sea ice transport through the Denmark Strait, and the SPG strength. Six simulations experience changes in the Labrador and Nordic seas larger than 50 % of the target; most simulations exhibit much weaker changes or none at all. Furthermore, only AllFor-7 shows changes in the SPG strength and Nordic Seas SST larger than the target, whereas none of the experiments shows changes in the Labrador Sea SSS larger than the target. The probability of simulating a shift in the North Atlantic climate at least as strong as the target is thus relatively small (10 %, or 2 over 20). The North Atlantic climate will most likely experience relatively small transitions even if it is following a trajectory close to those featuring an abrupt and strong shift between two regimes. This could explain why the SPG and, in general, the North Atlantic undergo a smooth transition throughout the past millennium instead of a strong shift in Past1000-R1 and R2 (Fig. 1a).

As described before, the target SPG shift in AllFor-1 (Past1000-R3) is driven by anomalous large export of freshwater from the Arctic into the Labrador Sea. This is illustrated by the 30-year trend between 1600 and 1629 in the total (liquid plus solid) FWTR through the Denmark Strait, of about 360 km3/year per decade (Fig. 6c). Similar trends are also found in the other five simulations featuring the largest SPG shifts (Fig. 6c), although these can be due to increases mainly in the solid component—as in the target—, in the liquid one, or in both components simultaneously (Fig. 6c). The relative importance of each component remains unclear and can relate, for example, to the region where deep convection preferentially occurs in the model and during the years when the discharge of freshwater occurs. Regardless of the case, and as shown for the target (Fig. 4), such increase in the Arctic freshwater export leads to a broad freshening of the upper Labrador Sea and, in turn, to a substantial, long-lasting SPG weakening, which drives further reduction in the salt transport along its northern rim. All simulated shifts are therefore triggered by an increase in the FWTR, independently of the type of experiment considered.

In those experiments with a small or absent shift, the 30-year trend in the total FWTR varies between approximately ±500 km3/year per decade (Fig. 6c). In some cases, small SPG shifts are even associated with a decrease in one or both components of the FWTR. For similar values of the FWTR trend, for example between 300 and 450 km3/year per decade, the change in SPG strength ranges between 0 and 0.78 Sv. This suggests, first, that the SPG is not too sensitive to the FWTR from the Arctic in our model, and, second, that at least the positive climate feedbacks described before must get activated after the discharge of freshwater into the Labrador Sea in order to trigger a SPG shift.

4 Discussion

4.1 Attribution of the SPG changes in the last millennium

The millennial-scale weakening of the SPG strength during the past millennium develops smoothly in Past1000-R1 and R2, and abruptly in R3. Such weakening must be driven by the imposed external forcing throughout the entire last millennium because the SPG strength in PiControl does not show any long-term trend (not shown). However, differences across the Past1000s suggest that any influence of the external forcing is shaped by internal climate variability. In none of the transient simulations, the SPG strength shows significant cross-correlation with the solar or the volcanic forcing (Fig. 1b for these two forcings; cross-correlation profiles not shown). A clear response to strong volcanic eruptions was found neither in the North Atlantic climate, nor in the SPG strength (Supplementary Material in Moreno-Chamarro et al. 2015). This contrasts with previous studies associating weaker SPG with periods of low solar irradiance (Moffa-Sánchez et al. 2014b) or anomalous (weaker or stronger) SPG to volcanic forcing (Zhong et al. 2011; Zanchettin et al. 2012; Schleussner and Feulner 2013). Zanchettin et al. 2013, nevertheless, showed that multiple response pathways to short-lived radiative perturbations like volcanic eruptions can exist, whose activation crucially depends on background climate conditions. By extension, a linear and unique response of the SPG to volcanic eruptions should not necessarily be expected. The simulated long-lasting trend in the SPG strength might also be attributed to varying orbital parameters (Kaufman et al. 2009). A detailed investigation on the impact of the external forcing on the SPG dynamics at millennial time scales is however beyond the scope of this study.

The sensitivity experiments highlight the role of internal climate variability at the onset of the SPG shift in Past1000-R3, albeit the volcanic forcing increases the likelihood of occurrence of a SPG shift. These results contrast with previous studies in which similar abrupt events were attributed to the volcanic forcing alone in climate simulations of the past millennium (Miller et al. 2012; Schleussner and Feulner 2013; Schleussner et al. 2015). The sensitivity experiments used here, however, still include changes in the other imposed external forcings (see Sect. 2). Among these, only solar activity and volcanic eruptions have been found relevant to decadal climate variability (Fernández-Donado et al. 2013). Were the SPG shift in Past1000-R3 driven by changes in solar activity, rather than by the volcanic forcing, it would have occurred during the maximum around 1600 (Fig. 1b), in opposition to earlier studies associating weaker SPG with solar minima in the past millennium (Moffa-Sánchez et al. 2014b). It is worth remarking, however, that the late SPG shift found in one of the experiments does not occur during this TSI maximum, but after the Mount Parker’s eruption in 1641 (All-For-8 in Fig. 5). TSI changes thus seem to be not relevant for the inception of the SPG shift.

4.2 The Medieval Climate Anomaly: Little Ice Age transition

Climate changes in the North Atlantic simulated after the shift in the SPG strength and reconstructed during the LIA are indeed in broad agreement. In particular, Copard et al. 2012], based on the Nd isotopic composition (εNd) of deep-sea corals, concluded that the SPG underwent a transition from an enhanced phase during the MCA to a weak phase during the LIA in the course of the thirteenth and fourteenth centuries (Fig. 7a). In accordance with such reconstructed SPG weakening, several SST records based on different proxies show a rapid cooling at the beginning of the LIA, for instance, in the Nordic Seas (e.g. Berner et al. 2011) (Fig. 7c), off North Iceland (e.g., Sicre et al. 2008; Ran et al. 2011) (Fig. 7b), or in the eastern Labrador Sea (Moffa-Sánchez et al. 2014a; Fig. 7e). Other SST reconstructions show, in contrast, warmer conditions during the LIA than during the MCA, for example, in the eastern subpolar North Atlantic (e.g., Miettinen et al. 2012; Fig. 7d) or off northeastern Newfoundland (Sicre et al. 2014). These SST reconstructions thus draw a characteristic pattern of cold Nordic and Labrador seas opposed to a warm eastern subpolar North Atlantic and North Atlantic Current that resembles the simulated post-shift SST anomaly pattern (Fig. 7h) and that, in the model, results from a broad redistribution in upper-ocean heat content by the SPG alone (see Sect. 3.1). The simulated reduction in the oceanic meridional heat transport at subpolar latitudes after the SPG shift (Fig. 3d) also corresponds well with the reconstructed reduction in the inflow of warm and saline Atlantic waters into the Nordic Seas during the MCA-LIA transition (Dylmer et al. 2013). Likewise, simulated post-shift cold and fresh upper Labrador Sea conditions (Figs. 2, 3) agree with those oceanographic reconstructions in which an increasing influence of the Arctic waters on subpolar latitudes was found (Jensen et al. 2004; Ran et al. 2011; Moffa-Sánchez et al. 2014a) (Fig. 7f). Last but not least, simulated expanded sea ice conditions after the SPG shift (Figs. 2f, 3e) match reconstructed changes in the sea ice-proxy IP25 from the North Icelandic Shelf after 1300 (Massé et al. 2008; Fig. 7g) as well as the abrupt increase in ice-cap growth over the Arctic Canada and Iceland (Miller et al. 2012). Yet agreement between simulated post-shift anomalies and those reconstructed in the LIA mainly concerns their magnitude and duration but not their timing of occurrence.

Fig. 7
figure 7

Reconstructions of the North Atlantic climate variability in the preindustrial past millennium: a “Nd values from deep-sea corals in the Rockall Trough (Copard et al. 2012). b Alkenone SST reconstruction off North Iceland (in  °C) (Sicre et al. 2008). c August SST reconstruction from the Vøring Plateau (in  °C) (Berner et al. 2011). d August SST reconstruction from the Reykjanes Ridge (in  °C) (Miettinen et al. 2012). e and f Respectively, T. quinqueloba δ18O and percentage of the polar species N. paquiderma (sinistral) (Nps) from the Eirik Drift (Moffa-Sánchez et al. 2014a). g IP25 sea ice proxy from off North Iceland (Massé et al. 2008). h Simulated after (1650–1749) minus before (1450–1549) SST anomalies, as in Fig. 2 (shading, with hatching indicating nonsignificant anomalies; in  °C). Locations of the core sites of the reconstructions shown in (ag) are also indicated

The inception of the simulated and reconstructed anomalies occur in the early seventeenth century and around the thirteenth and fourteenth centuries respectively. Dominance of internal climate variability can nevertheless explain such differences in the timing and, therefore, reconcile reconstructed and simulated variability. Schleussner and Feulner (2013) reported that the onset of the LIA can occur at different times in an ensemble of simulations with the Climber-3α model, which the authors attributed to the sensitivity of the subpolar variability to minor changes in the North Atlantic freshwater budget in their model. They were thus able to simulate an earlier LIA onset (around 1450) in the North Atlantic by forcing the Nordic Seas with a constant freshwater offset. It is unknown whether MPI-ESM-P features a similar sensitivity to the freshwater export from the Arctic, but we argue that an analogous temporal spread in the onset of the shift as in Schleussner and Feulner (2013) can be expected: only when a sufficient amount of freshwater is discharged toward the subpolar North Atlantic under the appropriate climate conditions, a SPG shift can occur. In fact, some temporal spread in the onset of the shift is found in the sensitivity experiments analyzed in Sect. 3.2.2 after 1600. However, this hypothesis was not tested for a shift before 1593.

In view of these results, could the regime shift described here be regarded as onset of a LIA-like episode in the North Atlantic? To answer this question, we can compare the transient climate evolution of the three Past1000 simulations. Their global and hemispheric mean near-surface air temperatures (SAT) exhibit similar millennial negative trends and cold excursions in response to the strongest volcanic eruptions, and are similar to the reconstructed ones for this period (Fig. 3h, for the Northern Hemisphere mean SAT). Compared to Past1000-R1 and R2, the SPG shift in Past1000-R3 therefore does not shape the global or the hemispheric evolution of the climate of the past millennium; rather, climatic changes associated to the shift arise more prominently on regional scales in the North Atlantic/Arctic sector, as discussed in Sect. 3.1. The pattern of SAT anomalies between the MCA and the LIA in the North Atlantic/Arctic for the three Past1000s are relatively alike (Fig. 8). Similarities between Past1000-R3 and R1 or R2 are illustrated by the relatively high spatial correlation coefficient between the corresponding anomaly patterns (numbers in Fig. 8a, b). Interestingly, the pattern of MCA-LIA anomalies in Past1000-R3 is mainly explained by climate changes due to the shift in the SPG strength (Fig. 8d, and associated spatial correlation coefficient). The three Past1000s therefore present similar cold conditions in the North Atlantic at the end of the pre-industrial period, which we referred to as the LIA. However, while the transition from a relatively warm MCA to a relatively cold LIA develops smoothly in Past1000-R1 and R2, it occurs abruptly in Past1000-R3, in both cases mainly driven by changes in the SPG strength alone. Hence, in our model world there exits two possible realizations of the climate evolution during the past millennium, so that the actual onset of the LIA described in the reconstructions would then just be another possible realization.

Fig. 8
figure 8

Anomalies in SAT (in  °C) of the LIA (1450–1849) with respect to the MCA (950–1250) in a Past1000-R1, b Past1000-R2, and c Past1000-R3. d SAT anomalies between the last two centuries of the LIA after the SPG shift (i.e., 1650–1849) with respect to the MCA in Past1000-R3. Only anomalies that exceed the 95 % confidence level with a two-sided t test in PiControl are shown. Spatial correlation coefficients between patterns in a, b, and d, with the one in c are also respectively shown

Our dynamical interpretation of the transition toward LIA-type episodes differs from previously suggested mechanisms for other similar abrupt climate transitions in the North Atlantic (Miller et al. 2012; Drijfhout et al. 2013; Lehner et al. 2013; Schleussner and Feulner 2013; Schleussner et al. 2015). These modeling studies primarily related the cold climate conditions in the North Atlantic to a slowdown of the meridional heat transport by the AMOC, in agreement with records describing a weak AMOC during the LIA (Lund et al. 2006; Wanamaker et al. 2012). This hypothesis, nonetheless, contrasts with a recently-proposed index for the AMOC strength over the past millennium in which no weakening was reconstructed before the twentieth century (Rahmstorf et al. 2015), more in accordance with our study. A particular case that here deserves special mention is the climate regime shift in the North Atlantic at the onset of the LIA described in climate simulations with the model CLIMBER-3α (Schleussner and Feulner 2013). The proposed mechanism included a volcanically-forced sea ice expansion in the Nordic Seas, leading to a decrease in the local deep water formation and, thereby, in the overflow through the Greenland Scotland Ridge. The latter resulted, on the one hand, in a baroclinic spin-up of the SPG due to increasing recirculation of subtropical water into the subpolar region and, on the other hand, in a LIA-like basin-wide cooling due to weakened heat transport by the AMOC. These results agree with SST reconstructions from the Nordic Seas and the eastern subpolar North Atlantic (Schleussner et al. 2015; these SST reconstructions are also shown in Fig. 7c, d), although their simulated SPG spin-up differs from reconstructed evidence of a weaker SPG during the LIA (e.g., Copard et al. 2012; Moffa-Sánchez et al. 2014a), which is more in accordance with our study. In addition, whereas the SPG shift in CLIMBER-3α appears as a transition from one stable circulation regime into another unconditionally stable mode, the shift in our model (MPI-ESM-P) can rather be regarded as a modulation of a continuous succession of rapid transition, with similar ranges of variability before and after the shift. The reasons of the different behaviors of the SPG simulated by these two models are unclear. One possible explanation concerns the models’ different complexity generating different dominant mechanisms, as shown, for instance, in the case of the bistability of SPG regimes in the multi-model investigation by Born et al. (2013).

On the other hand, Copard et al. (2012) associated the SPG weakening in their reconstructions to the shift between different NAO phases during the MCA-LIA transition (Trouet et al. 2009). Ortega et al. (2015), however, proposed a new NAO index for the past millennium in which such a transition is not found, more in agreement with our results for the NAO and previous modeling studies (e.g., Lehner et al. 2012; Gómez-Navarro and Zorita 2013). All in all, the latest evidence hence suggests that cold conditions reconstructed in the North Atlantic during the LIA were driven neither by an AMOC slowdown, nor by a shift between NAO phases, but by a weakening of the SPG strength and associated meridional heat transport.

5 Conclusions

We investigate the decadal to multicentennial SPG variability and underlying mechanisms in an ensemble of three climate simulations (Past1000s) covering the last millennium and conducted with the MPI-ESM-P model. In particular, we describe an abrupt shift in the SPG strength detected in one ensemble member around 1600 that separates two different regimes of the North Atlantic climate, the later one featuring LIA-characteristics. We also explore the mechanism that renders the later regime stable for at least 250 years as well as the trigger of the SPG shift. Our results are supported by sensitivity experiments in which the initial conditions are slightly modified and/or the volcanic forcing removed. Results indicate that:

  1. 1.

    A millennial-scale weakening in the SPG strength throughout the past millennium is robustly simulated in the three Past1000s. This weakening can be relatively smooth, through a long-term trend, or abrupt, through a decadal-scale weakening shift.

  2. 2.

    The climate regime simulated after the SPG shift is stable on multicentennial time scales due to long-lasting anomalies within the North Atlantic: relatively fresh conditions in the upper Labrador Sea, weak oceanic meridional heat transport in the subpolar North Atlantic, expanded sea ice in the Labrador and Nordic seas, and persistent blocking-like SLP anomalies over the subpolar North Atlantic. Simulated post-shift climatic anomalies resemble the reconstructed changes in the area during the LIA, except for the timing of the occurrence. The simulated SPG shift is not associated with a shift in the AMOC strength or phase of the NAO. These results describe therefore a new possible mechanism for the onset of LIA-type episodes in the North Atlantic driven by SPG changes alone.

  3. 3.

    The SPG shift is triggered by anomalously large freshwater export from the Arctic and can mainly be attributed to internal climate variability. Volcanic forcing, also in the form of a cluster of small-magnitude volcanic eruptions, can facilitate the inception of the SPG shift, but is not a necessary condition. This highlights the relevance of internal climate variability with respect to externally forced changes in order to explain the evolution of the North Atlantic climate throughout the past millennium.