Quantitative attribution of Northern Hemisphere temperatures over the past 2000 years

Abstract

Quantitative assessment of natural internal variability and externally forced responses of Northern Hemisphere (NH) temperatures is necessary for understanding and attributing climate change signals during past warm and cold periods. However, it remains a challenge to distinguish the robust internally generated variability from the observed variability. Here, large-ensemble (70 member) simulations, Energy Balance Model simulation, temperature ensemble reconstruction, and three dominant external forcings (volcanic, solar, and greenhouse gas) were combined to estimate the internal variability of NH summer (June–August) temperatures over the past 2000 years (1–2000 CE). Results indicate that the Medieval Climate Anomaly was predominantly attributed to centennial-scale internal oscillation, accounting for an estimated 104% of the warming anomaly. In contrast, the Current Warm Period is influenced mainly by external forcing, contributing up to 90% of the warming anomaly. Internal temperature variability offsets cooling by volcanic eruptions during the Late Antique Little Ice Age. These findings have important implications for the attribution of past climate variability and improvement of future climate projections.

Graphical abstract

Keywords

Common Era / Internal variability / Dark Ages Cold Period / Medieval Climate Anomaly / Current Warm Period

1 Introduction

Internal climate variability refers to inherent fluctuations caused by processes within the climate system itself such as the El Niño Southern Oscillation. Externally driven climate variability is influenced by factors outside the climate system include solar forcing, volcanic aerosols, and anthropogenic greenhouse gases. Quantification of the roles of internal variability and the externally forced responses is useful in the assessment and attribution of past and future atmospheric temperature fluctuations (Steinman et al., 2015). However, there are substantial uncertainties in their relative contributions, causing vigorous debate over attribution analysis of the observed temperature variability (Goosse et al., 2005).

The last two millennia represent ideal periods for the quantification of internal variability, as they span periods when natural forcing factors and human activities alternately dominated the climate system and attribution of its changes (Mann, 2007). Various methods are available for estimating the effects of internal variability on the climate system. For example, a proxy-reconstructed index of internal variability over the Common Era (CE) can be obtained by removing forced temperature responses attributed to volcanic and solar forcings through linear regression (Wang et al., 2017). However, nonlinear responses to these forcings cannot be fully mitigated by linear regression methods (Mann et al., 2022).

Climate model simulations provide another method for the estimation of internal variability, with two primary categories based on different simulation configurations. 1) The magnitude of internal variability can be estimated on the basis of a control experiment with no external forcing (Min et al., 2005); and 2) forced responses can be estimated by arithmetic averaging across multiple ensemble members using a single model with identical external forcings but slightly different initial conditions, followed by subtraction of the ensemble mean (the forced responses) from each ensemble member to deduce internal variability (Deser et al., 2020; Hawkins and Sutton, 2009; Suarez-Gutierrez et al., 2021; Tebaldi and Knutti, 2007). However, climate models are not expected to reproduce the phase of the observed internal variability because of the inherent randomness of internally generated variability (Marotzke and Forster, 2015).

Employment of a combination of observation data and model simulations is an optimal approach to the estimation of internal climate variability, considering that they encapsulate both internal variability and externally forced variability. When climate models yield accurate estimates of forced responses, the observed internal climate variability may be derived from the difference between proxy observations and externally forced variability, as simulated by climate models (Steinman et al., 2015). Coupled Model Intercomparison Project Phase 5 (CMIP5) climate simulations and various multi-proxy reconstructions have been applied to describe forced temperature variability over the past millennium (Mann et al., 2022; Schurer et al., 2013; Wang et al., 2022), albeit with inherent uncertainties due to the variable climate sensitivities of CMIP5 models (Meehl et al., 2005) and the imperfect proxy reconstructions (Mann et al., 2022).

We utilized an ensemble simulation of summer temperature anomalies in the Northern Hemisphere (NH) over the past 2000 years, the global LOch-Vecode-Ecbilt-CLio-agIsm Model-Large Common Era Ensemble (LOVECLIM-LCE) simulations including 70 ensemble members with full forcings (Shi et al., 2021), and produced an ensemble reconstruction of NH summer temperature anomalies based on the latest PAGES 2k temperature database, version 2 (PAGES2k Consortium, 2017). Forced responses and internal variability were explored through integration of model simulations and proxy reconstructions to further quantify their contributions to warm and cold periods in the NH summer temperature anomalies throughout the past two millennia.

2 Data and methods

2.1 Proxy-based reconstructions

The 1999 “hockey-stick” curve indicated a relatively flat temperature trend over the last millennium, followed by a sharp increase during the 20th century (Mann et al., 1999), prompting extensive academic interest. Subsequent NH temperature reconstructions have indicated that both the shape and magnitude of temperature changes over the last millennium have been influenced by the reconstruction method and the chosen proxy data (Christiansen and Ljungqvist, 2017). A fundamental challenge in temperature reconstruction revolves around the quality of proxy records. This involves addressing complexities arising from weak low-frequency signals and seasonal signal mixing.

In 2017 CE, the PAGES2k consortium released a new data set including 692 proxy records, representing the most extensive spatial coverage to date (PAGES2k Consortium, 2017). This data set was then employed in 2019 CE to reconstruct global annual mean temperature anomalies using seven different methods. Among these, the Optimal Information Extraction (OIE) method excelled in extracting low-frequency signals, allowing for clearer representation of significant multi-centennial warm/cold periods such as the Little Ice Age (LIA) and the Medieval Climate Anomaly (MCA) (Neukom et al., 2019a). This performance was achieved by considering the nonlinear process of proxy record response to climate, based on Bayesian theory (Shi et al., 2021).

Moreover, most (59%) of the 692 records were tree-ring chronologies responding mainly to growing-season temperature, which directly impacts photosynthesis and the formation of growth rings (Fritts, 1976). We employed this data set and the OIE method to reconstruct NH summer (June–August) mean temperature anomalies over the last two millennia. To account for uncertainty in the reconstruction process, we retained the calibration process iteratively 1000 times. We also employed traditional verification metrics such as the Reduction of Error (RE), Coefficient of Efficiency (CE), and Root Mean Square Error (RMSE) over an 1850–1950 CE calibration period and a 1951–2000 CE verification period to evaluate the accuracy of the reconstructions. Detailed information about the data set and the OIE method are provided by PAGES2k Consortium (2017) and Shi et al. (2021).

Furthermore, we screened previous NH temperature reconstructions to assess the accuracy and reliability of OIE reconstruction. Screening criteria included annual resolution and coverage of the last two millennia. Based on these criteria, four reconstructions were selected: Shi15 (Shi et al., 2015), Moberg05 (Moberg et al., 2005), Christiansen12 (Christiansen and Ljungqvist, 2012), and Büntgen21 (Büntgen et al., 2021).

2.2 Model-based simulations

The Last Millennium Ensemble (LME) simulation with 12 ensemble members was recently applied to estimate internal variability over the last millennium (Mann et al., 2022), based on the Community Earth System Model (CESM), version 1.1 (Otto-Bliesner et al., 2016). No other millennial model ensemble simulations have been made available to date with such a substantial number of members. However, the CESM-LME simulations, even with 12 members, may not be sufficiently comprehensive to accurately represent forced signals, with the recommended sample size for robust estimation usually involving more than 30 members (Milinski et al., 2020; Ramachandran and Tsokos, 2014).

We applied the LOVECLIM-LCE simulations with the latest LOVECLIM version 1.3, which is a fully coupled climate Earth-system model of intermediate complexity (Goosse et al., 2010). The simulation included 70 individual model runs spanning the past 2000 years, making it the largest single-model ensemble of this period to date. Consistent with Paleoclimate Modeling Intercomparison Project Phase 4 (PMIP4) recommendations, the simulation incorporated external forcings including volcanic activity, solar activity, greenhouse gases, and astronomical parameters. Each simulation began with a 1000-year model run with no external forcing to establish the initial spin-up. Further details of the LOVECLIM-LCE simulations can be found in Shi et al. (2021, 2022).

2.3 Separation of forced and internal components

The linear regression method is often used to remove forced responses from observed data, with the residual considered as internal variability (e.g., Schurer et al., 2013; Wang et al., 2017). Internal variability (T_{\mathrm{I}}) is thus calculated using Eq. (1), as follows:

(1)T_{\mathrm{R}}=T_{\mathrm{F}}+T_{\mathrm{I}},

where T_{\mathrm{R}} is the proxy-based ensemble reconstruction (OIE reconstruction) of NH temperature anomalies over the past 2000 years, containing the externally forced responses (T_{\mathrm{F}}) and the internal variability (T_{\mathrm{I}}).

One strategy for distinguishing the impact of external forcings is using a large ensemble from a single General Circulation Model (GCM) or multiple GCMs, with internal variability being reduced by averaging across different ensemble members (Deser et al., 2020; Schurer et al., 2013).

Another involves using the Energy Balance Model (EBM), which is a linear model accounting for linear responses to various forcings, based on fundamental physical principles. It simulates the balance of incoming and outgoing radiation at the Earth surface and resulting temperature changes (Mann, 2011; North et al., 1981) and has been found to be highly effective in discerning NH temperature responses to external radiative forcings (Mann et al., 2014). Factors such as volcanic eruptions, solar activity, and anthropogenic greenhouse gas have been incorporated within the EBM framework to calculate annual mean NH temperatures (Mann et al., 2014). We applied this model but with data derived specifically from the summer season rather than annual mean.

The third strategy was also applied, combining the first strategy and a multiple linear regression between main external forcings and observed temperatures, as derived from proxy reconstructions. The multiple linear regression gauges the extent to which changes in external forcings may fit temperature variability within a linear scaling framework (Knudsen et al., 2014; Wang et al., 2017). The externally forced responses (T_{\mathrm{F}}, Eq. (1)) was thus derived as follows:

(2)T_{\mathrm{F}1}={\alpha }_1{MEM}_{\mathrm{L}\mathrm{O}\mathrm{V}\mathrm{E}\mathrm{C}\mathrm{L}\mathrm{I}\mathrm{M}}^{}+{\alpha }_2,

(3)T_{\mathrm{F}2}={\beta }_{1}EBM+{\beta }_2,

(4)T_{\mathrm{F}3}={{\gamma }_{1}F}_{\mathrm{V}}+{\gamma }_2{F}_{\mathrm{S}}+{\gamma }_3{F}_{\mathrm{G}\mathrm{H}\mathrm{G}}^{}+{\gamma }_4{MEM}_{\mathrm{L}\mathrm{O}\mathrm{V}\mathrm{E}\mathrm{C}\mathrm{L}\mathrm{I}\mathrm{M}}^{}+{\gamma }_5.

The first estimate of the externally forced responses (T_{F1}, Eq. (2)) was derived from linear regression of the arithmetic Model Ensemble Mean (MEM) using the 70 members of the LOVECLIM-LCE simulations (Shi et al., 2022) ({MEM}_{\mathrm{L}\mathrm{O}\mathrm{V}\mathrm{E}\mathrm{C}\mathrm{L}\mathrm{I}\mathrm{M}}). The linear regression includes a regression coefficient ({\alpha }_1) and a constant term ({\alpha }_2, the intercept), which accounts for the offset between simulated and reconstructed temperatures.

The second estimate (T_{\mathrm{F}2}, Eq. (3)) was derived from linear regression based on a zero-dimensional EBM model including parameters {\beta }_1 and {\beta }_2 and was used instead of the direct EBM output to account for the offset between reconstructed and simulated temperatures. The parameters were consistent with those of Mann et al. (2014). External forcings in this model included volcanic (F_{\mathrm{V}}; Toohey and Sigl, 2017), solar (F_{\mathrm{S}}; Vieira et al., 2011), greenhouse gas (CO₂; F_{\mathrm{G}\mathrm{H}\mathrm{G}}; Meinshausen et al., 2017), and aerosol (F_{\mathrm{A}\mathrm{E}\mathrm{R}}; Crowley, 2000) forcings. The first three forcing data sets aligned with those utilized in the LOVECLIM-LCE simulations (Shi et al., 2022), and the aerosol forcing was the same as that used in the EBM model (Mann et al., 2014).

The third estimate (T_{\mathrm{F}3}, Eq. (4)) was generated by execution of multiple regressions onto the first three external forcings and the ensemble average from Eq. (2), representing a combination of empirical and model-based estimates. As with the previous methods, the constant term {\gamma }_5 was the intercept accounting for the offset between simulated and reconstructed temperatures. Aerosol forcing was excluded, as its effects on the EBM model from 1870 CE onward had no effect on climate prior to the industrial revolution (before 1850 CE). The factors {\gamma }_1, {\gamma }_2, {\gamma }_3, and {\gamma }_4 are regression coefficients between proxy-based reconstructions (T_{\mathrm{R}}) and the independent variables (F_{\mathrm{V}}, F_{\mathrm{S}}, F_{\mathrm{G}\mathrm{H}\mathrm{G}}, {MEM}_{\mathrm{L}\mathrm{O}\mathrm{V}\mathrm{E}\mathrm{C}\mathrm{L}\mathrm{I}\mathrm{M}}).

All three methods (T_{\mathrm{F}1}–T_{\mathrm{F}3}) were applied in estimations of forced temperature responses. The internal variability (T_{\mathrm{I}}) of NH temperature anomalies over the past 2000 years was then derived by subtracting the forced temperature signals (T_{\mathrm{F}1}–T_{\mathrm{F}3}, Eqs. (2)–4) from the proxy-based temperature reconstruction (T_{\mathrm{R}}, Eq. (1)).

The primary objective of this study was to isolate internal temperature variability from proxy-based temperature records. To achieve this, the three methods represented by Eqs. (2)–(4) were applied to estimate forced temperature signals to be subtracted from proxy-based temperatures. Among these, Eq. (4) was deemed optimal because it integrates both nonlinear influences of external forcings in the MEM simulation and linear influences of external forcings. A potential concern was the issue of multi-collinearity, which can arise if independent variables in multiple linear regressions are highly correlated, possibly distorting the interpretation of results. However, our analysis confirmed that correlations among the different independent variables of Eq. (4) never exceeded the threshold for collinearity (correlation coefficient, r < 0.7) specified by Dormann et al. (2013). Multicollinearity thus posed no substantial concern.

3 Reconstructed and simulated NH temperatures

Both reconstructed and simulated NH temperatures are required to obtain the observed internal variability. The OIE NH summer temperature reconstruction correlated well with other widely used reconstructions spanning the last two millennia (Fig.1). The robustness of this reconstruction is highlighted by its high scores for key metrics (RE = 0.79, CE = 0.57, and RMSE = 0.18).

Over the past 2000 years, the OIE NH summer temperature reconstruction revealed three distinct stages. The first millennium (1–1000 CE) was characterized by several centennial warm periods, and the second millennium was marked by several centennial cold periods (1001–1849 CE, prior to industrialization), followed by a rapid and unprecedented warming trend (1850–2000 CE, since the industrial revolution) (Fig.1(a)).

There were three historic notable warm periods, the Roman Warm Period (RWP, 1–300 CE; Ljungqvist, 2010), the MCA (850–1150 CE; Masson-Delmotte et al., 2013), and the Current Warm Period (CWP, 1919–2018 CE; Shi et al., 2021), and two cold periods, the Late Antique Little Ice Age (LALIA, 536–660 CE; Büntgen et al., 2016) and the LIA (1450–1850 CE; Masson-Delmotte et al., 2013). We assessed the uncertainties of the time spans of these warm and cold periods by combining the OIE reconstruction with the four previous reconstructions based on tree-ring (Shi15, Shi et al., 2015; and Büntgen21, Büntgen et al., 2021) and multi-proxy (Moberg05, Moberg et al., 2005; and Christiansen 12, Christiansen and Ljungqvist, 2012) records. Although three of the four reconstructions were calibrated to annual mean temperatures rather than NH summer mean temperatures, they were based mainly on tree-ring proxy records reflecting conditions during the active summer growing season. Thus, they should provide a reliable indication of summer temperature trends.

We applied a threshold of 0.5 standard deviation (SD) of temperature anomalies to define warm/cold periods for all NH reconstructions (Fig.2(a)–Fig.2(e)) and depicted uncertainties of the time spans of these periods using box plots (Fig.2(f)). Time spans of the LIA (1450–1850 CE) and CWP (1919–2018 CE) were consistent with previous studies, but the other periods had substantial uncertainties (Fig.2(f)). We defined the time spans of these warm/cold periods—i.e., the RWP (9–289 CE), MCA (851–1094 CE), and LALIA (516–609 CE)—according to the median values of the five reconstructions (shaded red and blue strips in Fig.2). No reconstruction supported the occurrence of the Dark Ages Cold Period (DACP, 400–765 CE; Helama et al., 2017).

The OIE NH summer temperature reconstruction was strongly related to other reconstructions over multidecadal-scale (Fig.1(b)), with correlation coefficients ranging from 0.45 (with Moberg05) to 0.76 (with Shi2015). Furthermore, the long-term cooling trend of the OIE NH temperature reconstruction before the industrial revolution was greater than that of all other reconstructions (Fig.1(c)). This cooling trend is possibly linked to solar insolation change associated with orbital configurations, as supported by evidence derived from regional tree-ring studies (Esper et al., 2012) and global integrated research (Neukom et al., 2019a).

We consider the ensemble NH temperature reconstructions as representative of the observed records, based on the comprehensive data sets, well-tested methods, and significant correlation with other reconstructions across the high-frequency band. Furthermore, previous studies have indicated that certain low-frequency signals may be filtered out or masked by noise within proxy data, leading to a phenomenon known as “the lost low-frequency climate signal” (Esper et al., 2002), and various statistical methods have been developed to counter this effect (e.g., Christiansen and Ljungqvist, 2017; McShane and Wyner, 2011). Note that, our findings confirm that the OIE method excels in preserving the low-frequency signal, as shown by Fig.1(c).

NH summer temperatures of the LOVECLIM-LCE simulations correlated significantly with the CESM-LME simulations (r = 0.38) during the 850–2000 CE period, albeit with a much weaker cooling trend over the second millennium (before industrialization) than that from the CESM-LME simulations (Fig.3(a)). LOVECLIM-LCE simulations based on an intermediate-complexity climate model thus have a strong performance in NH temperature simulation over the past millennium.

Simulated ensemble mean NH temperatures respond to external forcings such as volcanic, solar, and CO₂ forcings (Mann et al., 2014; Shi et al., 2022). Of these, volcanic forcing exhibited the strongest correlation with simulated temperature variability for both the LOVECLIM-LCE and CESM-LCE simulations over centennial-scale (Fig.3(b)). Solar activity had a negligible influence on ensemble mean NH temperatures of model simulations over centennial-scale (Fig.3(c)). The limited effect of solar activity has also been found over other high frequency bands (Schurer et al., 2014; Shi et al., 2022). Greenhouse gas forcing (Fig.3(d)) strongly impacted simulated NH temperatures of both models prior to the industrial revolution, although current climate models may have overestimated the climate response to greenhouse gas forcing (Hausfather et al., 2022).

4 Forced responses of NH temperatures

The multidecadal-scale has been identified as the most consistent between model-simulated and proxy-based temperature reconstructions over the past 2000 years (Neukom et al., 2019a) so we filtered out high-frequency information, leaving data from only above the multidecadal-scale (Fig.4). This helped maximize the extraction of external forcing responses from proxy-based reconstructions. Estimates of forced responses on a centennial-scale are shown in Fig.5. The correlation between proxy reconstruction and climate-model simulation is often used to assess model simulation performance (Bothe et al., 2015). Here, r² was used as a measure of the extent to which independent forcings explain the variance of dependent variables.

The ensemble mean of the LOVCLIM-LCE simulation, designated as the forced-response estimate, T_{\mathrm{F}1}, accounted for just 10% and 3% of the total variance of proxy-reconstructed NH summer mean temperature anomalies over the last two millennia at multidecadal and centennial timescales, respectively (Fig.4(b) and Fig.5(b)). The first estimate, T_{\mathrm{F}1}, primarily encompasses responses to volcanic eruptions and greenhouse gas concentrations, as represented in the LOVECLIM-LCE simulations. The lack of pronounced responses to solar forcings in this model (Fig.3(c)) likely contributes to the discrepancies observed between proxy-reconstructed and model-simulated temperatures (Ammann et al., 2007). Therefore, LOVECLIM-LCE model simulations may not accurately reproduce forced responses, particularly those related to solar activity (Shi et al., 2022; Feng et al., 2022).

The second estimate of forced responses, T_{\mathrm{F}2}, derived from the EBM model (Eq. (3)), accounted for 27% and 22% of the total variance of reconstructed NH summer mean temperatures over the last two millennia at multidecadal and centennial timescales, respectively (Fig.4(c) and Fig.5(c)) with more distinct centennial variability than the ensemble average of the LOVCLIM-LCE simulation (Fig.5(b) and Fig.5(c)). This could be due to the EBM model reflecting a relatively high contribution from solar activity, although the EBM estimate considered only linear responses of external forcings.

The third estimate (T_{\mathrm{F}3}, Eq. (4)) combined the model-based and linear response-based estimates and may be considered as best representing forced responses. Our results support this assumption, with the explained variances (36% and 58% of the total variance of reconstructed NH summer mean temperatures at multidecadal and centennial timescales, respectively) for the proxy-based NH temperatures over the last two millennia being the highest of the three sets of estimates (Fig.4(d) and Fig.5(d)). Moreover, the long-term cooling trend in T_{\mathrm{F}3} was more pronounced than the “hockey-stick” response of the ensemble mean of the LOVCLIM-LCE and EBM simulations (Fig.5(b)–Fig.5(d)). Current model simulations may thus underestimate natural externally forced NH temperature anomalies.

5 Contributions of unforced and forced variability to warm/cold periods

With the OIE proxy estimation of NH summer mean temperatures (Fig.1(a)) and the forced component of variability, T_{\mathrm{F}3}, derived from the combined method of model simulation and linear regression for forcings, we were able to obtain the unforced (internal) estimate of NH temperature variability by subtracting the forced estimate from the OIE reconstruction. Total NH temperatures (OIE reconstruction), internal variability, and forced temperature variability based on the three estimates are shown in Fig.6. Explained variances of each component during different warm and cold periods, based on T_{\mathrm{F}3}, are shown in Fig.7. Variance was calculated as the percentage ratio of unforced/forced variability to total variability (Mann et al., 2008).

The three different estimates of forced signals yielded the most notable discrepancies for the LIA, which exhibited significant cooling owing primarily to the internal variability of T_{\mathrm{F}1}. However, T_{\mathrm{F}2} and T_{\mathrm{F}3} indicated very similar contributions from forced and internal variability (Fig.6), consistent with the observation that T_{\mathrm{F}1}, based on model simulations, underestimated the impact of solar forcing that likely contributed to cooling during the LIA.

Although these three estimates may vary, they do share three common features as follows. 1) The centennial component of NH internal variability dominated the MCA (Fig.6), with internal variability explaining 104% of this warming anomaly based on the third estimate (Fig.7). 2) External forcings largely determine the CWP (Fig.6), with the T_{\mathrm{F}3} estimate accounting for 90% of this warming anomaly (Fig.7). 3) The centennial component of NH internal variability with two warm periods dominated the RWP (Fig.6), with the T_{\mathrm{F}3} estimate accounting for 60% of the warming anomaly (Fig.7). Internal variability thus appears to have played a major role during the RWP.

Regarding the LALIA, the centennial component of NH internal variability within a warm period offsets the cooling effects imposed by external forcings during this period (Fig.6), and the unforced estimate, based on T_{\mathrm{F}3}, explained 39% of the total variance (Fig.7), resulting in no apparent LALIA in proxy-based reconstructions (Fig.6(c)). It thus seems reasonable to conclude that there was no distinct LALIA/DACP, likely because of a combination of internal variability and external forcings.

6 Mechanism of the centennial-scale internal variability

Centennial-scale internal variability may be attributable to natural inherent oscillations in the climate system, such as changes in oceanic and atmospheric circulation patterns. A comparison of internal variabilities was undertaken using several indices, namely the Pacific decadal oscillation (PDO), Atlantic Multi-decadal Oscillation (AMO), North Atlantic Oscillation (NAO), and Atlantic meridional overturning circulation (AMOC). Results indicated significant correlation between centennial-scale internal variability and both the AMOC and AMO indices. Specifically, r² values of 48% were observed with the AMOC index during the industrial period (1900–2000 CE; Sun et al., 2015), and 32% with the AMO index over the last millennium (800–2000 CE; Wang et al., 2017) (Fig.8). This implies a close link between centennial-scale internal variability of NH summer temperatures and AMOC/AMO variability, driven mainly by the climate system itself.

The relationship between centennial-scale internal variability of NH summer temperatures and the AMOC/AMO circulation index may be explained by a delayed oscillator model (Sun et al., 2015), highlighting a robust dynamic interplay between the NAO and AMOC. In this relationship, a positive NAO aligns with an intensification of both the Icelandic Low and the Azores High, augmenting atmospheric forcing of the ocean. The cumulative effect of NAO atmospheric forcing results in multi-decadal AMOC variation. Subsequently, the momentum of an enhanced AMOC propels oceanic heat toward the subpolar/polar region causing subpolar warming, corresponding to a positive AMO phase and a weakening of westerly winds (Knight et al., 2005). This initiates a delayed negative feedback from AMOC to NAO, allowing the oscillation to proceed in the opposite direction. Overall, this model effectively elucidates the strong dynamical coupling between NAO and AMOC/AMO, and the impact of AMOC variations on NH temperatures (Sun et al., 2021). In addition, new evidences have established that the AMOC variability is self-sustained in the climate system over both multidecadal and multicentennial timescales (Wei and Zhang, 2022; Li and Yang, 2022).

7 Conclusions

We have quantitatively distinguished internal variability and externally forced responses of NH summer mean temperature anomalies over the last two millennia by integrating large-ensemble model simulations based on the LOVECLIM-LCE simulations and ensemble proxy reconstructions based on the OIE statistical method. The NH CWP is driven primarily by external forcings, which account for 90% of the warming anomaly. In contrast, the MCA is largely governed by internal variability, contributing up to 104% of the warming anomaly. The LALIA appears to be influenced mainly by external forcings on multidecadal-scale but is significantly offset by internal variability over centennial-scale. This likely explains the lack of a distinct DACP in proxy-based NH summer mean temperatures. Furthermore, the RWP and LIA were likely caused by a combination of internal variability and external forcings. Ultimately, this study indicates that the dominant internal climate factor driving centennial-scale fluctuations in NH summer temperatures over the last two millennia has been the AMOC, and that the primary external forcing agent is volcanic activity.

The distinct role of internal variability during the two past warm periods of the preindustrial CE lends support to the absence of naturally coherent warm periods in the NH over the past two millennia (Goosse et al., 2005; Neukom et al., 2019b). This study offers a robust tool for analysis of internal climate variability, although substantial uncertainties persist in both model simulations and proxy-based climate reconstructions. For example, the intermediate-complexity climate model has limitations compared with GCMs, and high-resolution proxy records present challenges in their low signal-to-noise ratio relative to the instrumental data.

References

[1] Ammann C M, Joos F, Schimel D S, Otto-Bliesner B L, Tomas R A (2007). Solar influence on climate during the past millennium: results from transient simulations with the NCAR Climate System Model.Proc Natl Acad Sci USA, 104(10): 3713–3718

[2] Bothe O, Evans M, Donado L F, Bustamante E G, Gergis J, Gonzalez-Rouco J F, Goosse H, Hegerl G C, Hind A, Jungclaus J H, Kaufman D S, Lehner F, Mckay N P, Moberg A, Raible C C, Schurer A P, Shi F, Smerdon J E, von Gunten L, Wagner S, Warren E, Widmann M, Yiou P, Zorita E (2015). Continental-scale temperature variability in PMIP3 simulations and PAGES 2k regional temperature reconstructions over the past millennium.Clim Past, 11(12): 1673–1699

[3] Büntgen U, Allen K, Anchukaitis K J, Arseneault D, Boucher É Bräuning A, Chatterjee S, Cherubini P, Churakova O V, Corona C, Gennaretti F, Grießinger J, Guillet S, Guiot J, Gunnarson B, Helama S, Hochreuther P, Hughes M K, Huybers P, Kirdyanov A V, Krusic P J, Ludescher J, Meier W J H, Myglan V S, Nicolussi K, Oppenheimer C, Reinig F, Salzer M W, Seftigen K, Stine A R, Stoffel M, St. George S, Tejedor E, Trevino A, Trouet V, Wang J, Wilson R, Yang B, Xu G, Esper J (2021). The influence of decision-making in tree ring-based climate reconstructions.Nat Commun, 12(1): 3411

[4] Büntgen U, Myglan V S, Ljungqvist F C, McCormick M, Di Cosmo N, Sigl M, Jungclaus J, Wagner S, Krusic P J, Esper J, Kaplan J O, de Vaan M A C, Luterbacher J, Wacker L, Tegel W, Kirdyanov A V (2016). Cooling and societal change during the Late Antique Little Ice Age from 536 to around 660 AD.Nat Geosci, 9(3): 231–236

[5] Christiansen B, Ljungqvist F C (2012). The extra-tropical Northern Hemisphere temperature in the last two millennia: reconstructions of low-frequency variability.Clim Past, 8(2): 765–786

[6] Christiansen B, Ljungqvist F C (2017). Challenges and perspectives for large-scale temperature reconstructions of the past two millennia.Rev Geophys, 55(1): 40–96

[7] Crowley T J (2000). Causes of climate change over the past 1000 years.Science, 289(5477): 270–277

[8] Deser C, Lehner F, Rodgers K B, Ault T, Delworth T L, DiNezio P N, Fiore A, Frankignoul C, Fyfe J C, Horton D E, Kay J E, Knutti R, Lovenduski N S, Marotzke J, McKinnon K A, Minobe S, Randerson J, Screen J A, Simpson I R, Ting M F (2020). Insights from Earth system model initial-condition large ensembles and future prospects.Nat Clim Chang, 10(4): 277–286

[9] Dormann C F, Elith J, Bacher S, Buchmann C M, Carl G, Carré G, Marquéz J R G, Gruber B, Lafourcade B, Leitão P J, Münkemüller T, McClean C J, Osborne P E, Reineking B, Schröder B, Skidmore A K, Zurell D, Lautenbach S (2013). Collinearity: a review of methods to deal with it and a simulation study evaluating their performance.Ecography, 36(1): 27–46

[10] Esper J, Cook E R, Schweingruber F H (2002). Low-frequency signals in long tree-ring chronologies for reconstructing past temperature variability.Science, 295(5563): 2250–2253

[11] Esper J, Frank D C, Timonen M, Zorita E, Wilson R J S, Luterbacher J, Holzkamper S, Fischer N, Wagner S, Nievergelt D, Verstege A, Büntgen U (2012). Orbital forcing of tree-ring data.Nat Clim Chang, 2(12): 862–866

[12] Feng S N, Liu X Q, Shi F, Mao X, Li Y, Wang J P (2022). Humidity changes and possible forcing mechanisms over the last millennium in arid Central Asia.Clim Past, 18(5): 975–988

[13] FrittsH (1976). Tree Rings and Climate. New Jersey: Blackburn Press

[14] Goosse H, Brovkin V, Fichefet T, Haarsma R, Huybrechts P, Jongma J, Mouchet A, Selten F M, Barriat P Y, Campin J M, Deleersnijder E, Driesschaert E, Goelzer H, Janssens I, Loutre M F, Morales Maqueda M A, Opsteegh T, Mathieu P P, Munhoven G, Pettersson E J, Renssen H, Roche D M, Schaeffer M, Tartinville B, Timmermann A, Weber S L (2010). Description of the Earth system model of intermediate complexity LOVECLIM version 1.2.Geosci Model Dev, 3(2): 603–633

[15] Goosse H, Renssen H, Timmermann A, Bradley R S (2005). Internal and forced climate variability during the last millennium: a model-data comparison using ensemble simulations.Quat Sci Rev, 24(12): 1345–1360

[16] Hausfather Z, Marvel K, Schmidt G A, Nielsen-Gammon J W, Zelinka M (2022). Climate simulations: recognize the ‘hot model’ problem.Nature, 605(7908): 26–29

[17] Hawkins E, Sutton R T (2009). The potential to narrow uncertainty in regional climate predictions.Bull Am Meteorol Soc, 90(8): 1095–1108

[18] Helama S, Jones P D, Briffa K R (2017). Dark ages cold period: a literature review and directions for future research.Holocene, 27(10): 1600–1606

[19] KnightJ R, AllanR J, FollandC K, Vellinga M, MannM E (2005). A signature of persistent natural thermohaline circulation cycles in observed climate. Geophys Res Lett, 32(20): L20708

[20] Knudsen M F, Jacobsen B H, Seidenkrantz M S, Olsen J (2014). Evidence for external forcing of the Atlantic Multidecadal Oscillation since termination of the Little Ice Age.Nat Commun, 5(1): 3323

[21] Li Y, Yang H J (2022). A theory for self-sustained multicentennial oscillation of the Atlantic Meridional Overturning Circulation.J Clim, 35(18): 5883–5896

[22] LjungqvistF C (2010). A new reconstruction of temperature variability in the extra-tropical Northern Hemisphere during the last two millennia. Geogr Ann, Ser A, 92(3): 339–351

[23] MannM E (2011). On long range dependence in global surface temperature series. Clim Change, 107(3–4): 267–276

[24] Mann M E (2007). Climate over the past two millennia.Annu Rev Earth Planet Sci, 35(1): 111–136

[25] Mann M E, Bradley R S, Hughes M K (1999). Northern hemisphere temperatures during the past millennium: inferences, uncertainties, and limitations.Geophys Res Lett, 26(6): 759–762

[26] Mann M E, Steinman B A, Brouillette D J, Fernandez A, Miller S K (2022). On the estimation of internal climate variability during the preindustrial past millennium.Geophys Res Lett, 49(2): e2021GL096596

[27] Mann M E, Steinman B A, Miller S K (2014). On forced temperature changes, internal variability and the AMO.Geophys Res Lett, 41(9): 3211–3219

[28] Mann M E, Zhang Z, Hughes M K, Bradley R S, Miller S K, Rutherford S, Ni F (2008). Proxy-based reconstructions of hemispheric and global surface temperature variations over the past two millennia.Proc Natl Acad Sci USA, 105(36): 13252–13257

[29] Marotzke J, Forster P M (2015). Forcing, feedback and internal variability in global temperature trends.Nature, 517(7536): 565–570

[30] Masson-DelmotteV, Schulz M, Abe-OuchiA, BeerJ, Ganopolski A, RoucoJ, JansenE, Lambeck K, LuterbacherJ, NaishT (2013). Information from Paleoclimate Archives. Climate Change 2013: The Physical Science Basis. In: Stocker TF, et al. eds. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change. Cambridge and New York: Cambridge University Press

[31] McShane B B, Wyner A J (2011). A statistical analysis of multiple temperature proxies: are reconstructions of surface temperatures over the last 1000 years reliable?.Ann Appl Stat, 5(1): 5–44

[32] Meehl G A, Covey C, McAvaney B, Latif M, Stouffer R J (2005). Overview of the coupled model intercomparison project.Bull Am Meteorol Soc, 86(1): 89–93

[33] Meinshausen M, Vogel E, Nauels A, Lorbacher K, Meinshausen N, Etheridge D M, Fraser P J, Montzka S A, Rayner P J, Trudinger C M, Krummel P B, Beyerle U, Canadell J G, Daniel J S, Enting I G, Law R M, Lunder C R, O’Doherty S, Prinn R G, Reimann S, Rubino M, Velders G J M, Vollmer M K, Wang R H J, Weiss R (2017). Historical greenhouse gas concentrations for climate modelling (CMIP6).Geosci Model Dev, 10(5): 2057–2116

[34] Milinski S, Maher N, Olonscheck D (2020). How large does a large ensemble need to be?.Earth Syst Dyn, 11(4): 885–901

[35] Min S K, Legutke S, Hense A, Kwon W T (2005). Internal variability in a 1000-yr control simulation with the coupled climate model ECHO-G — I. Near-surface temperature, precipitation and mean sea level pressure. Tellus A.Dynamic Meteor Oceanogr, 57(4): 605–621

[36] Moberg A, Sonechkin D M, Holmgren K, Datsenko N M, Karlén W, Lauritzen S E (2005). Highly variable Northern Hemisphere temperatures reconstructed from low- and high-resolution proxy data.Nature, 433(7026): 613–617

[37] Neukom R, Barboza L A, Erb M P, Shi F, Emile-Geay J, Evans M N, Franke J, Kaufman D S, Lücke L, Rehfeld K, Schurer A, Zhu F, Brönnimann S, Hakim G J, Henley B J, Ljungqvist F C, McKay N, Valler V, von Gunten L, PAGES 2k Consortium (2019a). Consistent multi-decadal variability in global temperature reconstructions and simulations over the Common Era.Nat Geosci, 12(8): 643–649

[38] Neukom R, Steiger N, Gómez-Navarro J J, Wang J, Werner J P (2019b). No evidence for globally coherent warm and cold periods over the preindustrial Common Era.Nature, 571(7766): 550–554

[39] North G R, Cahalan R F, Coakley J A Jr (1981). Energy balance climate models.Rev Geophys, 19(1): 91–121

[40] Otto-Bliesner B L, Brady E C, Fasullo J T, Jahn A, Landrum L, Stevenson S, Rosenbloom N, Mai A, Strand G (2016). Climate variability and change since 850 C.E.: an ensemble approach with the Community Earth System Model (CESM).Bull Am Meteorol Soc, 97(5): 735–754

[41] PAGES2k Consortium (Emile-Geay J, McKay N P, Kaufman D S, von Gunten L, Wang J H, Anchukaitis K J, Abram N J, Addison J A, Curran M A J, Evans M N, Henley B J, Hao Z X, Martrat B, McGregor H V, Neukom R, Pederson G T, Stenni B, Thirumalai K, Werner J P, Xu C X, Divine D V, Dixon B C, Gergis J, Mundo I A, Nakatsuka T, Phipps S J, Routson C C, Steig E J, Tierney J E, Tyler J J, Allen K J, Bertler N A N, Björklund J, Chase B M, Chen M T, Cook E R, de Jong R, DeLong K L, Dixon D A, Ekaykin A A, Ersek V, Filipsson H L, Francus P, Freund M B, Frezzotti M, Gaire N P, Gajewski K, Ge Q S, Goosse H, Gornostaeva A, Grosjean M, Horiuchi K, Hormes A, Husum K, Isaksson E, Kandasamy S, Kawamura K, Kilbourne K H, Koç N, Leduc G, Linderholm H W, Lorrey A M, Mikhalenko V, Mortyn P G, Motoyama H, Moy A D, Mulvaney R, Munz P M, Nash D J, Oerter H, Opel T, Orsi A J, Ovchinnikov D V, Porter T J, Roop H A, Saenger C, Sano M, Sauchyn D, Saunders K M, Seidenkrantz M S, Severi M, Shao X M, Sicre M A, Sigl M, Sinclair K, St George S, St Jacques J M, Thamban M, Kuwar Thapa U, Thomas E R, Turney C S M, Uemura R, Viau A E, Vladimirova D O, Wahl E R, White J W C, Yu Z C, Zinke J (2017). A global multiproxy database for temperature reconstructions of the Common Era. Scientific Data, 4: 170088

[42] RamachandranK M, TsokosC P (2014). Mathematical Statistics with Applications in R. London: Elsevier,

[43] Schurer A P, Hegerl G C, Mann M E, Tett S F B, Phipps S J (2013). Separating forced from chaotic climate variability over the past millennium.J Clim, 26(18): 6954–6973

[44] Schurer A P, Tett S F B, Hegerl G C (2014). Small influence of solar variability on climate over the past millennium.Nat Geosci, 7(2): 104–108

[45] Shi F, Lu H Y, Guo Z T, Yin Q Z, Wu H B, Xu C X, Zhang E L, Shi J F, Cheng J, Xiao X Y, Zhao C (2021). The position of the Current Warm Period in the context of the past 22000 years of summer climate in China.Geophys Res Lett, 48(5): e2020GL091940

[46] Shi F, Sun C, Guion A, Yin Q Z, Zhao S, Liu T, Guo Z T (2022). Roman Warm Period and Late Antique Little Ice Age in an Earth System Model Large Ensemble.J Geophys Res: Atmosph, 127(16): e2021JD035832

[47] ShiF, YangB, FengJ, Li J P, YangF M, GuoZ T (2015). Reconstruction of the Northern Hemisphere annual temperature change over the Common Era derived from tree rings. Quatern Sci, 35(5): 1051–1063 (in Chinese)

[48] Steinman B A, Mann M E, Miller S K (2015). Climate change. Atlantic and Pacific multidecadal oscillations and Northern Hemisphere temperatures.Science, 347(6225): 988–991

[49] Suarez-GutierrezL, Milinski S, MaherN (2021). Exploiting large ensembles for a better yet simpler climate model evaluation. Clim Dyn, 57(9–10): 2557–2580

[50] SunC, LiJ P, JinF F (2015). A delayed oscillator model for the quasi-periodic multidecadal variability of the NAO. Clim Dyn, 45(7–8): 2083–2099

[51] Sun C, Zhang J, Li X, Shi C M, Gong Z Q, Ding R Q, Xie F, Lou P X (2021). Atlantic Meridional Overturning Circulation reconstructions and instrumentally observed multidecadal climate variability: a comparison of indicators.Int J Climatol, 41(1): 763–778

[52] TebaldiC, KnuttiR (2007). The use of the multi-model ensemble in probabilistic climate projections. Philosoph Transact Royal Society A: Mathe, Phys Eng Sci, 365(1857): 2053–2075

[53] TooheyM, SiglM (2016). Ice core-inferred volcanic stratospheric sulfur injection from 500 BCE to 1900 CE. In: World Data Center for Climate (WDCC) at DKRZ

[54] Toohey M, Sigl M (2017). Volcanic stratospheric sulfur injections and aerosol optical depth from 500 BCE to 1900 CE.Earth Syst Sci Data, 9(2): 809–831

[55] Vieira L E A, Solanki S K, Krivova N A, Usoskin I (2011). Evolution of the solar irradiance during the Holocene.Astron Astrophys, 531(A6): A6

[56] Wang J L, Yang B, Ljungqvist F C, Luterbacher J, Osborn T J, Briffa K R, Zorita E (2017). Internal and external forcing of multidecadal Atlantic climate variability over the past 1200 years.Nat Geosci, 10(7): 512–517

[57] Wang Z Y, Wang J L, Jia J, Shi X Y, Wang S S, Pan C H (2022). Detection and attribution of summer temperature changes in China during the last millennium.Int J Climatol, 42(12): 6384–6402

[58] Wei X Y, Zhang R (2022). A simple conceptual model for the self-sustained multidecadal AMOC variability.Geophys Res Lett, 49(14): e2022GL099800

RIGHTS & PERMISSIONS

Higher Education Press