Abstract

This study evaluates different hypotheses of the origin of the Little Ice Age, focusing on the long-term response of Arctic sea ice and oceanic circulation to solar and volcanic perturbations. The authors analyze the Last Millennium Ensemble of climate model simulations carried out with the Community Earth System Model at the National Center for Atmospheric Research. The authors examine the duration and strength of volcanic perturbations, and the effects of initial and boundary conditions, such as the phase of the Atlantic multidecadal oscillation. They evaluate the impacts of these factors on decadal-to-multicentennial perturbations of the cryospheric, oceanic, and atmospheric components of the climate system. The authors show that, at least in the Last Millennium Ensemble, volcanic eruptions are followed by a decadal-scale positive response of the Atlantic multidecadal overturning circulation, followed by a centennial-scale enhancement of the Northern Hemispheric sea ice extent. It is hypothesized that a few mechanisms, not just one, may have to play a role in consistently explaining such a simulated climate response at both decadal and centennial time scales. The authors argue that large volcanic forcing is necessary to explain the origin and duration of Little Ice Age–like perturbations in the Last Millennium Ensemble. Other forcings might play a role as well. In particular, prolonged fluctuations in solar irradiance associated with solar minima potentially amplify the enhancement of the magnitude of volcanically triggered anomalies of Arctic sea ice extent.

1. Introduction

It is well known that large volcanic eruptions produce stratospheric sulfate aerosol clouds with a lifetime of 1–2 years and that these clouds scatter some of the incoming sunlight back to space, cooling Earth (Robock 2000). The aerosol layer circulates around the globe, warming the stratosphere through absorption of infrared radiation and resulting in a change to Earth’s energy balance. In addition, there are atmospheric and oceanic dynamical responses to large eruptions, producing characteristic regional and seasonal patterns of climatic perturbations. Although the role of volcanic eruptions as an important factor in the Little Ice Age (LIA) has been known for some time (Robock 1979; Crowley 2000), it has been suggested that a series of very large eruptions at the end of the thirteenth century, starting with the 1257 Samalas eruption (Lavigne et al. 2013), reduced North Atlantic oceanic heat flux into the Arctic so much that a feedback perpetuated this cool climate for centuries and started the LIA (Zhong et al. 2011; Miller et al. 2012), the 1250–1850 worldwide cooling of Earth’s climate. For example, in the climate model simulations of Zhong et al. (2011) and Miller et al. (2012), some cases of large volcanic eruptions have been followed by LIA-like anomalies. Additional long coupled climate simulations of this kind resulted in additional studies of the mechanisms by which this response could happen. For example, Lehner et al. (2013) found pronounced cooling of northern Europe to be associated with a nonlinear response to sea ice growth over the Nordic seas. Here, we take advantage of a new set of climate model simulations, the Last Millennium Ensemble (LME; Otto-Bliesner et al. 2016), to investigate whether centennial fluctuations of the climate, similar to events such as the historical LIA, arise in the simulated response of the coupled climate system, and in particular, how they depend on the external forcings that are prescribed to mimic historical activity of volcanic eruptions and solar minima during the last millennium.

The LME conducted a series of 1000-yr simulations with the National Center for Atmospheric Research (NCAR) Community Earth System Model version 1.1 (CESM; Hurrell et al. 2013), forced by changing insolation, volcanic eruptions, land surface, greenhouse gases, and aerosols, in combination and separately. We examined the results to see whether volcanic eruptions alone were enough to produce the LIA-like response and whether a combination of volcanic and solar forcings would significantly amplify the cooling.

Large eruptions can also produce decadal-scale shifts in the North Atlantic Ocean circulation (Otterå et al. 2010; Booth et al. 2012; Zanchettin et al. 2012, 2013). These oceanic changes are commonly studied as perturbations to the Atlantic meridional overturning circulation (AMOC; Buckley and Marshall 2016) and are often found to lead to sea surface temperature (SST) perturbations expressed statistically as the Atlantic multidecadal oscillation (AMO) and calculated by averaging SST over 20°–60°N. Swingedouw et al. (2015) showed that the response to volcanic eruptions depends on the phase of the AMO at the time of the eruption. Schleussner and Feulner (2013) suggested that an increase in Nordic sea ice extent on decadal time scales as a consequence of major eruptions would lead to a spinup of the subpolar gyre and a weakened AMOC, eventually causing a persistent, basinwide cooling. Here, we also examine the simulated response of the Atlantic meridional overturning circulation to large volcanic eruptions in the LME and how this might be linked to the long-term impacts on Arctic sea ice.

This paper is organized as follows. In section 2, we describe the model and dataset studied in this work and give a brief overview of its suitability for the subject of our study. In section 3, we define LIA-like analogs, and in particular, we examine Northern Hemisphere (NH) sea ice and its relation with North Atlantic circulation in the control run. Section 4 focuses on multidecadal fluctuations of the coupled climate system and examines their emergence in LME simulations with prescribed external forcings, focusing in particular on the role of volcanic eruptions. In section 5, we present analysis of centennial fluctuations of Northern Hemisphere climate and their relation to the LIA-like changes that follow strong volcanic eruptions. Section 6 provides plausible interpretations of our findings by comparing them with other modeling and observational studies. Section 7 discusses the relative role of other factors, in particular, solar minima, in triggering and sustaining LIA-like analogs. Section 8 provides the summary and conclusions.

2. The Last Millennium Ensemble

The climate model used for the LME, CESM-LME, consists of the Community Atmospheric Model, version 5 (Neale et al. 2012), the Los Alamos National Laboratory (LANL) Parallel Ocean Program, version 2 (Smith et al. 2010), the Community Land Model, version 4 (Lawrence et al. 2011), and the LANL Community Ice Code 4 (Hunke and Lipscomb 2008), with approximately 2° horizontal resolution in the atmosphere and land components and about 1° resolution in the ocean and sea ice components (Otto-Bliesner et al. 2016). The LME consists of 29 ensemble members, each 1000 years long, with random atmospheric perturbations imposed at the beginning of each simulation. Initial conditions correspond to the year 850 CE, with one simulation conducted as a control run (i.e., annual cycling of average incoming solar radiation with no external forcing). The remaining members include various combinations of time-dependent external forcings during the period 850–1850 CE. In particular, volcanic forcing and solar perturbations are imposed separately in five and four ensemble members, respectively. Moreover, prescribed greenhouse gases, orbital parameters, and land cover forcings are imposed separately in three ensemble members each. Another 10 ensemble members combine all of these external forcings at once. Volcanic and solar forcings are described briefly below and are characterized in more detail (along with the other external forcings) by Schmidt et al. (2011, 2012) and Otto-Bliesner et al. (2016).

Volcanic forcing is given by the monthly, zonally resolved distributions of stratospheric aerosols, as estimated by Gao et al. (2008) from ice cores. Approximately 70 eruptions occurred during the 1000-yr period of interest, varying in strength, duration, and latitude. The largest injection, 260 tons of SO2, occurred in 1257 and is associated with the eruption of the tropical Samalas volcano in the Rinjani volcanic complex (Lavigne et al. 2013). This event was followed by another three strong eruptions in 1268, 1275, and 1284. Other notable sulfate injections took place in 1452, 1600, and 1815, resulting from the tropical eruptions of Kuwae, Huaynaputina, and Tambora, respectively. Table 1 lists the 10 largest eruptions in this period, giving the actual dates as we now understand them and the dates that appeared in the forcing time series used for these simulations. The differences occur because of uncertainties in layer counting for the ice cores for the two oldest eruptions in the thirteenth century (difference of 1 year) and an error in the input dataset for Lakagigar (also known as Laki). In the rest of the paper, we refer to the eruptions by the dates actually used in the simulations.

Table 1.

The 10 largest volcanic eruptions during the 850–1850 LME simulations.

The 10 largest volcanic eruptions during the 850–1850 LME simulations.
The 10 largest volcanic eruptions during the 850–1850 LME simulations.

In the LME runs, the Gao et al. (2008) record of past volcanism was used with a linear relationship between volcanic emissions and radiative forcing. Some other simulations of the past millennium used the Crowley et al. (2008) volcanic record, in which a 2/3 power scaling was used for eruptions larger than the 1991 Pinatubo eruption (Crowley and Unterman 2013; Metzner et al. 2014). This scaling accounts for the growth of sulfate aerosols with larger emissions, which makes them less effective per unit mass than the constant smaller size distribution assumed with the linear relationship. Different events are captured in the two datasets, and common events often have different amplitudes based on the different methods of conversion (Sigl et al. 2014). As pointed out by Zambri et al. (2017), this means that the LME runs may have applied volcanic forcing that is larger than that being used now for such runs, which use the updated Sigl et al. (2015) dataset with scaling for large eruptions.

Solar forcing is given by spectrally resolved radiance changes matching the perturbations of total solar irradiance reconstructed by Vieira et al. (2011). This forcing is characterized by multidecadal perturbations that predominantly correspond to the Wolf (1280–1350), Sporer (1460–1550), and Maunder (1645–1715) sunspot minima. The corresponding amplitudes are on the order of 0.1 W m−2, in accordance with the calibration of climatic response against paleoclimatic proxies by Feulner (2011). Additionally, 11-yr perturbations are superimposed on the solar irradiance in an attempt to mimic the solar cycle, which is suppressed during solar minima and has an amplitude of up to 1 W m−2 otherwise.

As the Little Ice Age is the primary focus of this study, we are particularly fortunate that the LME was run with a climate model that has a good simulation of Arctic climate (Berdahl and Robock 2013). The NCAR Community Climate System Model, version 4 (CCSM4), used in the Paleoclimate Model Intercomparison Project Phase 3, a predecessor of CESM, was the only one of 10 models that was cold enough in the Arctic to simulate a sea ice climate like the observations. Moreover, the CESM1(CAM5) version of the model used for the LME simulations has been shown to provide a similarly good Arctic climate, albeit with sea ice sensitivity to the atmospheric temperatures seemingly enhanced (Kay et al. 2012; Holland and Landrum 2015). Additionally, CESM representation of the North Atlantic circulation has been shown to agree very well with observations (Danabasoglu et al. 2012; Yeager and Danabasoglu 2014). In particular, CCSM4’s AMOC exhibits low-frequency, nonperiodic fluctuations, with its frequency spectrum having a flat maximum within the range of 50–200 years that seems to agree well with the observations. Such AMOC characteristics agree much better with the observations than other models; some of them feature AMOC oscillating on much faster (i.e., 20 year) time scales. However, important biases remain in CESM, in particular a weaker-than-observed AMOC [attributed to a weaker North Atlantic Oscillation (NAO)], and this can impact the qualitative part of our analysis as well.

3. Sea ice cover as simulated in the coupled climate system

As LIA-like virtual analogs are of interest here, we first establish a suitable approach to define and study them in our dataset. No precise definition of the historical LIA exists; rather, this term refers loosely to the integrated knowledge (gained from numerous proxies gathered from all over the world; Mann et al. 2008; Ortega et al. 2015; Rahmstorf et al. 2015; Moreno-Chamarro et al. 2017) of anomalous cooling of multicentennial duration and global scale. Thus, to study analogs of the LIA, as simulated in the LME analyzed here, we focus primarily on the temporal evolution of total area of sea ice over the Arctic Ocean and adjacent oceans in northern latitudes. With LIA-like perturbations identified in this way, we subsequently provide their more detailed statistics. In particular, having established prolonged and positive anomalies of NH total sea ice area as our primary index of LIA-like analog here, we study its couplings (primarily by means of lead–lag regression and composites) with other parts of the simulated climate system. We also evaluate the concept that a sea ice–ocean feedback can explain the origin of the LIA by looking at two coupling mechanisms for Arctic sea ice–North Atlantic Ocean circulation that act on different time scales and have relative amplitudes that differ between internally driven and externally forced climate changes.

As shown in Fig. 1 for the control run, the correlation of NH sea ice area with poleward oceanic heat transfer, as averaged longitudinally over the North Atlantic sector, is systematically negative and can be as low as −0.6 for the latitude range of 60°–80°N. As shown in Fig. 2, spatially, these annual-average correlation coefficients are as high as 0.6 over the southern part of the Greenland Sea and as low as −0.6 in its northern part, as well as in the adjacent Norwegian Sea. These regions of strong correlation stretch southwestward along the eastern coast of Greenland, reaching as far as the Labrador Sea. These areas are known to be the oceanic corridors connecting the Atlantic basin with the Arctic Ocean (e.g., Lehner et al. 2013), where deep oceanic convection drives coupled system dynamics globally. Moreover, correlation of Arctic sea ice with poleward heat transfer extends south to the other regions of the Atlantic Ocean, south of 60°N and north of the equator, indicating coupling of the Atlantic meridional overturning circulation with Arctic climate. On the other hand, oceanic heat transfer between the Pacific and Arctic Oceans does not seem to play a significant role. Put together, this indicates, in agreement with other studies (e.g., Zhang 2015), that coupling of Arctic sea ice with poleward heat transfer in the Atlantic Ocean is pronounced in the internally driven climate of the control run.

Fig. 1.

Correlation coefficient (y axis) of zonally averaged, 10-yr smoothed, northward heat flux in the Atlantic Ocean and Arctic sea ice extent in the LME runs. Each line is one ensemble member forced as indicated.

Fig. 1.

Correlation coefficient (y axis) of zonally averaged, 10-yr smoothed, northward heat flux in the Atlantic Ocean and Arctic sea ice extent in the LME runs. Each line is one ensemble member forced as indicated.

Fig. 2.

Correlation coefficient of annual average northward heat flux for latitudes north of 50°N and Arctic sea ice extent for the LME control run.

Fig. 2.

Correlation coefficient of annual average northward heat flux for latitudes north of 50°N and Arctic sea ice extent for the LME control run.

Similar correlations of heat flux averaged over given latitudes and Arctic sea ice extent occur in the other LME members, albeit with a noticeable, systematic drop in magnitude (see Fig. 1). Here, the maximum decrease of correlation (amounting to 50% of the correlation coefficients) takes place for 0°–60°N and for the ensemble members with all the forcings prescribed simultaneously. This indicates that external perturbations impact the Arctic sea ice area (and, therefore, the LIA analogs studied here) in ways that do not involve, at least in the current version of CESM, sea ice–ocean feedbacks that operate in the same way as internal variabilities of the climate system. In particular, external forcing weakens the coupling of sea ice with poleward heat transfer within the North Atlantic (Danabasoglu et al. 2012; Kay et al. 2012; Holland and Landrum 2015). We revisit this idea later in the study.

As correlations of Arctic sea ice area and subpolar circulation of the Atlantic Ocean are relatively large for all LME members, here we look at them further in order to determine the more detailed picture of Arctic sea ice variability and, in particular, its relevant time scales. Figure 3 shows the lagged correlation for all ensemble members of LME of Arctic sea ice area with poleward heat transfer at 66°N, a representative latitude within the 60°–80°N range featuring strong correlations (shown in Fig. 1). The large negative correlation coefficient corresponds to simultaneous reduction of poleward oceanic heat transfer and expansion of sea ice area, a dependence that is found to be robust for different temporal filters. The minimum correlation occurs at lag zero, indicating a synchronous overturning circulation of the regional seas and thermodynamic and dynamical evolution of the surface waters, the latter having direct impact on the sea ice growth rate. Negative correlations are sustained for approximately 12 years in the case of the control run, indicating the multidecadal nature of internal variability of the coupled Arctic–North Atlantic climate. For simulations with external forcings, the average time scale lengthens up to 16 years. Importantly, a time scale of the same range is observed systematically for the entire range of latitudes less than 60°N, pointing to the planetary scale of oceanic circulation within the Atlantic basin (in CESM).

Fig. 3.

Correlation coefficient (y axis) of zonally averaged, annually averaged, northward heat flux in the Atlantic Ocean at 66°N and Arctic sea ice extent in the LME runs at different time lags. Each line is one ensemble member forced as indicated in Fig. 1. A negative lag means that heat flux leads sea ice. The black dash–dotted horizontal line shows the zero correlation level for reference. The red dash–dotted horizontal line shows the threshold for statistically significant negative correlation coefficients computed from a two-sided t test with N−2 degrees of freedom at significance level 0.01, where N = 1500 is equal to the number of decadal events in the LME dataset.

Fig. 3.

Correlation coefficient (y axis) of zonally averaged, annually averaged, northward heat flux in the Atlantic Ocean at 66°N and Arctic sea ice extent in the LME runs at different time lags. Each line is one ensemble member forced as indicated in Fig. 1. A negative lag means that heat flux leads sea ice. The black dash–dotted horizontal line shows the zero correlation level for reference. The red dash–dotted horizontal line shows the threshold for statistically significant negative correlation coefficients computed from a two-sided t test with N−2 degrees of freedom at significance level 0.01, where N = 1500 is equal to the number of decadal events in the LME dataset.

Given the strong coupling of Arctic climate with North Atlantic oceanic circulation in the LME coupled climate system, for the rest of our study, we apply an AMOC index whenever we want to study overturning circulation of the North Atlantic Ocean, particularly its response to volcanic eruptions or its association with centennial LIA-like anomalies. We define the AMOC index (shown in Fig. 4) as the Atlantic maximum overturning streamfunction, sampled at 35°N and subsequently low-pass filtered (10 years). Figure 4 shows that the results are not sensitive to the range of latitudes (20°–70°N) and that low-pass filtering is a suitable approach, given the dominance of the decadal time scale over higher frequencies.

Fig. 4.

AMOC index (maximum overturning streamfunction at any depth) for the LME control run for (a) annual average for 20°–70°N and for 35°N, (b) 13-yr running mean for 20°–70°N and for 35°N, (c) annual average for 20°–70°N, compared to 13-yr running mean for 35°N, and (d) annual average for 35°N, compared to 13-yr running mean for 35°N.

Fig. 4.

AMOC index (maximum overturning streamfunction at any depth) for the LME control run for (a) annual average for 20°–70°N and for 35°N, (b) 13-yr running mean for 20°–70°N and for 35°N, (c) annual average for 20°–70°N, compared to 13-yr running mean for 35°N, and (d) annual average for 35°N, compared to 13-yr running mean for 35°N.

In light of the LME statistics shown above and the decadal time scales at which Arctic sea ice is coupled with poleward heat transfer within the North Atlantic Ocean, here we focus on spatial patterns of SST evolution as obtained through lead–lag regression on the AMOC index. We focus on SST, as it is the part of the climate system located at the boundary of the two media of interest (sea ice and ocean) that determines the impact of overturning circulation of the North Atlantic Ocean on Arctic sea ice. As the upper layers of the ocean evolve with both oceanic and atmospheric influences, the particular relation of SST and AMOC remains unclear (Marini and Frankignoul 2014; Clement et al. 2015) and can differ significantly between models and observations. In the case of the LME, SST in both the Atlantic and Arctic basins appears strongly coupled with AMOC. To illustrate this, Fig. 5 shows 30-yr lead–lag regression patterns of SST against the AMOC index. The maximum value of AMOC (lag zero) is preceded by its slow increase on the multidecadal time scale, corresponding to intensification of the North Atlantic overturning circulation, with response of the North Atlantic and Arctic SST being particularly pronounced a few years later. In particular, Fig. 5 reveals that beginning at lag −5 years, SST starts to increase off the east coast of North America as well as over the Arctic Ocean. Anomalously warm waters found over the eastern Atlantic slowly move northward, leading to homogeneously warm waters for latitudes 50°–80°N over the next decade or so. The corresponding AMO index (not shown; defined as the mean SST over the North Atlantic latitudinal band of 20°–60°N) reaches a maximum with a few-year lag with respect to AMOC, with which it correlates as strongly as 0.6. In the meantime, positive Arctic SST anomalies intensify and propagate southward, accompanied by simultaneous melting of sea ice (not shown). A few years later, warm surface waters propagating north- and southward from the Atlantic and Arctic Oceans merge over subpolar latitudes, perturbing thermodynamic and dynamical conditions (i.e., buoyancy) of the water column there and subsequently suppressing deep convection taking place within regional seas. Over the next decade, SST and sea ice anomalies slowly return to climatological conditions. Continuous weakening of anomalously strong deep convection, as well as AMOC, leads to the opposite phase a few years later.

Fig. 5.

SST anomaly (K) composite for different phases of AMOC for the LME control run, with the SST lagging the AMOC peak by (a) −10, (b) −5, (c) 0, (d) 5, (e) 10, and (f) 15 years. Colors show SST anomalies of magnitude exceeding a threshold computed using a two-sided z test at significance level 0.01, with standard errors estimated using 50 degrees of freedom corresponding to the number of decadal events in the data.

Fig. 5.

SST anomaly (K) composite for different phases of AMOC for the LME control run, with the SST lagging the AMOC peak by (a) −10, (b) −5, (c) 0, (d) 5, (e) 10, and (f) 15 years. Colors show SST anomalies of magnitude exceeding a threshold computed using a two-sided z test at significance level 0.01, with standard errors estimated using 50 degrees of freedom corresponding to the number of decadal events in the data.

4. Multidecadal perturbations of climate and the response to volcanic eruptions

Internal fluctuations of Arctic sea ice and Atlantic circulation are significantly correlated, with positive anomalies of the former accompanied by negative ones of the latter. So far, we have discussed internal variability of Arctic sea ice and the North Atlantic Ocean operating on decadal time scales. Building on that, we now examine whether similar sea ice and ocean anomalies can last much longer (as in the case of LIA-like centennial expansion of Arctic sea ice) and if they can be triggered externally, in particular by volcanic eruptions.

For that, we analyze the five LME members with only volcanic forcing. The simulated response is expected to depend on the amplitude, timing, and location of the applied forcing, which vary significantly among the approximately 70 eruptions considered here. Additionally, for a given eruption, the decadal adjustment of the ocean (and, in general, the climate) can vary between different ensemble members, as volcanic eruptions impact predominantly atmospheric dynamics whose high-frequency chaotic behavior leads to distinctive adjustments of the rest of the climate system. Thus, we screen these five ensemble members to determine if the probability of a given AMOC phase increases after some eruptions, and if so, what the average time scale, phase, and amplitude of the preferred response are.

Figure 6 shows 20 years of the AMOC index as a function of time, with year zero corresponding to the moment of eruption for five ensemble members for the 10 most powerful eruptions (see Table 1). Before the eruptions, no AMOC phase can be identified as preferable among the 50 cases examined. Moreover, the AMOC amplitude evolves relatively slowly for each time series. A decade after the eruptions, however, a positive AMOC phase emerges slowly and steadily as the prevailing one. These results indicate that sufficiently strong eruptions induce a positive AMOC phase at decadal time scales irrespective of the initial phase. Such an AMOC response is in stark contrast to that reported by pioneering studies that looked at how volcanic eruptions could initiate the Little Ice Age and attributed it to an AMOC weakening (e.g., Zhong et al. 2011; Miller et al. 2012).

Fig. 6.

AMOC anomaly (Sv; 1 Sv ≡ 106 m3 s−1) following the 10 largest volcanic eruptions (see Table 1) for the five LME runs forced only with volcanic eruptions. Each line is for one eruption for one ensemble member. Each ensemble member is a different color. Year 0 is the year of the eruptions. The anomaly is calculated with respect to the 1000-yr mean for each run separately. While no phase dominates before the eruptions, a clear preference for positive phase is seen in the decade and a half following the eruptions.

Fig. 6.

AMOC anomaly (Sv; 1 Sv ≡ 106 m3 s−1) following the 10 largest volcanic eruptions (see Table 1) for the five LME runs forced only with volcanic eruptions. Each line is for one eruption for one ensemble member. Each ensemble member is a different color. Year 0 is the year of the eruptions. The anomaly is calculated with respect to the 1000-yr mean for each run separately. While no phase dominates before the eruptions, a clear preference for positive phase is seen in the decade and a half following the eruptions.

Figure 7 shows scatterplots of the AMOC phase at the eruption onset and 8 years after each eruption for all the eruptions in the volcano-only LME simulations. The relationship between the initial and final phases can be inferred either by examining the spread along both axes or through the gradient of the line of best fit (also shown in Fig. 7). The scatterplots support our previous conclusion: for the 15 strongest eruptions, relaxation toward the positive AMOC phase takes place irrespective of initial conditions (which have no sign preference, as indicated by the large and symmetric spread about zero). The next 15 strongest eruptions also result in a positive AMOC phase, although, as indicated by the gradient of the fitted line, with smaller amplitude and thus lower probability. The remaining eruptions are much weaker and do not have a significant impact on AMOC. The absence of a pronounced response is indicated by the symmetric scatterplot spread around zero and the quasi-straight fitted curve overlaying the individual observations. This indicates a normal oscillation with a 10–20-yr periodicity, so that if AMOC was in one phase at the time of the eruptions, it would likely be in the opposite phase 8 years later. In summary, our results indicate that volcanic eruptions can lead to a strengthening of AMOC on decadal time scales (and, presumably, reduction of Arctic sea ice, if corresponding dependencies found for the control run dominate as well in the volcanic simulations), but the odds of that happening increase with the strength of the given eruption. However, as we show later, this tendency is actually overwhelmed by the direct radiative forcing from large eruptions, resulting in the opposite response: an expansion of Arctic sea ice.

Fig. 7.

Scatterplot of the AMOC anomaly (Sv) 8 years after each eruption on the y axis, compared to the AMOC anomaly in the year of the eruption on the x axis, for 69 volcanic eruptions for the five LME runs forced only with volcanic eruptions. The anomaly is calculated with respect to the 1000-yr mean for each run separately for (top left) the 15 largest eruptions, (top right) the next 15 largest, (bottom left) the next 15 largest, and (bottom right) the 24 smallest eruptions. The red, black, pink, and blue stars are for each eruption in each run, and the green stars correspond to the best nonlinear fit (polynomial, of degree 3) for all the points.

Fig. 7.

Scatterplot of the AMOC anomaly (Sv) 8 years after each eruption on the y axis, compared to the AMOC anomaly in the year of the eruption on the x axis, for 69 volcanic eruptions for the five LME runs forced only with volcanic eruptions. The anomaly is calculated with respect to the 1000-yr mean for each run separately for (top left) the 15 largest eruptions, (top right) the next 15 largest, (bottom left) the next 15 largest, and (bottom right) the 24 smallest eruptions. The red, black, pink, and blue stars are for each eruption in each run, and the green stars correspond to the best nonlinear fit (polynomial, of degree 3) for all the points.

The existence (and, in particular, emergence or termination) of the multicentennial cooling and sea ice extent has been inferred from various proxies (e.g., Dobrowolski et al. 2002, 2016; Axford et al. 2011; Larsen et al. 2012; Gałka et al. 2014, 2016; Marcisz et al. 2015; Feurdean et al. 2015; Lamentowicz et al. 2015) and is typically associated with three specific eruptions, namely, those of 1258 (along with three succeeding eruptions), 1452, and 1600. To take a more detailed look at these periods, as well as the decades following other strong eruptions, we plot in Fig. 8 the AMOC time series for the volcanic ensemble members, separately for the five strongest eruptions included in Fig. 6. In the case of the 1258 Samalas and 1452 Kuwae eruptions, no AMOC phase appears to prevail in the respective year of the eruptions, yet for all five ensemble members, AMOC strengthens during the following decade. Similarly, mean background conditions for the Huaynaputina eruption do not appear anomalous, but a positive AMOC phase emerges in the aftermath in four out of five ensemble members. In the case of the Tambora eruption, all ensemble members are characterized by a negative AMOC phase before or during year 1815 and gradually shift toward a positive phase over the 15 years that follow.

Fig. 8.

AMOC anomaly (Sv) following the five largest volcanic eruptions for the five LME runs forced only with volcanic eruptions. Year 0 is the year of the eruptions. The anomaly is calculated with respect to the 1000-yr mean for each run separately. In general, the AMOC anomaly becomes positive no matter what the initial AMOC state.

Fig. 8.

AMOC anomaly (Sv) following the five largest volcanic eruptions for the five LME runs forced only with volcanic eruptions. Year 0 is the year of the eruptions. The anomaly is calculated with respect to the 1000-yr mean for each run separately. In general, the AMOC anomaly becomes positive no matter what the initial AMOC state.

5. Centennial perturbations of Arctic sea ice

We have shown so far that in the case of the last millennium climate simulated by CESM (and also other studies, e.g., Zhang 2015), perturbations of Arctic sea ice area are coupled strongly to decadal anomalies of North Atlantic SST, as they are strongly tied to AMOC. This variability acts internally within the coupled atmosphere–ocean climate system on the time scale of approximately 20 years, as demonstrated for the control run in the preceding sections. At the same time, interdecadal fluctuations of AMOC (and, thus, Arctic sea ice and the rest of the climate system) are also forced externally, as shown above for volcanic eruptions prescribed in the LME. Here, we follow on that finding and examine longer, centennial time scales. In particular, Fig. 9 shows millennial time series of perturbations of AMOC for runs with volcanic forcing and of total Arctic sea ice area for runs with volcanic, solar, and all forcings. For the volcanic runs, AMOC anomalies (Fig. 9b) are of similar magnitude to those in the control run (not shown). Moreover, these externally forced perturbations exhibit fluctuations at multidecadal time scales (a time span shorter than the one observed for Arctic sea ice extent discussed below). The centennial response of this coupled climate system to volcanic eruptions is tied to decadal oceanic fluctuations in a complex way, as already observed by Pausata et al. (2015). In particular, after reaching the peak intensity a few decades after the strongest eruption, prolonged weakening of AMOC begins, leading to a negative phase of AMOC of centennial extent, with pronounced multidecadal fluctuations of smaller amplitude occurring in the meantime. For Arctic sea ice, at the same time, significant expansion can be observed systematically for simulations with volcanic eruptions (see Fig. 9c). The maximum sea ice extent occurs abruptly a few years following each large eruption and, thus, cannot be attributed simply to coupling with oceanic AMOC anomalies. During the decade that follows, a slow decrease of Arctic sea ice (along with AMOC strengthening) is observed. This prolonged reduction of Arctic sea ice does not lead, however, as in the case of AMOC, to its opposite phase of centennial duration. Instead, relatively small but positive anomalies of Arctic sea ice prevail in simulations with volcanic forcing. Such a centennial regime of moderately enhanced Arctic sea ice extent and weakened North Atlantic oceanic circulation occurs after the Samalas, Kuwae, and Huaynaputina eruptions, following decadal adjustment of accelerated AMOC and abruptly expanded Arctic sea ice and entering possibly a new quasi-equilibrium state of the coupled sea ice–ocean–atmosphere system.

Fig. 9.

(a) Volcanic forcing (global average of volcanic aerosol mass in kg m−2). (b) Anomaly of the AMOC index (Sv) for the volcano-only runs. (c) Total Arctic sea ice area (×106 km2, smoothed with a 10-yr running mean). (d) As in (c), but for solar forcing only. (e) As for (c), but for all forcings. (f) As in (e), but for a linear combination of the runs with each separate forcing, including the forcings not plotted separately here, land cover, orbital, and greenhouse gases. Colored lines are for each of the ensemble members, and the black line is their mean. The horizontal dashed line in (c)–(f) is the mean from the control run.

Fig. 9.

(a) Volcanic forcing (global average of volcanic aerosol mass in kg m−2). (b) Anomaly of the AMOC index (Sv) for the volcano-only runs. (c) Total Arctic sea ice area (×106 km2, smoothed with a 10-yr running mean). (d) As in (c), but for solar forcing only. (e) As for (c), but for all forcings. (f) As in (e), but for a linear combination of the runs with each separate forcing, including the forcings not plotted separately here, land cover, orbital, and greenhouse gases. Colored lines are for each of the ensemble members, and the black line is their mean. The horizontal dashed line in (c)–(f) is the mean from the control run.

In summary, volcanically triggered centennial perturbations confirm that external forcings can, indeed, lead to a new regime for both AMOC and Arctic sea ice extent and that the details of this transition are more complex than the ones found in the case of internal variability and inferred from large lead–lag correlations of Arctic sea ice and Atlantic poleward heat transfer in the control run. The next section follows up with further discussion of this disparity between external and internal drivers of LIA-like perturbations and their occurrence in the Last Millennium Ensemble. More systematic analysis of this regime emergence, intensity and duration, and sensitivity to various factors (e.g., strength of volcanic and other forcings) would require an additional set of targeted climate simulations, which we plan to conduct and report in a future study.

Centennial trends show that the Arctic sea ice time series in the volcanic runs can be categorized into three main periods: the first four centuries (Warm Medieval Anomaly, 850–1250), characterized by negative anomalies of total Arctic sea ice; the following 350 years that correspond to the initial period of LIA-like regime of the LIA and span between the Samalas and Huaynaputina eruptions (1258–1600); and the two centuries that follow (1600–1800), with Arctic sea ice anomalies even more prominent. Moreover, the occurrence of LIA-like anomalies (and, presumably, other regional phenomena not discussed here) appears to emerge after a multidecadal period of AMOC adjustment, with volcanically triggered perturbations at these two time scales having distinct dynamics, as demonstrated above.

In summary, our findings indicate that although Arctic sea ice fluctuates along with AMOC on interdecadal time scales in the LME control run, externally forced perturbations of sufficiently strong amplitude can lead to more complicated climate response on the centennial time scale. We illustrate this by showing Arctic sea ice and northward heat transfer averaged over a 2-century time span during the periods that overlap in time with historical events of the Medieval Warm Period (MWP, 1050–1250) and the LIA (1615–1815) in Fig. 10. The Arctic sea ice is, as expected, negative for the first period and positive for the second period—a tendency systematically observed in case of simulations with volcanic forcings and even more pronounced in case of the all-forcing runs. At the same time, the AMOC strength does not differ, on average, for the MWP and LIA, except for large positive values for some of the all-forcing runs (as opposed to negative values that would be expected in the case of sea ice–ocean feedback dominance). In general, the AMOC amplitude is characterized by a large spread, with small differences among the periods that are typically associated with the historical events of the LIA and MWP. Quantitatively, correlation of AMOC and Arctic sea ice drops from the decadal time scales analyzed above to the centennial time scales analyzed here. Given that, one can infer that for climates driven by external forcings of the kind, strength, and frequency known from the past 1000 years of Earth’s history, parts of the climate system (e.g., the atmospheric or oceanic circulation in subpolar and polar latitudes) other than the North Atlantic Ocean (important in the case of internally driven dynamics) should be considered while studying LIA-like events.

Fig. 10.

Scatterplot of average Arctic sea ice area anomaly (×106 km2) and poleward oceanic circulation anomaly at 66°N (Sv), for runs with all the forcings (black, 10 runs), volcanic forcing (red, 5 runs), and solar forcing (green, 4 runs), during a portion of the LIA (1615–1815, marked with *) and a portion of the MWP (1050–1250, marked with +). The anomalies are with respect to the mean over the entire period for the given ensemble member. Both sea ice and heat anomalies are normalized by dividing by the absolute value of the maximum anomaly from all the ensemble members for the LIA period.

Fig. 10.

Scatterplot of average Arctic sea ice area anomaly (×106 km2) and poleward oceanic circulation anomaly at 66°N (Sv), for runs with all the forcings (black, 10 runs), volcanic forcing (red, 5 runs), and solar forcing (green, 4 runs), during a portion of the LIA (1615–1815, marked with *) and a portion of the MWP (1050–1250, marked with +). The anomalies are with respect to the mean over the entire period for the given ensemble member. Both sea ice and heat anomalies are normalized by dividing by the absolute value of the maximum anomaly from all the ensemble members for the LIA period.

6. Discussion of volcanic eruptions and their impact on LIA-like analogs

In this section, we proceed further with discussion of our LME statistics and the LIA-like centennial anomalies extracted earlier, with a particular focus on providing some relevant context based on other studies and their interpretation of LIA physics. Importantly, although centennial perturbations are our primary interest here, other volcanically triggered temporal scales have been shown to act at the same time, indicating that cross-scale temporal interactions could be relevant for a more complete explanation of LIA-like phenomena.

The prevailing hypothesis regarding the cause of the Little Ice Age that has been widely embraced by the climate community in recent years (Zhong et al. 2011; Miller et al. 2012; Lehner et al. 2013; and others) proposes a physical mechanism driven primarily by a positive feedback of sea ice and the ocean. This feedback relies heavily on Arctic sea ice coupling with poleward heat transfer in the North Atlantic Ocean—a feature that we robustly extracted for the internally evolving climate of the control run (see Fig. 1 and its discussion), in agreement with many other studies (Holland et al. 2006; Hwang et al. 2011; Zhong et al. 2011; Miller et al. 2012; Sandø et al. 2014). According to these studies, the historical event of the Little Ice Age was caused by strong forcings that acted externally on Earth’s internal climate in a way that mainly involved the sea ice–ocean feedback. Moreover, the physical mechanism that explains the LIA relies predominantly on simple amplification of perturbations of the parts of the climate system relevant to this feedback. Due to the timely occurrence of a prolonged solar minimum and powerful volcanic eruptions, these two forcings are given as the most probable causes of the LIA. In particular, abrupt atmospheric cooling imposed externally either by volcanic or solar forcing onto the region where the sea ice–ocean feedback operates would introduce a significant amount of Arctic sea ice. Perturbed in this way, the climate dynamics adjust accordingly by increasing the southward export of sea ice and the subsequent thermodynamic conditions (i.e., buoyancy driven by salinity and temperature) of waters merging in the Labrador Sea. The resulting reduction of deep oceanic convection subsequently weakens AMOC, decreases poleward heat transfer, cools the Arctic Ocean, and increases sea ice anomalies there.

In Fig. 9 and its subsequent discussion, we also find multidecadal periods of simultaneously enhanced Arctic sea ice extent and weakening of the oceanic circulation within the North Atlantic basin. Thus, our results appear consistent with the above discussion of the mechanisms proposed by other studies and with the hypothesis of LIA-like perturbations being sustained on centennial time scales through sea ice–ocean feedback. In particular, three LIA-like prolonged periods—1320–75, 1520–1610, and 1640–1740—followed the powerful volcanic eruptions of Samalas (1258), Kuwae (1452), and Huaynaputina (1600). However, a weakening of AMOC occurs with a multidecadal lag with respect to each eruption, while for the intervening few decades, AMOC is strengthened instead. This feature is observed robustly (given the limited ensemble analyzed here) yet has not been observed often in other studies, where a fairly steady decrease in poleward heat transfer has taken place within a decade after an eruption. Such a discrepancy could potentially have arisen because the AMOC strength is regulated by different factors on shorter, decadal time scales, after which an equilibrated AMOC and Arctic sea ice system emerges, with the sea ice–ocean feedback becoming the prevailing mechanism by then. Moreover, the AMOC transition, from its initially accelerated stage to its negative phase, is somewhat expected due to its deterministic and highly predictable nature established by other studies (Hurrell et al. 2006; Msadek et al. 2010; Tulloch and Marshall 2012; Branstator and Teng 2014). In particular, AMOC fluctuations appear to be fairly coherent (i.e., of the same sign and duration) in different LME members analyzed here for up to a century.

The AMOC strengthening that directly follows volcanic eruptions has not been robustly observed among other studies—on the contrary, its weakening has been reported more widely (see discussion above)—and thus, no robust understanding of it has yet been achieved by the climate community. To understand the physical causes of this particular transition of AMOC, the subpolar waters of the Labrador Sea appear to be a good candidate for study, as deep convection there has long been recognized as a crucial regulator of deep oceanic circulation within the North Atlantic basin. Such geographical preference is determined by the low temperatures, high salinity, and relatively high buoyancy of the upper-ocean layers there, compared with other regions. Various processes have been proposed as important drivers of this region’s oceanic conditions, some of them favorable to development and intensification of deep oceanic convection and subsequent AMOC strengthening. Atmospheric perturbations associated with high-frequency weather-related events are one well-known factor having strong impact on surface conditions of the Labrador Sea, with strong surface winds and subsequent invigoration in evaporation from the ocean surface leading to its substantially lower temperatures and higher salinity. The origins of some significant anomalies of negative salinity and subsequent shutdown of deep oceanic convection, such as the famous 1969–70 event observed in the Labrador Sea and referred to as the Great Salinity Anomaly (Dickson et al. 1988; Gelderloos et al. 2012), have been partially attributed to atmospheric forcing; an anomalously weak NAO and enhanced westerlies over the Labrador Sea resulted in a net reduction of surface heat loss, intense freshening, and a collapse of deep ocean convection. Our results, however, do not find support for this mechanism, as composites of surface air pressure do not reveal any consistent pattern of atmospheric circulation among ensemble members for this region at decadal time scales following volcanic eruptions. Similarly, another possible factor—anomalously positive influx of Arctic waters—is also one that can be excluded, as it works in the opposite direction (see the discussion above of sea ice–ocean feedbacks) and leads to AMOC reduction due to freshening of surface water and reduction of deep convection by increased export of sea ice from the Arctic Sea. Thus, another possible mechanism of adjusting the thermodynamic conditions within the Labrador Sea—the advection of waters carried by the Irminger Current, which consists of the northwest branch of the subpolar gyre—is also considered here (Hátún et al. 2005; Böning et al. 2006; Born et al. 2016). Here, indeed, we find systematic evidence of oceanic circulation anomalies that overlap in time with the AMOC strengthening and with the decades following the relevant eruptions. In particular, Fig. 11 shows the subpolar gyre strengthening as a distinct response of the oceanic circulation, with AMOC-driven AMO anomalies being negligible in the decade that directly followed the eruption. Moreover, Fig. 12d indicates prolonged (on decadal time scales) reduction of the barotropic streamfunction in the Labrador Sea that, in turn, corresponds to a strengthening of the subpolar gyre.

Fig. 11.

NH (a) diagnostic barotropic streamfunction (Sv) and (b) sea ice extent (%) anomalies, as compared to the control run. Composites result from averaging over the simulation period 1260–80 and over five LME simulations with just volcanic forcing. Colors show anomalies of magnitude exceeding a threshold computed using a two-sided z test at significance level 0.05, with standard errors estimated using 5 degrees of freedom, corresponding to the number of decadal events in the data.

Fig. 11.

NH (a) diagnostic barotropic streamfunction (Sv) and (b) sea ice extent (%) anomalies, as compared to the control run. Composites result from averaging over the simulation period 1260–80 and over five LME simulations with just volcanic forcing. Colors show anomalies of magnitude exceeding a threshold computed using a two-sided z test at significance level 0.05, with standard errors estimated using 5 degrees of freedom, corresponding to the number of decadal events in the data.

Fig. 12.

(a) Sea level pressure (hPa) difference anomalies for North Atlantic Ocean sectors (between the 60°–75°N, 30°–15°W region and the 40°–55°N, 25°–10°W region); (b) total precipitation anomalies (mm day−1) for the western Atlantic Ocean (averaged over 30°–45°N and 75°–45°W); (c) surface salinity anomalies (PSU) for mid-Atlantic sector (mean over 45°–60°N and 30°–15°W); and (d) barotropic streamfunction (Sv) anomalies for the Labrador Sea (averaged over 45°–65°N and 60°–30°W). Colored lines correspond to time series associated with running averages (annual means) of five LME simulations with just volcanic forcing. The black line is their mean. The Samalas eruption occurred at time 0. Anomalies are with respect to the mean for each simulation.

Fig. 12.

(a) Sea level pressure (hPa) difference anomalies for North Atlantic Ocean sectors (between the 60°–75°N, 30°–15°W region and the 40°–55°N, 25°–10°W region); (b) total precipitation anomalies (mm day−1) for the western Atlantic Ocean (averaged over 30°–45°N and 75°–45°W); (c) surface salinity anomalies (PSU) for mid-Atlantic sector (mean over 45°–60°N and 30°–15°W); and (d) barotropic streamfunction (Sv) anomalies for the Labrador Sea (averaged over 45°–65°N and 60°–30°W). Colored lines correspond to time series associated with running averages (annual means) of five LME simulations with just volcanic forcing. The black line is their mean. The Samalas eruption occurred at time 0. Anomalies are with respect to the mean for each simulation.

Many possible drivers of the subpolar gyre variability have been proposed to play a role by other studies. Besides the obvious one, which is AMOC itself, atmospheric (i.e., wind driven) forcings are often proposed (Trouet et al. 2009, 2012; Moffa-Sánchez 2014a,b), either in direct association with local conditions (e.g., jet stream instabilities, atmospheric blocking, and NAO) or remote teleconnections (e.g., from the Pacific; Gray et al. 2010). Again, our findings do not support this hypothesis, as we found no low-frequency anomalies of the atmospheric circulation over the subpolar regions during this period. In particular, we do not detect any decadal shifts for various atmospheric fields (e.g., western Atlantic precipitation or NAO variability; shown in Figs. 12a,b). The response of atmospheric dynamics is dominated by interannual time scales. In particular, years 2 to 4 seem to be characterized by intensification of western North Atlantic precipitation, which can be associated with El Niño and subsequent reduction of precipitation due to La Niña; both are known to have a systematic signal in the Last Millennium Ensemble (Otto-Bliesner et al. 2016), with teleconnections from the tropical Pacific to the Atlantic basin. Model results, and, in particular, a relatively weak NAO signal, seem to be inconsistent with recent paleo reconstructions (Ortega et al. 2015) and may be potentially attributed to CESM model biases, associated with poorly resolved stratosphere and its projection on couplings with the tropospheric dynamics and NAO signal.

Our findings lead us to argue that the intensification of the subpolar gyre is regulated, at least in the LME, by Gulf Stream perturbations introduced over the mid-Atlantic (see Fig. 11), which, in turn, can be traced back to volcanically induced cooling of the atmosphere, a reduction in precipitation, and an increase in salinity of oceanic currents. Additional statistics on Samalas, shown in Fig. 12, seem to support our hypothesis. In particular, the volcanic eruption is, indeed, followed by oceanic fluctuations (shown in Figs. 12c,d for salinity and the subpolar gyre) on decadal time scales. At the same time, such slow perturbations do not have their atmospheric counterparts (as demonstrated for pressure and precipitation anomalies over the Atlantic Ocean). However, volcanically triggered atmospheric dynamics perturb North Atlantic Ocean conditions strongly enough for them to shift into a positive AMOC. In particular, pronounced cooling in the year of the eruption is associated with global reduction in the precipitation (shown in Fig. 12b for the western Atlantic Ocean). This results in increased salinity of the central Atlantic Ocean a couple of years later (years 2–3), subsequent reduction in the barotropic streamfunction (years 6–9), and potential kickoff of the AMOC transition to its positive phase.

Our idea is not entirely new and is supported by several independent arguments put forward elsewhere. Recent studies argue that anthropogenically driven climate change and anomalous freshening of the North Atlantic, originating, for example, in a modified influx of precipitation and river discharge to the Gulf Stream or the melting of Greenland, can potentially impact AMOC (Zhang et al. 2011, 2014; Rahmstorf et al. 2015). Also, although solid evidence of anomalous cooling and sea ice extent exists for the historical period associated with the LIA (e.g., Christiansen and Ljungqvist 2011, 2012; Miller et al. 2012; and others), paleo proxies of oceanic circulation are spatially sparser and also harder to interpret (Kurahashi-Nakamura et al. 2014; Zhang et al. 2015; Christiansen and Ljungqvist 2017). Interestingly, contrary to common knowledge based partially on a number of proxy-related studies, some recent paleo and modeling studies do not find support of AMOC weakening during the LIA (Rahmstorf et al. 2015; Moreno-Chamarro et al. 2017) and, thus, provide a challenge to the proxy community as well as to the theories of LIA origin and the particular role of sea ice–ocean feedback in it (Wanamaker et al. 2012; Dylmer et al. 2013). Moreover, some modeling studies reach similar conclusions to ours; for example, Pausata et al. (2015), who studied the long-term impact of high-latitude eruptions on AMOC in 10 simulations conducted with the Norwegian Earth System Model NorESM1-M, report similar strengthening over the first few decades and gradual weakening over the next few decades following the imposed eruptions. Midlatitude freshening of the Atlantic, a reduction in precipitation, and enhancement of salinity content within the Gulf Stream are also hypothesized there as possible scenarios leading to LIA-like events.

Important questions remain: for example, what determines the opposite AMOC response that is found to arise among different models? To gain more quantitative insight into this question, the relevance of remote environmental conditions for deep convection in the Labrador Sea should be studied in a more systematic fashion, with a particular recommendation of sensitivity tests of factors such as the following: (i) the parameterization of oceanic and atmospheric convection, (ii) sea ice physics, (iii) its coupling with the ocean component of the system, and (iv) the impact of ocean resolution, in particular its topography, which serves as a bottom boundary condition. These types of questions have been raised before, for example, in the context of various limitations of models in simulating a realistic AMOC (Weaver et al. 2012; Buckley and Marshall 2016). In particular, Danabasoglu et al. (2012), Yeager and Danabasoglu (2012), and Yeager (2015) found that simulated AMOC properties differ strongly in different versions of the CESM model and attributed this to a range of possible factors. To name a few, they emphasized the role of subgrid-scale parameterization of oceanic eddies on density anomalies, parameterization of Nordic sea overflows, stratification of the Labrador Sea, and bottom topography.

Another interesting observation can be made based on the spatial map of NH sea ice anomalies shown in Fig. 11. As expected, short-lived (i.e., of interannual time scale) yet strong atmospheric cooling following the eruptions impacts Arctic sea ice directly and leads to its expansion (which persists for decades). Interestingly, though, the magnitude of this expansion exhibits pronounced regional differences, with the Pacific Ocean seemingly playing a large role. This is somewhat surprising and contradictory to the main perception that has been brought up earlier in regard to our LME statistics in the case of the control run and also other studies on LIA origins. According to this view, Arctic sea ice area is regulated primarily by the poleward transfer of heat in the Atlantic Ocean, with the Pacific Ocean having negligible impact. The Pacific Ocean’s importance in regulating the total amount of Arctic sea ice is reported less often, with the particular role of external forcings remaining relatively poorly understood. Besides differences in the thermodynamic properties of the atmosphere and ocean, and subsequently, the efficiency in sea ice formation, other factors could potentially play an indirect role, for example, by impacting the poorly understood oceanic circulation of the Bering Sea. For example, freshwater anomalies introduced either by anomalous precipitation over the region or indirectly through river discharge could be of significance. A targeted set of sensitivity tests would ultimately be needed to validate this or other potential factors, as the robustness of our results is inherently limited by the small sample of simulations employed here. Irrespectively of the nature of these regional differences, the predominant amplification of sea ice over the Bering Sea gives some important clues for interpretation of Fig. 1, and, in particular, our earlier remark that external forcings reduce significantly (up to half, in the case of simulations with all the forcings) correlation of Arctic sea ice with poleward heat transfer in the Atlantic Ocean.

To summarize, the first decade of abrupt advancement in Arctic sea ice extent and the more gradual AMOC strengthening tend to lead, over the course of many decades, to a new regime, characterized predominantly by negative AMOC phase, with Arctic sea ice anomalies remaining moderately positive. Weakening of AMOC and simultaneously enhanced Arctic sea ice is in broad agreement with previously discussed coupling mechanisms of sea ice–ocean feedback. Furthermore, while AMOC weakens, the Pacific Ocean preserves a positive sea ice anomaly. This can be explained by negative Pacific decadal oscillation anomalies that are triggered, potentially, by teleconnections from the Atlantic (Zhang and Delworth 2007; Okumura et al. 2009; Zhang and Zhao 2015) and its coupling with sea ice and atmospheric circulation over the northern seas of the Pacific Ocean (Cavalieri and Parkinson 1987; Rodionov et al. 2007; Danielson et al. 2011; Kaufman et al. 2011; Wendler et al. 2014; Harada et al. 2014).

7. Role of solar forcing

We have demonstrated how volcanic eruptions can perturb the climate for up to centennial periods. In particular, Arctic sea ice anomalies exhibit a consistently abrupt and significant increase over the first decade after an eruption that drops to a moderate level and remains quasi-equilibrated for many decades afterward. The oceanic circulation in the North Atlantic, in turn, first exhibits multidecadal strengthening, after which it transits to an anomalously weak regime. Volcanic eruptions appear to be the crucial component in inducing these perturbations, and the strongest of these eruptions consistently precede LIA anomalies that last for up to a century, in agreement with the historical record. It must be acknowledged, though, that other external forcings have also been shown to be of importance and, thus, the question of how perturbations introduced by them compare to those from just volcanic forcings deserves further consideration. Centennial-scale Arctic sea ice expansion has been shown to emerge in simulations with just volcanic forcing, and it appears that this behavior can be detected even more clearly in simulations with all forcings (see Fig. 9e), but not with just solar forcing (Fig. 9d). At the same time, when combining the Arctic sea ice responses from the mean of each of the sets of runs done with individual forcings (Fig. 9f) and comparing them to the results from all the forcings applied at once (Fig. 9e), the linear combination shows larger responses. This indicates that to explain the simulated amplitude for all forcings combined, nonlinear interactions are required in a coupled atmosphere–ocean–sea ice North Atlantic–Arctic system. In particular, as our results indicate, and as others have shown as well, solar forcing can potentially impact the emergence and sustainability of LIA-like regimes.

While it appears from Fig. 10 that the response to solar forcing is as large as that to volcanic forcing, those plots show averages of 200-yr periods, and comparing Figs. 9c–e, it is clear that the volcanic-only runs produce a jump into an LIA-like perturbation at the end of the thirteenth century that is not evident for the solar runs. The all-forcing runs also exhibit a pronounced spike, which can be significantly larger than in simulations with solar forcings, as in the case of Arctic sea ice for the 1600 eruption and the following few years. Also, Fig. 13 shows time series of the volcanic and solar forcing expressed in terms of changes of the global average top-of-atmosphere shortwave radiative flux. Clearly, the volcanic forcing dominates, in agreement with other studies, although there happens to be a small, negative solar forcing coincident with the three largest volcanic eruptions (e.g., Bauer et al. 2003). So while it seems that volcanic eruptions can trigger the LIA, why, in the case of centennial time scales, do the all-forcing runs produce a more robust LIA than the volcanic forcing runs only, and what is the role of solar forcing?

Fig. 13.

Volcanic and solar top of atmosphere radiative forcing (W m−2) applied during the LME runs: (a) annual average volcanic forcing, (b) annual average solar forcing, and (c) sum of volcanic and solar forcing. Orange curves are 40-yr running mean and use the y axis on the right.

Fig. 13.

Volcanic and solar top of atmosphere radiative forcing (W m−2) applied during the LME runs: (a) annual average volcanic forcing, (b) annual average solar forcing, and (c) sum of volcanic and solar forcing. Orange curves are 40-yr running mean and use the y axis on the right.

The solar forcing applied in the LME (Fig. 13) indicates that the first and third solar minima, known as the Wolf and Maunder Minima, occurred decades after the Samalas and Huaynaputina eruptions and lasted for many decades. Climate exposure to these prolonged reductions in incoming radiation, therefore, appears to occur with a multidecadal lag with respect to volcanically induced periods of sea ice expansion during periods when negative AMOC prevails already with ongoing expansion of Arctic sea ice. The climate response for one of these periods can be inferred from Fig. 14, which shows a spatial map of Arctic sea ice anomalies and NH sea level pressure anomalies as obtained for simulations with just solar forcing for a 50-yr period overlapping with the Maunder Minimum. Atmospheric cooling results in the expansion of Arctic sea ice, being particularly pronounced over the Nordic seas. The atmospheric dynamics also shift systematically across different ensemble members, an outcome not found for simulations with volcanic eruptions, with a positive anomaly of surface pressure over the Arctic and Atlantic sectors typically associated with NAO variability. An adjustment of atmospheric circulation follows, with anomalously anticyclonic circulation over the Nordic seas. Such a pattern tends to weaken the jet stream and the strength of midlatitude storm tracks, favoring in turn anomalous advection of moist air toward subpolar regions, a reduction of AMOC strength, sea ice migration southward, and anomalous sea ice growth in the Nordic seas. The positive correlation of the solar minima with surface pressure at polar latitudes, along with the given atmospheric circulation response, seems to agree with arguments provided by such simple thermal wind balance. The temporal evolution of atmospheric and oceanic indices, shown in Fig. 15, seems to support our hypothesis. In particular, one can observe that the Maunder Minimum subperiod associated with the strongest reduction of incoming solar radiation (i.e., years 1670–1700) is systematically associated with a negative NAO phase and weakened subpolar gyre (and, as discussed earlier, reduced AMOC). In summary, we suggest that solar minima following volcanic eruptions enhance Arctic sea ice significantly, at least in the LME; that a multidecadal lag of the former to the latter can, indeed, play a role; and that in agreement with other studies, changes in atmospheric circulation patterns are important in that process (Moffa-Sánchez et al. 2014a,b; Kobashi et al. 2015). In the future, we hope to validate our findings by performing a set of well-targeted simulations that would improve further the current constraints and statistics of the LME, focusing in particular on the mutual timing and strength of solar and volcanic forcings.

Fig. 14.

NH (a) surface pressure (hPa) and (b) sea ice extent (%) anomalies, as compared to the control run. Composites result from averaging over simulated years 1650–1700 and over four LME simulations with just solar forcing. Colors show anomalies of magnitude exceeding a threshold computed using a two-sided z test at significance level 0.05, with standard errors estimated using 4 degrees of freedom, corresponding to the number of decadal events in the data.

Fig. 14.

NH (a) surface pressure (hPa) and (b) sea ice extent (%) anomalies, as compared to the control run. Composites result from averaging over simulated years 1650–1700 and over four LME simulations with just solar forcing. Colors show anomalies of magnitude exceeding a threshold computed using a two-sided z test at significance level 0.05, with standard errors estimated using 4 degrees of freedom, corresponding to the number of decadal events in the data.

Fig. 15.

(a) Sea level pressure (hPa) difference (between anomalies averaged over 60°–75°N, 30°–15°W and 40°–55°N, 25°–10°W) and (b) barotropic streamfunction (Sv) anomalies for the Labrador Sea sector (averaged over 45°–65°N, 60°–30°W region). Colored lines correspond to time series associated with running averages (annual means) of four LME simulations with just solar forcing. The black line is their mean, while the horizontal dash–dotted brown line denotes the value of zero for the given anomaly.

Fig. 15.

(a) Sea level pressure (hPa) difference (between anomalies averaged over 60°–75°N, 30°–15°W and 40°–55°N, 25°–10°W) and (b) barotropic streamfunction (Sv) anomalies for the Labrador Sea sector (averaged over 45°–65°N, 60°–30°W region). Colored lines correspond to time series associated with running averages (annual means) of four LME simulations with just solar forcing. The black line is their mean, while the horizontal dash–dotted brown line denotes the value of zero for the given anomaly.

Besides the potential of solar minima for the amplification of volcanically triggered sea ice expansion, the sensitivity of the latter to initial conditions modified by the former is also of interest. The deepest and longest solar minima follow the strongest volcanic eruptions during the last millennium, and as these two forcings interact with each other nonlinearly, it is hard to predict such a modulation based on simulations with the historical configuration of mutual timing applied in the LME. Here, we find some indications that initial conditions, and, in particular, Arctic sea ice, can play a role. Figure 16 shows the Arctic sea ice simulations for the all-forcing, volcanic, and solar runs for the time of the Samalas and Huaynaputina eruptions, two of the largest in the simulations. In the case of the 1258 Samalas eruption and the decade that followed it, it is clear that the volcanic runs started in a warmer climate with less sea ice, yet the anomalous response was as large as in the all-forcing runs. At the same time, Arctic sea ice anomalies following the 1600 Huaynaputina eruption were more pronounced in the runs with all the forcings, as compared to the runs with just volcanic forcing. Such a different response may be because Arctic sea ice anomalies were systematically positive at the moment of the Huaynaputina eruption but did not have any distinctive perturbation at the time of the Samalas eruption in the all-forcing runs, indicating that initial conditions can play a role for at least the short-term (annual to decadal) Arctic sea ice response.

Fig. 16.

Arctic sea ice area (×106 km2) for the LME (top) all-forcing, (middle) volcanic forcing, and (bottom) solar forcing runs. Year 0 is the year of the (left) Samalas and (right) Huaynaputina eruptions, both associated with the beginning of the LIA. Colored lines are the annual mean for each of the ensemble members, and the black line is their mean. The horizontal and vertical dash–dotted red lines denote the mean from the control run and the time of the given eruption, respectively. It is clear that the volcanic runs started in a warmer climate with less sea ice, but anomalously icy initial conditions in the all-forcing runs (e.g., due to solar minima) enhanced the sea ice extent in those runs.

Fig. 16.

Arctic sea ice area (×106 km2) for the LME (top) all-forcing, (middle) volcanic forcing, and (bottom) solar forcing runs. Year 0 is the year of the (left) Samalas and (right) Huaynaputina eruptions, both associated with the beginning of the LIA. Colored lines are the annual mean for each of the ensemble members, and the black line is their mean. The horizontal and vertical dash–dotted red lines denote the mean from the control run and the time of the given eruption, respectively. It is clear that the volcanic runs started in a warmer climate with less sea ice, but anomalously icy initial conditions in the all-forcing runs (e.g., due to solar minima) enhanced the sea ice extent in those runs.

Another interesting observation of nonlinear climate response to the given forcing strength and/or background conditions is seen in Fig. 16 for a time lag of 10 years with respect to the Samalas eruption. There, a pronounced sea ice expansion is observed for all ensemble members with just solar forcing. Its time series indicates that the timing of the development of this particular anomaly coincides with a deep minimum of the 11-yr solar cycle. However, other periods that overlap with solar cycle minima and shown in Fig. 16 (e.g., lag of 4 years with respect to the Samalas eruption and a lag of 2 and 13 years with respect to the Huaynaputina eruption) do not show any systematic signal in this ensemble, which may be because the corresponding anomalies of solar forcing are weaker. Such a nonlinear response of Northern Hemisphere climate to the strength of solar fluctuations has been observed recently by other studies as well (e.g., Lockwood 2012; Chiodo et al. 2016). Moreover, a deep minimum of the same strength and timing does not result in comparably abrupt expansion of Arctic sea ice in simulations with all the forcings. Such a dramatically different response can be associated with nonlinear dependence of the solar forcing tendency and background conditions of the climate, in this case, strongly perturbed by the Samalas eruption that happened at lag zero (approximately a decade earlier).

In summary, our findings support the hypothesis that volcanic, as well as solar, forcings are needed to explain the LIA-like perturbations, with their historical timing leading to initial conditions that favored a large response. Given the limited archive of the historical record of Earth’s climate, it would be a worthwhile exercise to conduct climate model simulations with the CESM-LME with solar and volcanic forcings applied separately or in combination for a range of initial conditions and/or background conditions.

8. Conclusions

This paper examines the relevance of volcanic eruptions for the initiation of Little Ice Age–like anomalies of centennial cooling and sea ice expansion for the Northern Hemisphere. For that, we examined the Last Millennium Ensemble conducted with the NCAR CESM climate model (Otto-Bliesner et al. 2016). This ensemble provides, among others, five simulations with volcanic forcing, four simulations with solar forcing, and 10 simulations with volcanic, solar, land, greenhouse gases, and orbital forcings together. We identify the Little Ice Age–like periods as those with positive anomalies of Arctic sea ice and study the decadal-to-centennial response to various internal and external perturbations of the coupled climate system simulated by the CESM. In particular, we revisit the concept of sea ice–ocean feedback and its relevance for the origin of the LIA. For that, we find that the simulated climate response requires two different mechanisms that act on different (decadal vs centennial) time scales and with strengths that differ between internally driven and externally forced climate mechanisms.

First, we examine internal drivers of Arctic sea ice variability and find that in agreement with previous studies, poleward oceanic heat transfer is significantly coupled with Arctic sea ice extent. The Atlantic meridional overturning circulation has internal fluctuations with a multidecadal time scale, and it induces pronounced Arctic sea ice anomalies in a control run and, to a smaller extent, in simulations with external forcings prescribed.

Our analysis of the LME runs with volcanic forcing indicates that sufficiently strong volcanic eruptions can significantly perturb both Arctic sea ice and the Atlantic meridional overturning circulation. Both sea ice and oceanic circulation are preferentially enhanced for a few decades following large volcanic eruptions. This suggests that internally coupled fluctuations of Arctic sea ice and AMOC, although dominant in the control run, are overwhelmed by other mechanisms relevant to sea ice extent on shorter time scales. In fact, the direct radiative forcing from volcanic eruptions and solar reduction may cool the North Atlantic and Pacific strongly enough to force the Arctic into a cold climate that persists after the forcings subside. We propose that such temporal (multidecadal) decoupling of AMOC and the Arctic can be associated with the abrupt emergence of such AMOC drivers as reduction of mid-Atlantic freshening, subpolar gyre strengthening, and intensification of deep convection within the Labrador Sea. In particular, our analysis indicates two separate mechanisms that act on different (decadal or centennial) time scales and are relevant for coupled Arctic sea ice and North Atlantic Ocean circulation evolution in different ways for externally forced and internally driven climates. Guided by our findings, we plan to apply more advanced data analysis techniques and pursue more qualitative analysis of centennial perturbation of North Atlantic/Arctic climate and its relation with shorter time scales and other regions (with particular focus on the Indo-Pacific; Slawinska and Giannakis 2017; Giannakis and Slawinska 2018).

Volcanic eruptions seem necessary for simulated initiation of the centennial cooling and sea ice expansion. In particular, after large volcanic eruptions and a few decades of AMOC strengthening, a prolonged period of AMOC weakening happens, leading to a new regime of mean overturning circulation that lasts for over a century and is significantly weaker than the mean for the control run. Positive Arctic sea ice anomalies caused by radiative cooling are reduced significantly during the time of AMOC strengthening, but in the end, moderately positive Arctic sea ice anomaly prevails along with a weaker AMOC, in agreement with the sea ice–ocean feedback mechanism.

Moreover, we find that on a centennial time scale, this model produces a Little Ice Age–like response with all forcings and with only volcanic forcing, but with a reduced magnitude for only volcanic forcing. This difference can be partially attributed to the coincidence of a small drop in solar insolation that mimics historical events of prolonged solar minima at the same time as the Little Ice Age. Interestingly, prolonged solar minima, as opposed to volcanic eruptions, seem to significantly modify atmospheric circulation in the subpolar North Atlantic and, through that, enhance volcanically triggered positive anomalies of Arctic sea ice.

The detailed attribution of volcanic and solar forcing to the emergence of Little Ice Age–like events, as well as their strength and duration, cannot be inferred here from the simulations with all the forcings; each one of the corresponding 10 simulations is only 1000 years long and forced by approximately 70 eruptions of various strength and a few solar minima of various duration, each one of them having random occurrence with respect to each other. Guided by our findings of the importance of solar and volcanic forcing, we plan to conduct a systematic range of sensitivity tests with CESM and various combinations of these two factors, identifying the dependence on their relative strength as well as mutual timing. In particular, we plan to study Little Ice Age–like events and their onset sensitivity and termination sensitivity to background conditions, as modified due to solar minima happening before the eruption. Another question to be revisited is the importance of solar minima in sustaining and determining the duration of the Little Ice Age–like anomalies from decades to centuries after an eruption.

Acknowledgments

This work is supported by NSF Grant AGS-1430051. We thank Bette Otto-Bliesner and her group for conducting and making available the Large Millennium Ensemble simulations. We thank Mitch Bushuk, Rong Zhang, Krzysztof Bartoszek, Ray Yamada, and Monika Barcikowska for useful discussions and editor Karen Shell and three reviewers for very useful comments.

REFERENCES

REFERENCES
Axford
,
Y.
, and Coauthors
,
2011
:
Do paleoclimate proxies agree? A test comparing 19 late Holocene climate and sea-ice reconstructions from Icelandic marine and lake sediments
.
J. Quat. Sci.
,
26
,
645
656
, https://doi.org/10.1002/jqs.1487.
Bauer
,
E.
,
M.
Claussen
,
V.
Brovkin
, and
A.
Huenerbein
,
2003
:
Assessing climate forcings of the Earth system for the past millennium
.
Geophys. Res. Lett.
,
30
,
1276
, https://doi.org/10.1029/2002GL016639.
Berdahl
,
M.
, and
A.
Robock
,
2013
:
Northern Hemisphere cryosphere response to volcanic eruptions in the Paleoclimate Modeling Intercomparison Project 3 last millennium simulations
.
J. Geophys. Res. Atmos.
,
118
,
12 359
12 370
, https://doi.org/10.1002/2013JD019914.
Böning
,
C. W.
,
M.
Scheinert
,
J.
Dengg
,
A.
Biastoch
, and
A.
Funk
,
2006
:
Decadal variability of subpolar gyre transport and its reverberation in the North Atlantic overturning
.
Geophys. Res. Lett.
,
33
, L21S01, https://doi.org/10.1029/2006GL026906.
Booth
,
B. B. B.
,
N. J.
Dunstone
,
P. R.
Halloran
,
T.
Andrews
, and
N.
Bellouin
,
2012
:
Aerosols implicated as a prime driver of twentieth-century North Atlantic climate variability
.
Nature
,
484
,
228
232
, https://doi.org/10.1038/nature10946.
Born
,
A.
,
T. F.
Stocker
, and
A. B.
Sandø
,
2016
:
Transport of salt and freshwater in the Atlantic Subpolar Gyre
.
Ocean Dyn.
,
66
,
1051
1064
, https://doi.org/10.1007/s10236-016-0970-y.
Branstator
,
G.
, and
H.
Teng
,
2014
:
Is AMOC more predictable than North Atlantic heat content?
J. Climate
,
27
,
3537
3550
, https://doi.org/10.1175/JCLI-D-13-00274.1.
Buckley
,
M. W.
, and
J.
Marshall
,
2016
:
Observations, inferences, and mechanisms of the Atlantic Meridional Overturning Circulation: A review
.
Rev. Geophys.
,
54
,
5
63
, https://doi.org/10.1002/2015RG000493.
Cavalieri
,
D. J.
, and
C. L.
Parkinson
,
1987
:
On the relationship between atmospheric circulation and the fluctuations in the sea ice extents of the Bering and Okhotsk Seas
.
J. Geophys. Res.
,
92
,
7141
7162
, https://doi.org/10.1029/JC092iC07p07141.
Chiodo
,
G.
,
R.
García-Herrera
,
N.
Calvo
,
J. M.
Vaquero
,
J. A.
Añel
,
D.
Barriopedro
, and
K.
Matthes
,
2016
:
The impact of a future solar minimum on climate change projections in the Northern Hemisphere
.
Environ. Res. Lett.
,
11
, 034015, https://doi.org/10.1088/1748-9326/11/3/034015.
Christiansen
,
B.
, and
F. C.
Ljungqvist
,
2011
:
Reconstruction of the extratropical NH mean temperature over the last millennium with a method that preserves low-frequency variability
.
J. Climate
,
24
,
6013
6034
, https://doi.org/10.1175/2011JCLI4145.1.
Christiansen
,
B.
, and
F. C.
Ljungqvist
,
2012
:
The extra-tropical Northern Hemisphere temperature in the last two millennia: Reconstructions of low-frequency variability
.
Climate Past
,
8
,
765
786
, https://doi.org/10.5194/cp-8-765-2012.
Christiansen
,
B.
, and
F. C.
Ljungqvist
,
2017
:
Challenges and perspectives for large-scale temperature reconstructions of the past two millennia
.
Rev. Geophys.
,
55
,
40
96
, https://doi.org/10.1002/2016RG000521.
Clement
,
A.
,
K.
Bellomo
,
L. N.
Murphy
,
M. A.
Cane
,
T.
Mauritsen
,
G.
Rädel
, and
B.
Stevens
,
2015
:
The Atlantic multidecadal oscillation without a role for ocean circulation
.
Science
,
350
,
320
324
, https://doi.org/10.1126/science.aab3980.
Crowley
,
T. J.
,
2000
:
Causes of climate change over the past 1000 years
.
Science
,
289
,
270
277
, https://doi.org/10.1126/science.289.5477.270.
Crowley
,
T. J.
, and
M. B.
Unterman
,
2013
:
Technical details concerning development of a 1200 yr proxy index for global volcanism
.
Earth Syst. Sci. Data
,
5
,
187
197
, https://doi.org/10.5194/essd-5-187-2013.
Crowley
,
T. J.
,
G.
Zielinski
,
B.
Vinther
,
R.
Udisti
,
K.
Kreutz
,
J.
Cole-Dai
, and
E.
Castellano
,
2008
:
Volcanism and the Little Ice Age
.
PAGES News
, Vol. 16, PAGES International Project Office, Bern, Switzerland, 22–23, http://www.pages.unibe.ch/download/docs/newsletter/2008-2/Special_Section/Science_Highlights/Crowley_2008-2(22-23).pdf.
Danabasoglu
,
G.
,
S.
Yeager
,
Y.
Kwon
,
J.
Tribbia
,
A.
Phillips
, and
J.
Hurrell
,
2012
:
Variability of the Atlantic meridional overturning circulation in CCSM4
.
J. Climate
,
25
,
5153
5172
, https://doi.org/10.1175/JCLI-D-11-00463.1.
Danielson
,
S.
,
E.
Curchitser
,
K.
Hedstrom
,
T.
Weingartner
, and
P.
Stabeno
,
2011
:
On ocean and sea ice modes of variability in the Bering Sea
.
J. Geophys. Res.
,
116
, C12034, https://doi.org/10.1029/2011JC007389.
Dickson
,
R. R.
,
J.
Meincke
,
S.
Malmberg
, and
A. J.
Lee
,
1988
:
The “great salinity anomaly” in the northern North Atlantic 1968–1982
.
Prog. Oceanogr.
,
20
,
103
151
, https://doi.org/10.1016/0079-6611(88)90049-3.
Dobrowolski
,
R.
,
T.
Durakiewicz
, and
A.
Pazdur
,
2002
:
Calcareous tufas in the soligenous mires of eastern Poland as an indicator of the Holocene climatic changes
.
Acta Geol. Pol.
,
52
,
63
73
, https://geojournals.pgi.gov.pl/agp/article/download/10065/8595
Dobrowolski
,
R.
,
K.
Bałaga
,
A.
Buczek
,
W. P.
Alexandrowicz
,
M.
Mazurek
,
S.
Hałas
, and
N.
Piotrowska
,
2016
:
Multi-proxy evidence of Holocene climate variability in Volhynia Upland (SE Poland) recorded in spring-fed fen deposits from the Komarów site
.
Holocene
,
26
,
1406
1425
, https://doi.org/10.1177/0959683616640038.
Dylmer
,
C. V.
,
J.
Giraudeau
,
F.
Eynaud
,
K.
Husum
, and
A.
De Vernal
,
2013
:
Northward advection of Atlantic water in the eastern Nordic Seas over the last 3000 yr
.
Climate Past
,
9
,
1505
1518
, https://doi.org/10.5194/cp-9-1505-2013.
Feulner
,
G.
,
2011
:
Are the most recent estimates for Maunder Minimum solar irradiance in agreement with temperature reconstructions?
Geophys. Res. Lett.
,
38
, L16706, https://doi.org/10.1029/2011GL048529.
Feurdean
,
A.
, and Coauthors
,
2015
:
Last millennium hydro-climate variability in central–eastern Europe (northern Carpathians, Romania)
.
Holocene
,
25
,
1179
1192
, https://doi.org/10.1177/0959683615580197.
Gałka
,
M.
,
K.
Tobolski
,
A.
Górska
,
K.
Milecka
,
B.
Fiałkiewicz-Kozieł
, and
M.
Lamentowicz
,
2014
:
Disentangling the drivers for the development of a Baltic bog during the Little Ice Age in northern Poland
.
Quat. Int.
,
328–329
,
323
337
, https://doi.org/10.1016/j.quaint.2013.02.026.
Gałka
,
M.
,
I.
Tanţău
,
V.
Ersek
, and
A.
Feurdean
,
2016
:
A 9000 year record of cyclic vegetation changes identified in a montane peatland deposit located in the eastern Carpathians (central-eastern Europe): Autogenic succession or regional climatic influences?
Palaeogeogr. Palaeoclimatol. Palaeoecol.
,
449
,
52
61
, https://doi.org/10.1016/j.palaeo.2016.02.007.
Gao
,
C.
,
A.
Robock
, and
C.
Ammann
,
2008
:
Volcanic forcing of climate over the last 1500 years: An improved ice core-based index for climate models
.
J. Geophys. Res.
,
113
, D23111, https://doi.org/10.1029/2008JD010239.
Gelderloos
,
R.
,
F.
Straneo
, and
C.
Katsman
,
2012
:
Mechanisms behind the temporary shutdown of deep convection in the Labrador Sea: Lessons from the Great Salinity Anomaly years 1968–71
.
J. Climate
,
25
,
6743
6755
, https://doi.org/10.1175/JCLI-D-11-00549.1.
Giannakis
,
D.
, and
J.
Slawinska
,
2018
:
Indo-Pacific variability on seasonal to multidecadal time scales. Part II: Multiscale atmosphere–ocean linkages
.
J. Climate
,
31
,
693
725
, https://doi.org/10.1175/JCLI-D-17-0031.1.
Gray
,
L. J.
, and Coauthors
,
2010
:
Solar influences on climate
.
Rev. Geophys.
,
48
,
RG4001
, https://doi.org/10.1029/2009RG000282.
Harada
,
N.
,
K.
Katsuki
,
M.
Nakagawa
,
A.
Matsumoto
,
O.
Seki
,
J. A.
Addison
,
B. P.
Finney
, and
M.
Sato
,
2014
:
Holocene sea surface temperature and sea ice extent in the Okhotsk and Bering Seas
.
Prog. Oceanogr.
,
126
,
242
253
, https://doi.org/10.1016/j.pocean.2014.04.017.
Hátún
,
H.
,
A. B.
Sandø
,
H.
Drange
,
B.
Hansen
, and
H.
Valdimarsson
,
2005
:
Influence of the Atlantic Subpolar Gyre on the thermohaline circulation
.
Science
,
309
,
1841
1844
, https://doi.org/10.1126/science.1114777.
Holland
,
M. M.
, and
L.
Landrum
,
2015
:
Factors affecting projected Arctic surface shortwave heating and albedo change in coupled climate models
.
Philos. Trans. Roy. Soc. London
,
373A
, 20140162, https://doi.org/10.1098/rsta.2014.0162.
Holland
,
M. M.
,
C. M.
Bitz
, and
B.
Tremblay
,
2006
:
Future abrupt reductions in the summer Arctic sea ice
.
Geophys. Res. Lett.
,
33
,
L23503
, https://doi.org/10.1029/2006GL028024.
Hunke
,
E. C.
, and
W. H.
Lipscomb
,
2008
: CICE: The Los Alamos sea ice model documentation and software, version 4.0. Los Alamos National Laboratory Tech. Rep. LA-CC-06-012, 76 pp., http://www.cesm.ucar.edu/models/ccsm4.0/cice/.
Hurrell
,
J. W.
, and Coauthors
,
2006
:
Atlantic climate variability and predictability: A CLIVAR perspective
.
J. Climate
,
19
,
5100
5121
, https://doi.org/10.1175/JCLI3902.1.
Hurrell
,
J. W.
, and Coauthors
,
2013
:
The Community Earth System Model: A framework for collaborative research
.
Bull. Amer. Meteor. Soc.
,
94
,
1339
1360
, https://doi.org/10.1175/BAMS-D-12-00121.1.
Hwang
,
Y.-T.
,
D. M. W.
Frierson
, and
J. E.
Kay
,
2011
:
Coupling between Arctic feedbacks and changes in poleward energy transport
.
Geophys. Res. Lett.
,
38
, L17704, https://doi.org/10.1029/2011GL048546.
Kaufman
,
C. A.
,
S. F.
Lamoureux
, and
D. S.
Kaufman
,
2011
:
Long-term river discharge and multidecadal climate variability inferred from varved sediments, southwest Alaska
.
Quat. Res.
,
76
,
1
9
, https://doi.org/10.1016/j.yqres.2011.04.005.
Kay
,
J.
,
M.
Holland
,
C.
Bitz
,
E.
Blanchard-Wrigglesworth
,
A.
Gettelman
,
A.
Conley
, and
D.
Bailey
,
2012
:
The influence of local feedbacks and northward heat transport on the equilibrium Arctic climate response to increased greenhouse gas forcing
.
J. Climate
,
25
,
5433
5450
, https://doi.org/10.1175/JCLI-D-11-00622.1.
Kobashi
,
T.
,
J. E.
Box
,
B. M.
Vinther
,
K.
Goto-Azuma
,
T.
Blunier
,
J. W. C.
White
,
T.
Nakaegawa
, and
C. S.
Andresen
,
2015
:
Modern solar maximum forced late twentieth century Greenland cooling
.
Geophys. Res. Lett.
,
42
,
5992
5999
, https://doi.org/10.1002/2015GL064764.
Kurahashi-Nakamura
,
T.
,
M.
Losch
, and
A.
Paul
,
2014
:
Can sparse proxy data constrain the strength of the Atlantic meridional overturning circulation?
Geosci. Model Dev.
,
7
,
419
432
, https://doi.org/10.5194/gmd-7-419-2014.
Lamentowicz
,
M.
,
M.
Gałka
,
Ł.
Lamentowicz
,
M.
Obremska
,
N.
Kühl
,
A.
Lücke
, and
V. E. J.
Jassey
,
2015
:
Reconstructing climate change and ombrotrophic bog development during the last 4000 years in northern Poland using biotic proxies, stable isotopes and trait-based approach
.
Palaeogeogr. Palaeoclimatol. Palaeoecol.
,
418
,
261
277
, https://doi.org/10.1016/j.palaeo.2014.11.015.
Larsen
,
D. J.
,
G. H.
Miller
,
Á.
Geirsdóttir
, and
S.
Ólafsdóttir
,
2012
:
Non-linear Holocene climate evolution in the North Atlantic: A high-resolution, multi-proxy record of glacier activity and environmental change from Hvítárvatn, central Iceland
.
Quat. Sci. Rev.
,
39
,
14
25
, https://doi.org/10.1016/j.quascirev.2012.02.006.
Lavigne
,
F.
, and Coauthors
,
2013
:
Source of the great A.D. 1257 mystery eruption unveiled, Samalas volcano, Rinjani Volcanic Complex, Indonesia
.
Proc. Natl. Acad. Sci. USA
,
110
,
16 742
16 747
, https://doi.org/10.1073/pnas.1307520110.
Lawrence
,
D. M.
, and Coauthors
,
2011
:
Parameterization improvements and functional and structural advances in version 4 of the Community Land Model
.
J. Adv. Model. Earth Syst.
,
3
, M03001, https://doi.org/10.1029/2011MS00045.
Lehner
,
F.
,
A.
Born
,
C. C.
Raible
, and
T. F.
Stocker
,
2013
:
Amplified inception of European Little Ice Age by sea ice–ocean–atmosphere feedbacks
.
J. Climate
,
26
,
7586
7602
, https://doi.org/10.1175/JCLI-D-12-00690.1.
Lockwood
,
M.
,
2012
:
Solar influence on global and regional climates
.
Surv. Geophys.
,
33
,
503
534
, https://doi.org/10.1007/s10712-012-9181-3.
Mann
,
M. E.
,
Z.
Zhang
,
M. K.
Hughes
,
R. S.
Bradley
,
S. K.
Miller
,
S.
Rutherford
, and
F.
Ni
,
2008
:
Proxy-based reconstructions of hemispheric and global surface temperature variations over the past two millennia
.
Proc. Natl. Acad. Sci. USA
,
105
,
13 252
13 257
, https://doi.org/10.1073/pnas.0805721105.
Marcisz
,
K.
,
W.
Tinner
,
D.
Colombaroli
,
P.
Kołaczek
,
M.
Słowiński
,
B.
Fiałkiewicz-Kozieł
,
E.
Łokas
, and
M.
Lamentowicz
,
2015
:
Long-term hydrological dynamics and fire history over the last 2000 years in CE Europe reconstructed from a high-resolution peat archive
.
Quat. Sci. Rev.
,
112
,
138
152
, https://doi.org/10.1016/j.quascirev.2015.01.019.
Marini
,
C.
, and
C.
Frankignoul
,
2014
:
An attempt to deconstruct the Atlantic multidecadal oscillation
.
Climate Dyn.
,
43
,
607
625
, https://doi.org/10.1007/s00382-013-1852-3.
Metzner
,
D.
,
S.
Kutterolf
,
M.
Toohey
,
C.
Timmreck
,
U.
Neimeier
,
A.
Freundt
, and
K.
Krüger
,
2014
:
Radiative forcing and climate impact resulting from SO2 injections based on a 200 000-year record of Plinian eruptions along the Central American Volcanic Arc
.
Int. J. Earth Sci.
,
103
,
2063
2079
, https://doi.org/10.1007/s00531-012-0814-z.
Miller
,
G. H.
, and Coauthors
,
2012
:
Abrupt onset of the Little Ice Age triggered by volcanism and sustained by sea-ice/ocean feedbacks
.
Geophys. Res. Lett.
,
39
, L02708, https://doi.org/10.1029/2011GL050168.
Moffa-Sánchez
,
P.
,
A.
Born
,
I. R.
Hall
,
D. J. R.
Thornalley
, and
S.
Barker
,
2014a
:
Solar forcing of North Atlantic surface temperature and salinity over the past millennium
.
Nat. Geosci.
,
7
,
275
278
, https://doi.org/10.1038/ngeo2094.
Moffa-Sánchez
,
P.
,
I. R.
Hall
,
S.
Barker
,
D. J. R.
Thornalley
, and
I.
Yashayaev
,
2014b
:
Surface changes in the eastern Labrador Sea around the onset of the Little Ice Age
.
Paleoceanogr.
,
29
,
160
175
, https://doi.org/10.1002/2013PA002523.
Moreno-Chamarro
,
E.
,
P.
Ortega
,
F.
González-Rouco
, and
M.
Montoya
,
2017
:
Assessing reconstruction techniques of the Atlantic Ocean circulation variability during the last millennium
.
Climate Dyn.
,
48
,
799
819
, https://doi.org/10.1007/s00382-016-3111-x.
Msadek
,
R.
,
K. W.
Dixon
,
T. L.
Delworth
, and
W.
Hurlin
,
2010
:
Assessing the predictability of the Atlantic meridional overturning circulation and associated fingerprints
.
Geophys. Res. Lett.
,
37
, L19608, https://doi.org/10.1029/2010GL044517.
Neale
,
R. B.
, and Coauthors
,
2012
:
Description of the NCAR Community Atmosphere Model (CAM 5.0). NCAR Tech. Note NCAR/TN-486+STR, 289 pp., http://www.cesm.ucar.edu/models/cesm1.0/cam/docs/description/cam5_desc.pdf
.
Okumura
,
Y.
,
C.
Deser
,
A.
Hu
,
A.
Timmermann
, and
S.
Xie
,
2009
:
North Pacific climate response to freshwater forcing in the subarctic North Atlantic: Oceanic and atmospheric pathways
.
J. Climate
,
22
,
1424
1445
, https://doi.org/10.1175/2008JCLI2511.1.
Ortega
,
P.
,
F.
Lehner
,
D.
Swingedouw
,
V.
Masson-Delmotte
,
C. C.
Raible
,
M.
Casado
, and
P.
Yiou
,
2015
:
A model-tested North Atlantic Oscillation reconstruction for the past millennium
.
Nature
,
523
,
71
74
, https://doi.org/10.1038/nature14518.
Otterå
,
O. H.
,
M.
Bentsen
,
H.
Drange
, and
L.
Suo
,
2010
:
External forcing as a metronome for Atlantic multidecadal variability
.
Nat. Geosci.
,
3
,
688
694
, https://doi.org/10.1038/ngeo955.
Otto-Bliesner
,
B. L.
, and Coauthors
,
2016
:
Climate variability and change since 850 CE: An ensemble approach with the Community Earth System Model
.
Bull. Amer. Meteor. Soc.
,
97
,
735
754
, https://doi.org/10.1175/BAMS-D-14-00233.1.
Pausata
,
F. S. R.
,
L.
Chafik
,
R.
Caballero
, and
D. S.
Battisti
,
2015
:
Impacts of high-latitude volcanic eruptions on ENSO and AMOC
.
Proc. Natl. Acad. Sci. USA
,
112
,
13 784
13 788
, https://doi.org/10.1073/pnas.1509153112.
Rahmstorf
,
S.
,
J. E.
Box
,
G.
Feulner
,
M. E.
Mann
,
A.
Robinson
,
S.
Rutherford
, and
E. J.
Schaffernicht
,
2015
:
Exceptional twentieth-century slowdown in Atlantic Ocean overturning circulation
.
Nat. Climate Change
,
5
,
475
480
, https://doi.org/10.1038/nclimate2554.
Robock
,
A.
,
1979
:
The “Little Ice Age”: Northern Hemisphere average observations and model calculations
.
Science
,
206
,
1402
1404
, https://doi.org/10.1126/science.206.4425.1402.
Robock
,
A.
,
2000
:
Volcanic eruptions and climate
.
Rev. Geophys.
,
38
,
191
219
, https://doi.org/10.1029/1998RG000054.
Rodionov
,
S. N.
,
N. A.
Bond
, and
J. E.
Overland
,
2007
:
The Aleutian Low, storm tracks, and winter climate variability in the Bering Sea
.
Deep-Sea Res. II
,
54
,
2560
2577
, https://doi.org/10.1016/j.dsr2.2007.08.002.
Sandø
,
A. B.
,
Y.
Gao
, and
H. R.
Langehaug
,
2014
:
Poleward ocean heat transports, sea ice processes, and Arctic sea ice variability in NorESM1-M simulations
.
J. Geophys. Res. Oceans
,
119
,
2095
2108
, https://doi.org/10.1002/2013JC009435.
Schleussner
,
C. F.
, and
G.
Feulner
,
2013
:
A volcanically triggered regime shift in the subpolar North Atlantic Ocean as a possible origin of the Little Ice Age
.
Climate Past
,
9
,
1321
1330
, https://doi.org/10.5194/cp-9-1321-2013.
Schmidt
,
G. A.
, and Coauthors
,
2011
:
Climate forcing reconstructions for use in PMIP simulations of the last millennium (v1.0)
.
Geosci. Model Dev.
,
4
,
33
45
, https://doi.org/10.5194/gmd-4-33-2011.
Schmidt
,
G. A.
, and Coauthors
,
2012
:
Climate forcing reconstructions for use in PMIP simulations of the last millennium (v1.1)
.
Geosci. Model Dev.
,
5
,
185
191
, https://doi.org/10.5194/gmd-5-185-2012.
Sigl
,
M.
, and Coauthors
,
2014
:
Insights from Antarctica on volcanic forcing during the Common Era
.
Nat. Climate Change
,
4
,
693
697
, https://doi.org/10.1038/nclimate2293.
Sigl
,
M.
, and Coauthors
,
2015
:
Timing and climate forcing of volcanic eruptions for the past 2500 years
.
Nature
,
523
,
543
549
, https://doi.org/10.1038/nature14565.
Slawinska
,
J.
, and
D.
Giannakis
,
2017
:
Indo-Pacific variability on seasonal to multidecadal time scales. Part I: Intrinsic SST modes in models and observations
.
J. Climate
,
30
,
5265
5294
, https://doi.org/10.1175/JCLI-D-16-0176.1.
Smith
,
R.
, and Coauthors
,
2010
: The Parallel Ocean Program (POP) reference manual: Ocean component of the Community Climate System Model (CCSM) and Community Earth System Model (CESM). Los Alamos National Laboratory Tech. Rep. LAUR-10-01853, 141 pp., http://www.cesm.ucar.edu/models/cesm1.0/pop2/doc/sci/POPRefManual.pdf.
Swingedouw
,
D.
,
P.
Ortega
,
J.
Mignot
,
E.
Guilyardi
,
V.
Masson-Delmotte
,
P. G.
Butler
,
M.
Khodri
, and
R.
Séférian
,
2015
:
Bidecadal North Atlantic ocean circulation variability controlled by timing of volcanic eruptions
.
Nat. Commun.
,
6
,
6545
, https://doi.org/10.1038/ncomms7545.
Trouet
,
V.
,
J.
Esper
,
N. E.
Graham
,
A.
Baker
,
J. D.
Scourse
, and
D. C.
Frank
,
2009
:
Persistent positive North Atlantic Oscillation mode dominated the Medieval Climate Anomaly
.
Science
,
324
,
78
80
, https://doi.org/10.1126/science.1166349.
Trouet
,
V.
,
J. D.
Scourse
, and
C. C.
Raible
,
2012
:
North Atlantic storminess and Atlantic meridional overturning circulation during the last millennium: Reconciling contradictory proxy records of NAO variability
.
Global Planet. Change
,
84–85
,
48
55
, https://doi.org/10.1016/j.gloplacha.2011.10.003.
Tulloch
,
R.
, and
J.
Marshall
,
2012
:
Exploring mechanisms of variability and predictability of Atlantic meridional overturning circulation in two coupled climate models
.
J. Climate
,
25
,
4067
4080
, https://doi.org/10.1175/JCLI-D-11-00460.1.
Vieira
,
L. E. A.
,
S. K.
Solanki
,
N. A.
Krivova
, and
I.
Usoskin
,
2011
:
Evolution of the solar irradiance during the Holocene
.
Astron. Astrophys.
,
531
,
A6
, https://doi.org/10.1051/0004-6361/201015843.
Wanamaker
,
A. D.
,
P. G.
Butler
,
J. D.
Scourse
,
J.
Heinemeier
,
J.
Eiríksson
,
K. L.
Knudsen
, and
C. A.
Richardson
,
2012
:
Surface changes in the North Atlantic meridional overturning circulation during the last millennium
.
Nat. Commun.
,
3
,
899
, https://doi.org/10.1038/ncomms1901.
Weaver
,
A. J.
, and Coauthors
,
2012
:
Stability of the Atlantic meridional overturning circulation: A model intercomparison
.
Geophys. Res. Lett.
,
39
, L20709, https://doi.org/10.1029/2012GL053763.
Wendler
,
G.
,
L.
Chen
, and
B.
Moore
,
2014
:
Recent sea ice increase and temperature decrease in the Bering Sea area, Alaska
.
Theor. Appl. Climatol.
,
117
,
393
398
, https://doi.org/10.1007/s00704-013-1014-x.
Yeager
,
S.
,
2015
:
Topographic coupling of the Atlantic overturning and gyre circulations
.
J. Phys. Oceanogr.
,
45
,
1258
1284
, https://doi.org/10.1175/JPO-D-14-0100.1.
Yeager
,
S.
, and
G.
Danabasoglu
,
2012
:
Sensitivity of Atlantic meridional overturning circulation variability to parameterized Nordic Sea overflows in CCSM4
.
J. Climate
,
25
,
2077
2103
, https://doi.org/10.1175/JCLI-D-11-00149.1.
Yeager
,
S.
, and
G.
Danabasoglu
,
2014
:
The origins of late-twentieth-century variations in the large-scale North Atlantic circulation
.
J. Climate
,
27
,
3222
3247
, https://doi.org/10.1175/JCLI-D-13-00125.1.
Zambri
,
B.
,
A.
LeGrande
,
A.
Robock
, and
J.
Slawinska
,
2017
:
Northern Hemisphere winter warming and summer monsoon reduction after volcanic eruptions over the last millennium
.
J. Geophys. Res. Atmos.
,
122
,
7971
7989
, https://doi.org/10.1002/2017JD026728.
Zanchettin
,
D.
,
C.
Timmreck
,
H.-F.
Graf
,
A.
Rubino
,
S.
Lorenz
,
K.
Lohmann
,
K.
Krüger
, and
J. H.
Jungclaus
,
2012
:
Bi-decadal variability excited in the coupled ocean–atmosphere system by strong tropical volcanic eruptions
.
Climate Dyn.
,
39
,
419
444
, https://doi.org/10.1007/s00382-011-1167-1.
Zanchettin
,
D.
,
O.
Bothe
,
H. F.
Graf
,
S. J.
Lorenz
,
J.
Luterbacher
,
C.
Timmreck
, and
J. H.
Jungclaus
,
2013
:
Background conditions influence the decadal climate response to strong volcanic eruptions
.
J. Geophys. Res. Atmos.
,
118
,
4090
4106
, https://doi.org/10.1002/jgrd.50229.
Zhang
,
L.
, and
C.
Zhao
,
2015
:
Processes and mechanisms for the model SST biases in the North Atlantic and North Pacific: A link with the Atlantic meridional overturning circulation
.
J. Adv. Model. Earth Syst.
,
7
,
739
758
, https://doi.org/10.1002/2014MS000415.
Zhang
,
L.
,
L.
Wu
, and
J.
Zhang
,
2011
:
Simulated response to recent freshwater flux change over the Gulf Stream and its extension: Coupled ocean–atmosphere adjustment and Atlantic–Pacific teleconnection
.
J. Climate
,
24
,
3971
3988
, https://doi.org/10.1175/2011JCLI4020.1.
Zhang
,
L.
,
C.
Wang
, and
S. K.
Lee
,
2014
:
Potential role of Atlantic warm pool-induced freshwater forcing in the Atlantic meridional overturning circulation: Ocean–sea ice model simulations
.
Climate Dyn.
,
43
,
553
, https://doi.org/10.1007/s00382-013-2034-z.
Zhang
,
R.
,
2015
:
Mechanisms for low-frequency variability of summer Arctic sea ice extent
.
Proc. Natl. Acad. Sci. USA
,
112
,
4570
4575
, https://doi.org/10.1073/pnas.1422296112.
Zhang
,
R.
, and
T. L.
Delworth
,
2007
:
Impact of the Atlantic multidecadal oscillation on North Pacific climate variability
.
Geophys. Res. Lett.
,
34
, L23708, https://doi.org/10.1029/2007GL031601.
Zhang
,
X.
,
M.
Prange
,
U.
Merkel
, and
M.
Schulz
,
2015
:
Spatial fingerprint and magnitude of changes in the Atlantic meridional overturning circulation during marine isotope stage 3
.
Geophys. Res. Lett.
,
42
,
1903
1911
, https://doi.org/10.1002/2014GL063003.
Zhong
,
Y.
,
G. H.
Miller
,
B. L.
Otto-Bliesner
,
M. M.
Holland
,
D. A.
Bailey
,
D. P.
Schneider
, and
A.
Geirsdóttir
,
2011
:
Centennial-scale climate change from decadally-paced explosive volcanism: A coupled sea ice-ocean mechanism
.
Climate Dyn.
,
37
,
2373
2387
, https://doi.org/10.1007/s00382-010-0967-z.

Footnotes

© 2018 American Meteorological Society. For information regarding reuse of this content and general copyright information, consult the AMS Copyright Policy (www.ametsoc.org/PUBSReuseLicenses).