Multi-centennial Holocene climate variability in proxy records and transient model simulations

on


a b s t r a c t
Variability on centennial to multi-centennial timescales is mentioned as a feature in reconstructions of the Holocene climate.As more long transient model simulations with complex climate models become available and efforts have been made to compile large proxy databases, there is now a unique opportunity to study multi-centennial variability with greater detail and a large amount of data than earlier.This paper presents a spectral analysis of transient Holocene simulations from 9 models and 120 proxy records to find the common signals related to oscillation periods and geographic dependencies and discuss the implications for the potential driving mechanisms.Multi-centennial variability is significant in most proxy records, with the dominant oscillation periods around 120e130 years and an average of 240 years.Spectra of model-based global mean temperature (GMT) agree well with proxy evidence with significant multi-centennial variability in all simulations with the dominant oscillation periods around 120e150 years.It indicates a comparatively good agreement between model and proxy data.A lack of latitudinal dependencies in terms of oscillation period is found in both the model and proxy data.However, all model simulations have the highest spectral density distributed over the Northern hemisphere high latitudes, which could indicate a particular variability sensitivity or potential driving mechanisms in this region.Five models also have differentiated forcings simulations with various combinations of forcing agents.Significant multi-centennial variability with oscillation periods between 100 and 200 years is found in all forcing scenarios, including those with only orbital forcing.The different forcings induce some variability in the system.Yet, none appear to be the predominant driver based on

Introduction
In paleoclimate reconstructions, variability on (multi-)centennial time scales is often mentioned as a feature of the climate system, but the details behind it remain elusive and a matter of debate (Andresen et al., 2004a(Andresen et al., , 2004b;;Sepp€ a et al., 2009;Wanner et al., 2011).With the outlook on future climate change, variability on this timescale is of particular interest as the estimated response to anthropogenic forcing is often simulated until 2100 CE, for example, in the Coupled Model Intercomparison Project (CMIP6) (Eyring et al., 2016).Uncertainties in estimations of future climate change are related to both the difficulty of predicting emission scenarios and our understanding of the climate system and, thereby, how well the models can simulate it (Lee et al., 2021).If oscillations on multi-centennial timescales are a significant feature of the climate system, it is essential for future mitigation and adaptation strategies to understand the contribution of natural variability (Lohmann, 2018).
Examining climate variability on multi-centennial timescales requires long and well-resolved time series of paleoclimate data from either proxies or climate models, as even the longest instrumental records are too short.The time series need to be detailed enough to resolve the variability signals and long enough to include several oscillations, which has been a challenge for general circulation models (GCMs) as they require large amounts of computational power to simulate long timespans.Previous studies with proxy-model comparisons of paleoclimate variability using multimodel ensembles from GCMs have been restricted to simulations of the last millennium (Fern andez-Donado et al., 2013;Laepple and Huybers, 2014b;Dee et al., 2017;Parsons et al., 2017;Ljungqvist et al., 2019) and only in recent studies have longer transient simulations been used (Zhu et al., 2019;Sun et al., 2022).
Recently, more long Holocene transient simulations are becoming available through the Paleoclimate Modelling Intercomparison Project 4 (PMIP4) (Otto-Bliesner et al., 2017;Kageyama et al., 2018) and other independent works.As a result, there is now an opportunity to study multi-centennial variability (referring to both single century and multi-centennial variability in this work) in greater detail than hitherto possible.In this study, we perform a spectral analysis on temperature in transient Holocene model simulations from 9 different models to test whether they reproduce significant oscillations in the multi-centennial frequency band and if it emerges as a global or regional signal.
As the motivation for examining multi-centennial variability originates from paleoclimate reconstructions, analyzing proxy records alongside the model data is relevant.Therefore, this work will include spectral analysis of 120 proxy records allowing for a model vs. proxy comparison of the variability signals.
Proposed driving mechanisms of multi-centennial scale variability have been debated in the scientific literature.Several internal and external drivers have been suggested (Stuiver et al., 1995;Wanner et al., 2008;Marchitto et al., 2010).One external forcing often hypothesized is solar irradiance, as both proxy reconstructions (Beer, 2000;Bond et al., 2001;De Jager, 2005;Asmerom et al., 2007;Usoskin et al., 2016;Wu et al., 2018) and modelling of solar variability (Wilson, 2013) show oscillations close to or coinciding with those from paleoclimate reconstructions.Others suggested internal mechanisms such as the Atlantic Meridional Overturning Circulation (AMOC) (Stuiver et al., 1995;McDermott et al., 2001) and Arctic sea ice (Müller et al., 2012;H€ orner et al., 2017) and generally, it is challenging to distinguish between internal variability and externally forced variations (Laepple and Huybers, 2014a).Additional single/differentiated forcing simulations from five models will be analyzed to discuss the potential drivers further.
In summary, this work aims to examine the Holocene climate variability in transient simulations from 9 GCMs to 120 proxy records with an emphasis on the period of the oscillations, regional dependencies and discuss the implications for potential drivers.

Proxy records
The proxy records are collected from the "Arctic Holocene Proxy Climate Database" (AHPC) created by Sundqvist et al. (2014) and the "Temperature 12k Database" (Temp12k) created by Kaufman et al. (2020).From these databases, we select records with a resolution of at least 50 years (except for a few representing glacial parameters).After removing double entries, the proxy ensemble totals 120 records with an average temporal resolution of 28 years, distributed over the climate variables seen in Table 1 in the Appendix.
The proxy records have global coverage with a bias towards the Northern hemisphere (NH) mid-to high-latitudes, underrepresentation in Siberia, South America and Africa and no records from Australia (see Fig. 1).The details of the individual records are listed in table 4. As this work aims to study multicentennial variability in the Holocene, the records extending further back than the Holocene have the part older than 10,000 BP excluded from the analysis to avoid mixing up different climate regimes and potentially biasing the results.In the quality control done by Kaufman et al. (2020), a few of the records included in Sundqvist et al. (2014) are commented on concerning uncertainties in age estimations and local factors.These records are marked in table 4 but are still included in the analysis as all except one are still accepted by Kaufman et al. (2020).There are also comments for other records in Kaufman et al. (2020), although the comments do not necessarily provide critique as some just give additional information.For the specific details on the comments, we refer to the data published by Kaufman et al. (2020).

Transient model simulations
The model data from the nine different models are shown in Table 2 with the experimental setup for each simulation.Because the simulations are transient, they start from boundary conditions of the past climate (exact starting time varies) and then simulate the climate up through the Holocene, based on timevaried values of different forcings.The simulations are not done specifically for this collaborative work, so as Table 2 shows, some differences exist in the transient forcings included in the different models, so a complete forcing scenario or full forcing for each model is not the same across the ensemble.Time varied solar irradiance is, for example, only included in CESM1, HadCM3, MPI-ESM, and NNU 12k, while the rest use a fixed value.
Temperature is the only variable analyzed from the model data, so the temperature proxies are the most relevant group for the model-proxy comparison.Proxies mainly reflect temperature conditions near the surface, so we attempt to make the model data more comparable by using surface temperature (ts) for AWI-ESM, CESM1, LOVECLIM, NNU12k, and TraCE.We use the 2 m near surface air temperature (tas) for the other models.As we focus on climate variability, this is not expected to impact the results in the frequency domain.
All simulations cover at least 5900 years, with CESM1 being the shortest with 5900 years, and HadCM3, LOVECLIM1.3,NNU 12k and TraCE cover at least 10.000 years.NNU 12k and TraCE extend back to 12,000-and 22,000-years BP, respectively, but only the Holocene part is used for the spectral analysis.AWI-ESM consists of two separate simulations with a 500-year overlap, and we treat them as two different outputs throughout the analysis.The annual mean data from the model output is analyzed for global mean temperature (GMT) and six latitude bands HadCM3, LOVECLIM, MPI-ESM, NNU 12K, and TraCE are the five models with additional differentiated forcing simulations with a total of 14 extra simulations, and for these, we only consider GMT.As these are not explicitly created to be compared within this work, they have some differences in the simulation setup.TraCE uses a single forcing approach while the others have differentiated forcing setups in the individual simulations.We show the details about which forcings are included in the different simulations in Table 3.Since the TraCE simulations originally were more extended, the single forcing scenarios fix the other forcings at conditions from around 20.000 BP (exact time varies for different forcings).Hence, the TraCE single forcing simulations might not be as representative of the Holocene as the other three models, e.g., large NH ice sheets are still present.NNU 12k does not have a simulation with all its forcings included, so the GHG simulation is chosen to represent this model alongside the full forcing simulations from the other models.For HadCM3, the þSolVolcKK10 experiment is used as the full forcing scenario, and the þSolVolcHYDE includes the same forcings but with a different source for anthropogenic land use.MPI-ESM has a simple scenario, OrbGHG, with only orbital and greenhouse gas forcing and a full scenario (labelled all) which also includes solar irradiance, volcanic eruptions, and anthropogenic land use, although land use is only partially included in the last 1000 years (850-1849 CE) of the simulation.

Spectral analysis
As multi-centennial variability is the focus of this work, we have removed the high-frequency noise using a low-pass Butterworth filter in both the proxies and model data.As the Butterworth filter works on regularly spaced time series, the proxy data are interpolated with an Akima interpolation (Akima, 1970) to a 10 year, evenly spaced resolution as some of the original records have irregular sampling intervals.As the orininal resolution of the proxy records spans from 1 to 50 year, some records are downsampled by the resolution while others are upsampled by the interpolation.After the interpolation, the butterworth filter is used to filter out oscillations with frequencies higher than 20 years.This procedure is applied to all the records, regardless of the archive or proxy type to provide a uniform handling of the records in the pre-processing stage and effectively mainly filter the records where it is necessary (where the original resolution is better than 20 years).The procedure does open up for potential interpolation biases, but the alternative of using different approaches to different proxy types would not be free of uncertainties either and when the aim of this study is to compare a large and diverse proxy group (archive and proxy type, location, record length, etc.) the uniform approach is chosen.
The model data generally contains more high-frequency variability than the proxies due to the better temporal resolution and because some proxies naturally filter climate signals (Evans et al., 2013;Dee et al., 2015).Therefore, we set the filtering value to 60 years instead of 20 years for the model data, which is found to be more suitable even though multi-centennial variability can still be seen with the 20-year filter, as too much spectral density is distributed over the multi-decadal range, potentially clouding some of the slower multi-centennial oscillations.The upper part of the multi-decadal range is still preserved with the 60-year filter and the detail in the multi-centennial range is better.Therefore the 60-year filter is chosen for the analysis of the model data.
The spectral analysis is done with the program 'REDFIT' developed by Mudelsee (2002) in a version adapted to the R package by Bunn et al. (2021), which also includes a significance test against red noise at the 95% significance level.A large amount of data generate many spectra, and we therefore apply kernel density distributions to help illustrate the results in the different proxy and model groups.The kernel density distribution is calculated with a Fourier transform and linear approximation to estimate the density at different points (Sheather and Jones, 1991) and thereby how the significant periods are distributed over the centennial range.

Multi-centennial variability in proxy records
The spectral analysis of the proxy data reveals significant variability cycles in most records.Only 10 out of the 120 records contain no significant spectral peaks.Fig. 1 shows the location of the proxy records and a colour grouping based on the period of the significant oscillation.Only the peak with the highest spectral density is included on the map.Fig. 1 shows no apparent distinction in the oscillation period based on the geographical location (on a global scale at least).The same can be said of the 10 records with no significant oscillations.The different oscillations are not confined to specific regions, and we find the shortest (50e100 years) and longest groups (601e1000 years) in both the Northern and Southern hemispheres.The same is true for the 10 records that did not show any significant oscillations.Fig. 1 also shows the number of records in the different groups of the oscillation period.The two groups between 101 and 300 are the largest, with 31 and 27 records between 101-200 and 201e300, respectively.We provide the details of the individual records and their periods in Table 2.
As many records have more than one significant spectral peak, Fig. 2 shows the distribution of the oscillation periods found in the proxies, and includes up to three significant peaks per record.Considering all proxies in one group (Fig. 2a), the most common oscillation period at multi-centennial timescales is between 100 and 250 years, and the average is 240 years.There are also many records in the upper part of the multi-decadal range (50e100 years), but the number of significant oscillations longer than 300 years is few in comparison.
The two big groups of climate variables, temperature and SST (Fig. 2b and c), generally show the same distribution when considering them individually.However, SST is slightly flatter and less concentrated around the 100e250 years period.The other climate variables are too few to make any conclusions on their distribution, but they are not confined to any specific part of the multi-centennial range.
As the temperature records are the most extensive group, they are also made from the highest number of different archives, so they can be further sub-grouped based on the archive type, as shown in Fig. 2d.The grouping indicates that the different kinds of archives contain the multi-centennial oscillation signal, so there is no indication that the few outlying long oscillations are caused by a specific archive type, as the four longest oscillations found are from four different types (speleothems, midden, ice core and peat).
As mentioned earlier, there is a bias toward NH mid-to-highlatitudes.Despite this, the lack of dependencies on region, climate variable and archive type for the oscillation period indicates that multi-centennial variability has indeed been a large-scale feature of the Holocene climate, according to the proxy evidence.

Multi-centennial variability in transient model simulations
The evolution of the global mean temperature in transient simulations is shown in Fig. 3.There is a large spread between the models of almost 5 C at 10.000 BP and 4 C at present.To this point, it should be noted that some are surface temperature (ts) while others are near surface air temperature (tas), but neither of the two groups is consistently warmer or colder.The amplitude of variability in the different time series also varies between the models, with EC-Earth showing the highest amplitude up to ~0.5 C and down to ~0.1 C in LOVECLIM and CESM1.So, despite all being transient simulations, there are some apparent distinctions between the models, which could be related to different forcings in the simulations, and overall model performance, including model resolution.
The spectra from the full forcing simulation GMT (Fig. 4) show that all models exhibit significant spectral power in the multicentennial range, with most at frequencies corresponding to periods of 100e200 years.There are some longer oscillations up to 400 years in HadCM3 and CESM1, while AWI-ESM has the peak with the highest spectral density below 100 years.The other models have the highest spectral peak in the 100e200-year range.We also observe significant multi-decadal peaks, indicating the filter does not dampen the oscillations too much.The exact timing of the spectral peaks varies between the models alongside the spectral density.There are differences in the variability signals they produce, as indicated by their time series in Fig. 3.The overall results show a good agreement among the models regarding the presence of significant oscillations of 100e200 years in the GMT.
The spectra for the different latitude bands show significant peaks at frequencies corresponding to 100e200 years (Fig. 5), indicating the multi-centennial variability has no obvious dependencies on latitude.Comparing the latitude bands within each    and HadCM3, it still follows the same general trajectories, so there are discrepancies between how the models reproduce the Antarctic variability.This is also the case when looking at the spectral density relative to the Northern hemisphere bands, where it can be either higher or lower than the 60 o Ne30 o N band, depending on the model.However, it consistently has the highest spectral density in the SH, which could potentially be induced by the bipolar seesaw modulated by AMOC strength, if the variability is NH driven.
The models support the proxies regarding the absence of any notable latitude dependency upon having significant oscillations at multi-centennial timescales.There is an indication of the Arctic being a region of interest for feedback mechanisms or potential drivers as more spectral density is distributed over this latitude band.Further investigation of multi-centennial variability in this region is a highly relevant topic.To this point, it should also be noted that the cryospheric feedbacks in the Arctic could be a potential source of this increased variability, so these results might not be surprising.
The spectra from the differentiated forcing simulations in Fig. 6 all show significant peaks at frequencies corresponding to a 100e200-year period.This is also the case in the most basic transient setup using only orbital forcing.The other forcings all induce some variability to the system and, in some cases, alter the exact oscillation period or the trajectory of the spectra, but the multicentennial variability is always present.Anthropogenic land use seems to have a lesser impact than the other forcings, which could be due to minor change before ~4000 BC.Solar irradiance and volcanic eruptions alter the spectra, impacting the variability signals produced as they contribute to multi-decadal variability.In HadCM3, the volcanic eruptions clearly have a larger impact on multi-decadal variability than solar irradiance, so this appears to be an important forcing for the variability in this frequency range, which is also supported by Mann et al. (2021).However, the results do not support any of them being the main driver on multicentennial timescales.For example, not all the models include solar irradiance, but they still produce significant multi-centennial variability.

Proxy-model comparison
To see how the variability in the proxies compares to the model results, we plot the kernel density distribution of the proxies with those from the model latitude bands and GMT in Fig. 7 (differentiated/single forcing simulations not included).We created two kernel density distributions from the proxies, with the dark red being all proxy records and the pink being only those representing temperature, as it is the most comparable to the model results.However, Fig. 7 shows barely any difference between these two, even though the pink distribution only includes temperature and the dark red have all the climate variables from the proxy ensemble, implying that this distinction is not important.
The kernel density distributions from the proxies peak at ~120 years, so they are situated very close to those from the model data, which peak around 100 years on average.The average period for the proxies is ~240 years, and hence clearly longer than those from the models, which have an average period of around 130e160 years.In the analysis of the proxies, we also find significant oscillations with periods of up to 1000 years, so this will naturally make the average higher compared to the models, where only one oscillation longer than 400 years is found (a 410-year oscillation for 60 o Ne30 o N in HadCM3).
The kernel density distributions from the proxies are flatter than those from the models, but this is logical when considering the more diverse origin of the proxy records.Even though the climate models have their differences in both the model cores and the forcings included in these simulations, they still share the same principles that render them more homogeneous.The same cannot be said for the heterogeneous collection of proxies, where the natural processes that capture the climatic signal, the way they are reconstructed, the climatic variable they represent, and their temporal resolution all differ in the analyzed records.Some discrepancies therefore exist between the proxy and model data, as the two distributions in Fig. 7 still have their differences, so it is still possible that some modes of natural variability are not produced by the models even with solar irradiance and volcanic forcing.Combining the oscillations from the models into individual kernel density distributions (three oscillations from each latitude band and GMT) in Fig. 8 shows that HadCM3 is the model which distribution matches best with the proxies, as this is the model that most consistently produces longer periods in the different latitude bands.The distribution from MPI-ESM is also oriented more towards longer periods and is the second closest to the proxies in mean value after HadCM3.As these two models also have the complete forcing scenarios (see Table 2), with volcanic eruptions, solar irradiance, and land use, it indicates that these external forcings still are important to the variability signals, despite not finding clear evidence of any being the main driver in the differentiated forcing simulations.
Despite some discrepancies between the proxies and individual models, there is still a good agreement that there is significant multi-centennial temperature variability globally, with the majority of the oscillations having periods around 100e200 years and a ~100 longer average for the proxies due to the slightly flatter distribution.

Proxy records
Most of the proxy records and all the model data examined with spectral analysis show significant variability on multi-centennial timescales, despite having different origins and characteristics.The proxies show more variety in terms of period lengths with a more spread-out kernel density distribution in Fig. 7 and longer periods than in the model data.The proxy records represent a more diverse data group with the different proxy types and reconstructions methods used to produce the records, which may alter how variability is captured.Although this is logical reasoning, it is not an inevitable fact, as no clear distinction is found between the different groupings used in the analysis.Also the precipitation and moisture records only come from the Arctic compilation of Sundqvist et al. (2014), we cannot attest to the universal nature of multi-centennial variability in all climate variables.For example, it may be the case that monsoon rainfall shows no multi-centennial variability and that the Arctic variability only occurs due to the strong temperature dependence of snowfall in the region.Further work would be needed to explore that possibility.
The generic uncertainties related to dating, local factors and seasonality bias are potential sources of error in this analysis.Age uncertainty is not assessed in this study, but the analysis relies on the quality assessment done by Sundqvist et al. (2014) and Kaufman et al. (2020).As the analysis only depends on the relative age sequence and not the absolute values and synchronicity between the records, it is not expected to be considerable uncertainty in the results.Weaknesses in the chronology would, in this analysis, mainly obscure the results by adding noise, as we only examine the existence of low frequencies and not precisely when they occur and whether they are synchronous or not.Future work needs to consider the potential temporal changes in the variability signals in more detail, and here chronological and synchronicity uncertainties will be even more important to account for.In the original publications for the individual proxy records, multi-centennial variability is often mentioned as an apparent feature in the constructed time series, and it has only been examined with spectral analysis in a few cases.Generally, the significant peaks in the original papers are similar to the ones found in this analysis, but sometimes without complete alignment, as different methods and approaches are used for the analysis (Staubwasser et al., 2002;Sarnthein et al., 2003;Helama et al., 2010;Berner et al., 2008;Ersek et al., 2012;Yu, 2013;Kuhnert et al., 2014).Some papers on the proxy records have noted multi-centennial variability but did not quantify it with statistical analysis (Lamy et al., 2002;Andresen et al., 2004aAndresen et al., , 2004b;;Kim et al., 2004Kim et al., , 2007;;Williams et al., 2005;Anne de Vernal et al., 2005;Weldeab et al., 2007;Ojala et al., 2008;Sepp€ a et al., 2009;Belt et al., 2010;Berner et al., 2010;Martin-Chivelet et al., 2011;Larsen et al., 2012;Balascio and Bradley, 2012;Zhao et al., 2013;Salzer et al., 2014;Boldt et al., 2015;Shuman and Marsicek, 2016;Baker et al., 2017;Mezgec et al., 2017).Our analysis now confirms their initial qualitative observations.In addition, the periods in this analysis agree with a proxy study by Taricco et al. (2015), who finds common oscillations of 100e200 years in 26 proxy records across different latitudes.

Model evidence
The model results in this analysis agree with the proxy-based study by Taricco et al. (2015) regarding period length and a lack of latitudinal dependency in showing significant oscillations.However, the studies based on long transient model simulations to compare with are limited, as the simulations only recently have become available.Model performance over such long timescales is thereby also a relatively new issue.Some of the model differences in their time series and the spectral density could be related to this.It could be associated with the different experimental setups, because these simulations are not created to be a uniform ensemble.Organised coordination with specific protocols, such as standard setups with orbital, GHG and ice forcings and a realistic scenario with volcanic and solar irradiance included, would benefit studies like this one, as model ensembles are essential to mitigate model biases and uncertainties.
The TraCE simulations have been available for the longest time.The simulations have been used in a study by Zhu et al. (2019), which is the only study with spectral analysis of a transient simulation longer than the past millennium with a GCM.The study by Zhu et al. (2019) used the entire length of the TraCE simulation and the spectral peak they identified at 170e180 years is close to 160 years in our analysis (for GMT).Variability at multi-centennial timescales is also observed (without spectral analysis) in Braconnot et al. (2019) and Bader et al. (2020), and this analysis then confirms its presence.

Potential drivers
Finding the drivers of variability can be difficult, as it is hard to distinguish between forced variations projected on natural mechanisms and internal variability (Laepple and Huybers, 2014a).The different forcings in the differentiated forcing simulations do alter the distribution of spectral density slightly (some forcings more than others) in the individual models.However, no major differences are observed in having significant oscillations at multicentennial timescales.Instead, the exact period lengths change with the differentiated forcing simulations.This is also the case with the simulations that only include the slow orbital forcing, although the amount of spectral density is often situated slightly lower.The spectral analysis of the differentiated forcing simulations implies that the different forcings induce some variability to the climate system.Still, none is the primary driver of the multicentennial variability, as it is always present.
Internally driven centennial variability is in contrast to some papers from proxy reconstructions, including some used for this analysis that suggests solar irradiance be the primary driver (e.g.Sarnthein et al. (2003); McKay and Kaufman (2014); Salzer et al. (2014)).The hypothesis of solar-driven variability often stems from period alignments between climate reconstructions and solar irradiance.The Suess/De Vries solar cycle, which oscillates with 190e210 years Stefani et al. (2021), is also in good alignment with the periods found in this analysis, and it is a well-accepted feature of the sun's behaviour from both proxy studies (Beer et al., 1994;Beer, 2000;Bond et al., 2001;Ogurtsov et al., 2002;Stefani et al., 2021) and modelling of the solar system dynamics (Wilson, 2013), so it is still possible that it impacts the variability signals (Sun et al., 2022).
Solar irradiance is, however, not the only proposed driver as other internal mechanisms, such as the AMOC (Stuiver et al., 1995;McDermott et al., 2001;Rimbu et al., 2004;Wirtz et al., 2010) and Arctic sea ice (Müller et al., 2012;H€ orner et al., 2017;Lapointe et al., 2020) is suggested as potential drivers or at least related feedbacks in other proxy studies.The higher amount of spectral density in the NH high latitudes found in this analysis indicates the high relevance of driving mechanisms and/or feedbacks concentrated in this region.It is plausible that some internal mechanisms operate on multi-centennial timescales.Several climate modelling studies using long control simulations attribute the source of multicentennial variability to fluctuations of AMOC (Vellinga and Wu, 2004;Park and Latif, 2008;Friedrich et al., 2010;Delworth and Zeng, 2012;Jiang et al., 2021).Vellinga and Wu (2004) pointed out that AMOC drives the northward shift of the intertropical convergence zone (ITCZ), leading to the tropical salinity anomaly and propagating to the subpolar North Atlantic.Park and Latif (2008) and Delworth and Zeng (2012) emphasized the importance of freshwater anomalies from the South Atlantic transporting northward alongside AMOC.Some studies argue that salinity anomalies come from the Arctic Ocean dominated by freshwater exchanges between the Arctic and North Atlantic regions (Hawkins and Sutton, 2007;Yin et al., 2021).Recently, using simple conceptual models, Li and Yang (2022) identified theoretically a selfsustained multi-centennial mode in the AMOC by introducing an enhanced mixing mechanism in the subpolar ocean.
The multi-centennial variability may consist of external forcing and internal mechanisms.If the timing and period lengths do not always align, it will leave a complex signal when the contributors go in and out of phase with each other.If that is the case, the variability signal would be irregular through time.Analyzing the temporal dependency of the variability would help uncover this further, as that has not been a part of this analysis and some proxy time series are characterized by a change of frequency during the late Holocene (Allan et al., 2018;Bae and Chatzidakis, 2022).

Conclusions
This study has first attempted to analyze coherent variability signals at multi-centennial timescales in transient Holocene model simulations from GCMs and a compilation of 120 proxy records.This was done with the aim of finding the overall conclusions related to period length and potential latitudinal dependencies, as well as discussing the potential drivers through differentiated forcing simulations.
Broad scale multi-centennial variability with significant oscillation periods of 100e200 years is evident across both the model simulations and proxy data, and we found that it is a global signal.There is, however, more spectral density distributed across the NH high latitudes, which could be related to extra sensitivity or potential driving mechanisms concentrated in the region.
None of the external forcings is found to be the sole driver of the variability, which further highlights the importance of future studies of the Arctic to find potential internal drivers and/or essential feedbacks.As the external forcings still induce some variability to the system and variations in solar irradiance are a well-documented feature of the sun at multi-centennial timescales, a set of drivers may have impacted the multi-centennial variability in the Holocene.Volcanic forcing does also contribute to the variability and increases the spectral density, when included in the model, so further research is needed to examine potential drivers more closely and thoroughly.In addition, both the model and proxy data vary in their age coverage despite all being located in the Holocene and the proxies are not synchronized, so further research into any potential changes in the variability pattern through time with, e.g., wavelet analysis, is a relevant aspect of both future research into potential drivers and testing/elaborating on the results of the analysis of this paper.
(90 o Ne60 o N, 60 o Ne30 o N, 30 o Ne0 o , 30 o Se0 o , 60 o Se30 o S, and 90 o Se60 o S) to test for any potential dependencies on latitude.

Fig. 3 .
Fig. 3. Global Mean Temperature time series of the full forcing transient Holocene simulations filtered with a 60-year low-pass Butterworth filter.Age before present (BP) is relative to 1950.

Fig. 2 .
Fig. 2. Length of the significant oscillations found in the proxy data.Up to three spectral peaks are included per record.All records are included on a, b and c are separated into climate variables, and d is air temperature separated into archive type.

Fig. 4 .
Fig. 4. Spectra of model GMT from transient Holocene simulations filtered with a 60-year low-pass Butterworth filter.The red Lines show the 95% significance level from red noise.(For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Fig. 6 .
Fig.6.Spectra of GMT from differentiated forcing simulations.The dash line marks the 95% significance level relative to a red noise spectrum.The included forcings in the individual simulations can be seen in Table3.(For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Fig. 7 .
Fig. 7. Kernel density distribution of the oscillations found in the proxy and model data for the GMT and different latitude bands (up to three spectral peaks per time series).The dash line marks the mean of the group.All the groups not named proxy are from the model data.

Fig. 8 .
Fig. 8. Kernel density distribution of the oscillations found in the proxy and model data, where all latitudes bands and GMT are combined for each model (up to three spectral peaks per time series).The dash line marks the mean of the group.The upper part of the proxy distribution ( > 600 years) is not included, see Fig. 7 instead.

Table 1
Climate variables represented in the proxy records.

Table 2
Details of the models used.* Only data from 10.000 to 0 BP is used in the analysis.BP is relative to 1950.