Contrasting State-Dependent Effects of Natural Forcing on Global and Local Climate Variability

Natural forcing from solar and volcanic activity contributes significantly to climate variability. The post-eruption cooling of strong volcanic eruptions was hypothesized to have led to millennial-scale variability during Glacials. Cooling induced by volcanic eruption is potentially weaker in the warmer climate. The underlying question is whether the climatic response to natural forcing is state-dependent. Here, we quantify the response to natural forcing under Last Glacial and Pre-Industrial conditions


10.1029/2022GL098335
2 of 12 and precipitation extremes (Wang et al., 2021).These eruptions could induce a somewhat weaker response in warmer climates (Hopcroft et al., 2018), but volcanism will continue to play an important role in future variability (Bethke et al., 2017).These studies, however, do not examine the dependency of forced variability on the mean climate because they rely on future projections or the responses to large eruptions.
The paleoclimate record is crucial to assess whether a colder planet is more sensitive to natural forcing than a warmer one.Yet, temperature variability shows a mismatch between paleoclimate simulations and proxy data on the decadal-to-multicentennial scale (Ellerhoff & Rehfeld, 2021;Laepple & Huybers, 2014a).Despite major efforts, such as the Paleoclimate Modeling Intercomparison Project (PMIP; Braconnot et al., 2012;Kageyama et al., 2018), this discrepancy remains unresolved.While PMIP experiments successfully demonstrated the influ ence of natural and anthropogenic forcing on temperature variability over the last millennium (Otto-Bliesner et al., 2016), similar studies for Glacial states are missing.Transient paleoclimate simulations for the Last Glacial Maximum (LGM) have mostly been performed without high-frequency solar and volcanic forcing (Liu et al., 2009;Smith & Gregory, 2012).This lack of time-dependent forcing could potentially explain discrepancies between reconstructed and simulated variability.Additional uncertainty remains about the mechanisms of local, long-term variability (Franzke et al., 2020;Fredriksen & Rypdal, 2017;Huybers & Curry, 2006).
Separating internal and external variability has improved the understanding of climate dynamics and processes (Frankcombe et al., 2015;Haustein et al., 2019;Mann et al., 2022;Schurer et al., 2013).The approach should allow to identify drivers of local, decadal-to-multicentennial variability in cold and warm climates.This requires the comparison of unforced and forced climate simulations under Glacial and Interglacial conditions, and their validation against paleoclimate data over a wide range of timescales.Studying contributions to surface climate variability of system components that bridge internal and external factors is also necessary.Sea ice, for example, follows in extent the mean state.Natural forcing could, however, also drive multidecadal variability via modulation of the Atlantic Meridional Overturning Circulation (AMOC; Halloran et al., 2020).This highlights contributions to variability from climate components and mechanisms that bridge intrinsic and external factors.
Here, we contrast unforced and naturally forced simulations under LGM and PI conditions in an ensemble using the Hadley Centre Coupled Model Version 3.4 (HadCM3 (Gordon et al., 2000;Pope et al., 2000;Reichler & Kim, 2008;Stott et al., 2000);).We examine the mean local response of the surface climate to volcanism in the two climate states (Section 3.1).Spectral analysis (Section 3.2) further quantifies the state-and timescale-dependent effects of natural forcing on local, zonal, and global scales.It confirms a robust response to natural forcing across climate states, but a mean decline in local temperature variability with warming.To aid interpretation of the spectra, we investigate sea ice dynamics as it appears a main contributor to local, long-term variability.We validate simulated variances using proxy data (Section 3.3) to confirm that the addition of natural forcing significantly reduces the model-data mismatch on multidecadal and longer timescales.Thus, the inclusion of natural forcing provides a more accurate representation of climate variability, needed for climate simulations.
The simulations are monthly resolved and of millennial length.The boundary conditions (orography, orbital parameters, greenhouse gas concentrations) defining the mean state are constant over these runs.All runs start from the same LGM/PI spin-up simulation at consecutive years.Six runs are akin to control simulations without transient forcing.The remaining runs feature time-varying solar and volcanic forcing.We mark these with a star (*).Thus, three runs exist for each mean state and each forcing scenario (Table S2 in Supporting Information S1).Unless otherwise specified, our results represent average values of these sub-ensembles.Temperature, precipitation, sea level pressure, and wind fields are shown in Figure S2 in Supporting Information S1.The Last Glacial GMST is decidedly colder (9.5 ± 1.4) °C and the global mean precipitation rate (GMPR) is lower (935 ± 20) mm yr −1 , with a steeper equator-to-pole temperature gradient than the Pre-Industrial with (15.1 ± 1.3) °C and (1,048 ± 21) mm yr −1 , respectively.
We use transient volcanic and solar forcing (Figure S1 in Supporting Information S1), following the PMIP3 protocol for the last millennium (850-1850 CE; Schmidt et al., 2012) and updated every 10 days.We apply the same forcing in both states, as no reconstructions of solar and volcanic forcing for the LGM exist to date.This also makes comparing forced variability between the states easier.Total solar irradiance (from Steinhilber et al. (2009) and Wang et al. (2005)) is time dependent, with a superposed 11-year cycle (Schmidt et al., 2012).Volcanic forcing (Crowley & Unterman, 2013) is supplied as Aerosol Optical Depth (AOD) at four equal-area latitude bands (90-30°S, 30°S-0, 0-30°N, 30-90°N).It describes the attenuation of incoming radiation by volcanic aerosols at a wavelength of 0.55 μm and is converted into an aerosol mass loading factor (Schmidt et al., 2012).
Figure 1 shows the distribution of simulated GMST anomalies for the Last Glacial and Pre-Industrial.Forced scenarios (LGM*, PI*) exhibit larger fluctuations.In both states the GMST standard deviation increases by a factor of approximately 1.65 compared to unforced runs.There is no strong difference in the GMST distribution attributable to the mean climate.

Observations and Paleoclimate Reconstructions
We use observations and paleoclimate reconstructions to validate the variance from model simulation on interannual to multicentennial scales (2-5, 5-50, 50-200, and 200-500 years).We consider proxy records from Rehfeld et al. (2018) and the PAGES2k-Consortium (2017), and observations from the Met Office Hadley Centre's sea surface temperature data set (HadISST downloaded 11/2019;Rayner et al. (2003)).We focus on sea surface temperature observations as much of our proxy data is of marine origin.We select records that (a) are published and calibrated to temperature, (b) contain more than 50 data points, (c) cover at least three times the largest period of interest, and (d) have a mean sampling frequency of twice the highest frequency considered (Ellerhoff & Rehfeld, 2021).We exclude records with gaps larger than five times the required resolution.Ice core records are not considered on timescales below 50 years, where signal-to-noise ratios are low (Casado et al., 2020;Laepple et al., 2018).Our ensemble consists of 41 observations and 115 proxy records from six archive types (Data Set S1, Figures S9 and S10 in Supporting Information S1).

Effect Analysis
We analyze the global and local state-dependent effects of natural forcing in time and spectral domain.Following Swingedouw et al. (2017), we quantify local effects of moderate to large-magnitude volcanic eruptions using the mean standardized anomaly (MSA).The MSA represents the average value of the standardized anomalies across ensemble members.It is computed for 12-month means surrounding periods with high aerosol imprint (AOD > 0.13, corresponding to approx.−2.6 W/m 2 (Forster et al., 2021)) as follows of each individual gridbox time series X.The index i specifies the 12 months of the time series X corresponding to the set of periods T j for run j of each climate state.The normalization to the local variability σ allows detecting forced variations caused by volcanic eruptions.We test for statistical significance by bootstrapping using 400 block samples of X with a fixed length of 48 months.
We quantify the timescale-dependent variance of surface air temperature using the power spectral density (PSD, denoted spectrum).We use the multitaper method (Percival & Walden, 1993) with three windows and chi-square distributed uncertainties.The required assumption of weak stationarity (Davies & Chatfield, 1990) is reasonably fulfilled, given that we linearly detrend all time series (Fredriksen & Rypdal, 2016;Laepple & Huybers, 2014b;Nilsen et al., 2016).Spectra are smoothed logarithmically using a Gaussian kernel.Following Huybers and Curry (2006), we compute mean spectra after interpolation to the lowest resolution and binning into equally spaced log-frequency intervals.
We use variance ratios, as in Laepple and Huybers (2014b), Rehfeld et al. (2018) and Ellerhoff and Rehfeld (2021), to compare model simulations and observational data.First, observation and proxy data are interpolated onto an equidistant time axis of the same mean resolution as the raw signal.We compute the spectrum and obtain the variance by integration over the considered timescale (2-5, 5-20, 50-200, 200-500 years).Finally, we calculate the variance ratio by dividing the simulated by the reconstructed variance.Confidence intervals (CI) are obtained from a F-distribution, based on the degrees of freedom of the variance estimates.The "lgm3" and "pi2" (Table S2 in Supporting Information S1) runs are excluded for the longest timescale (200-500 years) as they are shorter than 1,000 years.Changes in variance ratios between forced and unforced runs are quantified by the areaweighted mean of the improvement factor (Appendix A).

Mean Response to Volcanic Forcing
Volcanic eruptions cause mean temperature decline at almost every location (Figure 2) as expected (Robock, 2000).The mean response, quantified by MSA, is weaker over the oceans than over land.The response is stronger between 30°N and 30°S than in high-latitude regions, following the mean AOD imprint (Figure 2c).The strongest cooling (up to three standard deviations) occurs over the Southeast Asian Archipelago (Figure 2b).These patterns are largely robust against changes in the mean climate.This also applies to precipitation, sea level pressure, and 500mbar wind speed (Figure S3 in Supporting Information S1).
The zonally averaged MSA (Figure 2 (c)) reveals small differences between the states during LGM* and PI* around the equator, 60°S, 50°N, and toward the North Pole.We identify differences in Southeast Asia, the Antarctic Ocean, over the Northern Hemisphere (NH) ice sheets, and the Barents Sea (Figures 2a and 2b).In Southeast Asia, the enhanced PI* surface climate response could be linked to the high AOD imprint from strong tropical volcanic eruptions (Fasullo et al., 2017), such as the 1,257 Samalas eruption.Changes in the land-sea mask in the region could alter the local coupling between ocean and atmosphere.In the LGM*, cooling in response to eruptions is enhanced at the Antarctic sea ice edge and in the Barents Sea.Both regions feature a higher amount of sea ice cover during the LGM.The variations in MSA extend toward the Arctic Ocean and Northern North Atlantic.Differences between the states could therefore be related to the potential for sea ice formation, likely amplifying the local response to volcanic eruptions (Timmreck, 2012).Remaining small differences are found in

Spectral Response at the Global and Local Scale
Examining power spectra for the global and local scale highlights the timescale-dependent impact of natural forcing.Global mean spectra of simulated temperatures (Figure 3a) are predominately determined by natural forcing.Including the forcing increases the power, and thus variance, on all timescales.At multidecadal scales, the forced GMST shows approximately five times more variance than unforced runs.State-dependent effects of the forced response are not discernible in these spectra.There are no pronounced spectral peaks.Enhanced centennial-scale variability in the Pre-Industrial could be attributed to a more variable AMOC (Figure S8 in Supporting Information S1).
Local mean spectra (Figure 3a) are characteristic for the mean state and less affected by natural forcing.They point to a greater temperature variance during the LGM.Differences between the states are the strongest on interannual scales, where LGM ( * ) variance is higher by a factor of approximately two compared to PI ( * ) .Zonal mean spectra (Figures 3b and 3c) reveal that the decrease in variability with warming is greatest at mid-, and especially high-latitudes, supporting a potential link to sea ice dynamics.Tropical variability widely agrees across states.Differences between forced and unforced local and zonal mean spectra are within uncertainties, but most pronounced for high-latitude, long-term variability.The PSD of global mean sea ice concentration is smaller under Pre-Industrial than Glacial conditions (Figure 3d).Above decadal scales, sea ice variability is significantly higher in forced compared to unforced scenarios.Global sea ice d

Comparison of Observed and Modeled Variability
We validate the simulated variability against observational and paleoclimate data and revisit the local, long-term variability mismatch (Ellerhoff & Rehfeld, 2021;Laepple & Huybers, 2014a;Rehfeld et al., 2018).Figure 4 shows the model-data mismatch as variance ratios.Proxy variance is increasingly larger on longer timescales compared to simulations.There is no major difference in the variance ratios between unforced and naturally forced runs on short timescales (2-5 and 5-50 years; Figures 3a and 3b).This can be explained by internal processes dominating simulated local variability at these scales.The PI simulation slightly overestimates interannual variability in the mid and high latitudes compared to sea surface temperature observations.Beyond periods of 50 years (Figures 3c-3d), simulated local variance is consistently smaller than proxy-based reconstructions.Including natural forcing in simulations decreases the mismatch for the majority of proxy sample sites.On periods of 50-200 years, the ratio bias is decreased by a factor (local mean improvement, Appendix A) of f = 1.38 (1.12, 1.71, 90% CI).The local mean improvement increases toward multicentennial scales, reducing the discrepancy.On periods of 200-500 years, the mismatch is reduced by a factor of 2.22 (1.75, 2.81).This is not sufficient to achieve consistency between modeled and proxy variance, but the mismatch is significantly smaller.

Discussion
We confirm that including natural forcing promotes temperature variability in model simulations across a range of timescales.In contrast to some experiments in the literature, we find that the modeled response of GMST does not strongly depend on the mean climate (Figures 1 and 3).Locally, weak state-dependent effects occur (Figures 2  and 3).Considering natural forcing reduces the model-data mismatch on local temperature variability, in particular on multidecadal and multicentennial scales (Figure 4).

Forced and Unforced Variability
Previous studies suggested state-dependent effects of volcanic forcing on global and hemispheric climate (Berdahl & Robock, 2013;Muthers et al., 2015;Swingedouw et al., 2017;Zanchettin et al., 2016).These results were obtained using ensembles of large volcanic eruptions.State dependency in these has been primarily linked to nonlinear processes and initial conditions (Zanchettin et al., 2013).We argue that the response to individual volcanic eruptions may well depend on the climate state.However, globally averaged effects from changes in response mechanisms are small considering realistic forcing scenarios, in line with a linear relation between GMST and external forcing (Fredriksen & Rypdal, 2017;Geoffroy et al., 2013;MacMynowski et al., 2011).In our ensemble, the GMST response to an eruption of the size of Samalas 1,257 shows no difference between LGM* and PI* (Figure S7a in Supporting Information S1).Global precipitation and sea ice concentration is only slightly enhanced in the LGM* (Figures S7b and S7c in Supporting Information S1).Interannual GMST variability and equilibrium climate sensitivity (ECS) can also be linked in a linear feedback framework (Cox, 2001).Assuming the above framework, and that HadCM3 simulates all relevant feedbacks, the similar forced response could indicate that ECS in LGM and PI are not strongly different.

Internal Variability Across States
The question of state-dependent variability has long motivated studies of past (Ditlevsen et al., 1996;Rehfeld et al., 2018;Shao & Ditlevsen, 2016) and future (Huntingford et al., 2013;Olonscheck et al., 2021;Rehfeld et al., 2020) climate.Our results reveal a decrease in mean local variability with warming (Figure 3).Decreasing sea ice dynamics and a smaller meridional temperature gradient are suggested as major causes.In line with other studies (Bathiany et al., 2018;Berdahl & Robock, 2013;Bethke et al., 2017;Brown et al., 2017;Olonscheck et al., 2021;Rehfeld et al., 2018), we find a clear zonal pattern, with greater reduction of variability in the mid and high latitudes (Figures 3b and 3c).This is corroborated by the small discrepancy between short-term variability from observations and simulations in the mid and high latitudes (Figure 4a).The sea surface temperature observations contain the recent global warming trend and sea ice retreat, our PI(*) simulations do not.This could contribute to the decrease in local, high-latitude variability.The mean climate also changes AMOC variability in HadCM3 (Figure S8 in Supporting Information S1).It is smaller in the LGM on multidecadal and multicentennial scales (Jackson & Vellinga, 2013).Under LGM conditions, the AMOC strength and correlation length is also increased by natural forcing (Figure S8 in Supporting Information S1).Potential mechanisms of the intensification are debated (Iwi et al., 2012;Mignot et al., 2011), and they could contribute to state-dependent enhancement of long-term regional variability through natural forcing.

Mechanisms Leading to Long-Term Variability
Across our experiments, sea ice variability and regions with varying sea ice extent, primarily the Southern Oceans and Barents Sea, are most affected by natural forcing.This is further supported by mean standardized anomalies of precipitation, sea level pressure, and wind speed over the North American ice sheet, the North Atlantic Ocean, Antarctica, and the Southern Oceans (Figure S3 in Supporting Information S1).Comparing simulations with the two-dimensional energy balance model (TransEBM; Ziegler and Rehfeld (2021); Figure S5 in Supporting Information S1) adds support to the role of sea ice in forced temperature variability.TransEBM is a fairly linear model with no atmospheric and oceanic dynamics.We use it to differentiate the contribution from deterministic forcing and sea ice to the variance.In TransEBM experiments we prescribe sea ice changes from HadCM3.Forming the ratio of the local mean TransEBM and HadCM3 variability (Figure S6 in Supporting Information S1) supports the strong sea ice contribution to interannual variability (sea also Figure 3a).The contribution remains significant on decadal and longer timescales, promoting sea ice variations as a key mechanism of local, long-term variability.
Our results provide crucial insights into the discrepancy between modeled and reconstructed local, long-term variability (Ellerhoff & Rehfeld, 2021;Laepple & Huybers, 2014a).Internal variability dominates the local temperature variance on annual to decadal scales (Goosse et al., 2005), but contributions from natural forcing are detectable beyond decadal timescales (Figure 3), increasing variance on longer timescales.This is supported by increased scaling coefficients (Figure S4 in Supporting Information S1), and, hence, increased persistence of forced temperatures on periods of 50-500 years, similar to (Vyushin et al., 2004).Including natural forcing in simulations improved model-data agreement of local variability on multidecadal and multicentennial scales (Figure 4).This is perhaps surprising given that the forcing has no centennial scale variability (Ellerhoff & Rehfeld, 2021).There is no change in agreement from interannual timescales, implying that the gain from forcing on local temperatures is small on these short timescales.Hence not only the integrated response to strong (Timmreck, 2012) but also to weak natural forcing contributes to long-term variability.Time-varying forcing appears thus beneficial for reliable simulations of global mean (Figure 3) and local, long-term variability.
Consistent with previous arguments (Bethke et al., 2017), our results challenge the common usage of external forcing that shows no time-varying changes (O'Neill et al., 2016).

Limitations and Potential
We used the same volcanic forcing reconstruction to drive simulations in LGM and PI, in an idealized setup.The volcanic record shows that large eruptions occurred throughout the Glacial (Brown et al., 2014).We do not account for the possibly lower rate than during the last millennium, which could also be due to undersampling.Furthermore, we may miss feedbacks in HadCM3 that are relevant for local climate variability.This could explain the underestimation of local variability compared to proxy data (Figure 4a).Sea ice dynamics, stratospheric and cloud-related feedbacks are key nonlinear mechanisms that can alter the response to volcanic forcing in a warmer climate (Aubry et al., 2021;Fasullo et al., 2017;Hopcroft et al., 2018).Projections of tropical eruptions with a newer model show enhanced (dampened) radiative forcing from strong (moderate) eruptions (Aubry et al., 2021).Cloud-related feedback, likely underestimated in HadCM3, is generally weaker than feedback from sea ice, but may be enhanced with warming (Hopcroft et al., 2018).
Future work could examine the response in simulations with models that show a higher ECS (Gettelman et al., 2019;Tatebe et al., 2019;Voldoire et al., 2019;Wu et al., 2019) and better sea-ice (Guarino et al., 2020).Insufficient sea ice and vegetation cover changes may significantly alter the response in extreme warming scenarios.Future studies could test long-term impacts of volcanism and local state dependency with more advanced climate models, including better representation of radiative-chemical feedback, aerosol indirect, stratospheric and sea ice processes.Moreover, future research could apply probabilistic eruption projections (Bethke et al., 2017) in larger ensembles (Zanchettin et al., 2016) to study forcing scenarios with localized eruptions at high and low latitude.This will aid understanding of long-term Earth-system sensitivity and the state-dependent response of multidecadal modes to natural forcing (Swingedouw et al., 2017).

Conclusion
Presenting the first millennial-length, naturally forced simulation for the LGM, we investigated state-dependent effects of volcanic and solar forcing on global and local climate variability.The modeled global temperature response shows no dependency on the mean climate.Weak local differences resulted primarily from sea ice dynamics, providing a key mechanism of long-term variability.Including natural forcing in climate model simulations improves the agreement between modeled and observed variability and, thus, calls into question the lack of time-dependent volcanic forcing in projections and model-data comparison.The robust temperature response suggests that findings on the ability of models to simulate past variability can help constrain forced variability across spatial and temporal scales.

Appendix A: Variance Ratio Improvement
We quantify the change in variance ratio r from unforced and naturally forced simulations to proxy records using the logarithmic measure l( x .We convert the logarithmic distance to the variance ratio improvement f = 10 Δl and estimate confidence intervals using the area-weighted mean of the error propagation We ensure a conservative coverage of the CIs by using the upper limit on   ( * )  from F-distributed uncertainties of the variance ratio estimate.This research used the Archer UK National Supercomputing Services.It benefited from discussions within the CVAS working group, a working group of the Past Global Changes (PAGES) project.We thank M. Casado, T. Laepple and A. Schurer for discussion, and C. Wirths for setting up volcanic forcing over latitude intervals in TransEBM.Funding is acknowledged by the PalMod project (www.palmod.de,subproject no.01LP1926C), the Deutsche Forschungsgemeinschaft (German Research Foundation) -project no.395588486 and no.316076679, and by the Heinrich-Böll-Stiftung.Open access funding enabled and organized by Projekt DEAL.

Figure 1 .
Figure1.Probability density (unitless) of simulated yearly global mean surface temperature anomalies from all Pre-Industrial and Last Glacial Maximum runs.Forced scenarios are marked with a (*).The ratio of the distributions' standard deviations is given by r σ .

Figure 2 .
Figure 2. (a) Mean standardized anomalies (MSA) of surface air temperature in the Last Glacial Maximum* (b) and the Pre-Industrial* state after volcanic eruptions.Dots indicate insignificant anomalies within the 99% quantile range of local variability.Gray shaded crosses show land ice.Hatched areas indicate areas with >50% yearly sea ice coverage.(c) Zonally averaged MSA and Aerosol Optical Depth (AOD) (black dashed).MSA and AOD are unitless quantities.

Figure 3 .
Figure 3. (a) Local and global power spectral density for simulated temperature using Hadley Centre Coupled Model Version 3.4.Global spectra are computed from global mean surface temperature.Local refers to the area-weighted average of all local spectra.(b and c) Area-weighted average of local spectra by climate zone, given by the tropics (−23.5 to 23.5°N), mid (23.5-66.5°),and high latitudes (>66.5°) for Last Glacial Maximum and Pre-Industrial.(d) Global spectra (units in % 2 yr) of global mean sea ice concentration, defined as percentage of the globe covered in sea ice.Lines show logarithmically smoothed (0.08 dB) mean spectra with shaded 95% CIs.

Figure 4 .
Figure 4. Ratio r(sim./obs.) of simulated to observed variance over latitude for unforced (black) and naturally forced (green) Hadley Centre Coupled Model Version 3.4 simulations for the Pre-Industrial.Model data is bilinearly interpolated to the location of the observation.We show the ratio of simulated PI temperature to observations for periods of 2-5 years (a), and to proxies spanning the last 8,000 years on interannual to decadal (b), multidecadal (c), and multicentennial (d) timescales.Symbols indicate the variance ratio and vertical lines their 90% CI.The local mean improvement f of the variance ratio is given in the lower left of each panel, with CIs in parentheses (Appendix A).
at site i, with (*) denoting the climate state.The distance = * − ( ) gives the change of the variance ratio bias between the forced and unforced simulation.For a set of N sites, we quantify the mean change  Δ = 1  ∑    with local area weights w i derived from the HadCM3 grid  ( ∑    = 1 )