Orbital Forcing Strongly Influences Seasonal Temperature Trends During the Last Millennium

Insolation changes caused by the axial precession induce millennial trends in last millennium temperature, varying with season and latitude. A characteristic seasonal trend pattern can be detected in both insolation and modeled surface temperature response. In the extratropical Northern Hemisphere, the maximum insolation trend occurs around April/May, while the minimum trend occurs between July and September. The temperature trend lags behind insolation trend by around a month. Hence orbital forcing potentially affects long‐term trends in proxy data, which are often sensitive to a distinct seasonal window. We find that tree‐ring reconstructions based on early growing season dominated records show different millennial trends from those for late summer dominated proxies. The differential response is similar to that seen in pseudo proxy reconstructions when considering proxy seasonality. This suggests that orbital forcing has influenced long‐term trends in climate proxies. It is therefore vital to use seasonally homogeneous data for reconstructing multicentennial variability.

example, Kaufman et al. (2009) found a pervasive cooling in high-latitude proxy reconstructions of summer temperature and inferred that it was caused by orbital forcing. McKay & Kaufman (2014) obtained an even higher negative trend in reconstructed annual temperature of the same area, although they questioned the suitability of their reconstruction to determine long-term variability. Esper et al. (2012) found consistent trends around −0.2 to −0.3 K/ky in Scandinavian latewood density (late summer) proxies, in agreement with model data, and attributed these to orbital forcing. However, they could not find similar trends in data based on ring width, and suggested that ring width data could be limited in preserving long-term trends.  found that simulations which are only subjected to orbital forcing could not replicate the long-term trends previously observed by Esper et al. (2012) and Kaufman et al. (2009) for global annual mean temperature. They suggest that seasonal or geographical biases in proxy reconstructions could potentially be responsible for the larger trends, as proxy reconstructions tend to be more representative of higher latitudes and some records are biased toward summer. They obtained stronger long-term trends in their simulations, when taking these biases into account. Klippel et al. (2020) investigated the magnitude of the trend displayed by proxy records and their inherent large-scale properties including latitude (mid vs. high) or seasonal sensitivity (summer vs. annual temperature), but found no significant relationship. However, they found that millennial trends in ice cores, marine, and lake sediments were stronger than in tree-rings. This discrepancy between proxy data and model simulations, but also between different types of proxy data, has so far remained unresolved.
While authors have compared mainly tree-ring width against latewood density data, more specific characteristics inherent to the data sources have not been distinguished. In this study, we attempt to group proxy data according to their seasonal sensitivity to find an explanation for the aforementioned discrepancies. The exact seasonality of large-scale proxy reconstructions is often difficult to establish, as data sources with different monthly sensitivity are mixed and the exact seasonal sensitivity of an individual record is uncertain. Most reconstructions of past surface temperatures have been calibrated to either summer (e.g. Anchukaitis et al., 2017;Briffa, 2000;Kaufman et al., 2009;Schneider et al., 2015;Wilson et al., 2016) or annual (e.g. D'Arrigo et al., 2006Frank et al., 2007;McKay & Kaufman, 2014;Neukom et al., 2019) instrumental temperature data. Esper et al. (2005) found that the target seasonality plays a role for the amplitude of the reconstruction and that amplitudes are higher when scaling to annual temperature compared to summer temperatures. The decreased amplitude for summer could be explained by early instrumental biases (Böhm et al., 2010). Several studies have suggested that reconstructions of seasonal versus annual temperatures are otherwise statistically indistinguishable on multidecadal and longer timescales (Cook et al., 2004;Esper, 2002;D'Arrigo et al., 2006). However, given that tree growth is related to growing season conditions, it was recommended to calibrate tree-ring reconstructions to summer temperature (Wilson et al., 2007(Wilson et al., , 2016. In this study, we investigate the effect of the seasonality of insolation trends due to orbital forcing on simulated and reconstructed temperature over the preindustrial last millennium (AD 850-1850). We find that orbital forcing induces a characteristic seasonal trend pattern in climate models, which lags behind the insolation by around a month. In contrast, negative trends arising from other forcings are fairly constant throughout the whole year. A similar seasonal trend pattern is found in proxy reconstructions when accounting for their different seasonal sensitivity, which could be explained by orbital forcing. Our results could not only resolve certain model-proxy discrepancies, but also provide an explanation for discrepancies between different proxy types, such as ring width and density data.

Orbital Data
We use the numerical solutions for the eccentricity e, the climatic precession ψ, and the obliquity ϵ and insolation quantities provided by Laskar et al. (2004). The insolation is given as the monthly mean on a longitudinal grid of 5° resolution. To first order, e and ϵ vary at periods of around 400,000 and 441,000 years (Berger & Loutre, 1991) and are in decline over the entire Holocene ( Figure S1). Given their large periodicity, e can be assumed to be approximately constant for our analysis. ϵ is responsible for a small decrease in global annual insolation (around 0.01 W m −2 ), and introduces a zonal dependence to the millennial insolation trend. In contrast, ψ induces major changes in seasonal insolation on shorter timescales, with a periodicity of around 24,000 years to first order (Berger & Loutre, 1991). The global seasonal insolation exhibits the same periodicity as the precession, following the precession with different phase shifts.
We compare these results with single forcing simulations from the CESM-LME. The simulations use a version of CESM-CAM5_CN (1.9 × 2.5_gx1v6), with a resolution of ∼2° in atmosphere and land components and ∼1° resolution in ocean and sea ice components. Besides the 13 simulations including all forcings, we use two unforced control simulations to calculate the model drift and several experiments including single forcings. Given that the drift is very small, we do not subtract it from the simulation before calculating longterm trends. We compare long-term trends from orbital forcing only with the influence of other forcings, calculated as the sum of volcanic, solar, greenhouse gas, and land use/land cover changes (three simulations for each forcing).

Instrumental Data
For the calibration of the proxy reconstructions, we use CRUTEM4 , a gridded data set of global historical near-surface air temperature anomalies over land with a resolution of 5° from 1850 to 2015. To avoid biases from insufficient instrumental shielding (Böhm et al., 2010) and poor geographic coverage, we follow Wilson et al. (2016) and exclude data prior to 1880 for calibrating mid-latitudinal data and prior to 1900 for high-latitudinal data. Uncertainty may still arise from varying coverage even for later periods, which decreases toward higher latitudes and is generally higher in Western Europe and North America.

Mid Latitudes
In order to reconstruct mid latitudinal temperatures, we use data provided by the Northern Hemisphere Tree-Ring Network Development (N-TREND) consortium as published by Wilson et al. (2016) and Anchukaitis et al. (2017), containing carefully selected high quality tree-ring chronologies and local reconstructions. It includes 54 tree-ring records, from which 11 are derived from ring-width (TRW) measurements, 18 are derived from maximum latewood density (MXD) and 25 are mixed records (MIX). The latter combine reconstructions derived from different measurement methods, such as TRW, MXD, and blue intensity (BI) data. MXD and BI provide similar information on relative wood density (see Björklund et al. (2014); Rydval et al. (2014) for more information). The records vary in length, starting between AD 750 and 1710. Each record has an individual optimal target season, reflecting the trees' growing season, identified by Wilson et al. (2016), correlating the individual records against the local CRU TS 3.2 (Harris et al., 2014) mean temperature grid. The target seasons differ widely, ranging from a single month to several (Figures 1a-1c). We classify the data into two groups: early summer (sensitive up to July) and late summer (sensitive from July onwards). Early summer includes only ring-width and mixed data, while late summer includes width, density and mixed proxies. Given that total ring-width is dominated by earlywood growth, and latewood density type records use only the final cells of the latewood formation (Björklund et al., 2017), our classification is sensible also from a biological point of view.

High Latitudes
We use multiproxy data from the revised PAGES (Past Global Changes) Arctic 2k data set (v1.1.1, McKay & Kaufman (2014)), containing selected and updated records from the first PAGES Arctic 2k database (2k Consortium, 2013). In order to replicate the reconstruction method used for mid latitudinal data, we include only annually resolved data extending into the twentieth century, restricting the data set to 32 proxy records (12 ice cores, 13 tree rings, and 7 lake sediments). We refer to this data as PAGES Arctic 2k*. The data set overlaps with a few tree-ring records from N-TREND, which in some cases contains updated chronologies from PAG-ES. Based on the seasonal sensitivity stated in the original publications, the data includes 18 summer proxies (JJA), 11 annual, 2 winter (DJF), and 1 February to August proxy. Besides the complete data set, we create an additional seasonally homogeneous subset including all summer records (summer only) (Figures 1d-1f).

Proxy Reconstruction Method
We create an ensemble of hemispheric reconstructions for land surface temperature as in Lücke et al. (2019), based on iterative nesting and scaling to instrumental data (D'Arrigo et al., 2006;Wilson et al., 2016Wilson et al., , 2007. We target the mid-latitudinal band between 40 and 75° N for N-TREND, and 60-85°N for PAGES Arctic  2k*. All subsets of N-TREND are calibrated over the period 1880-1988, maximizing data availability of both proxy data and observations. Arctic 2k* subsets are calibrated between 1900 and 1967, reflecting the decreased data availability at higher latitudes. To gain an estimate of calibration and coverage uncertainty, we vary the calibration period and bootstrap the data set to eventually gain a reconstruction ensemble of 100 ensemble members per data set. In addition to the full calibration period, we split the calibration period for N-TREND into three windows of each 60, 70, and 80 years (40, 50, and 60 for Arctic 2k*). Each reconstruction is again bootstrapped by removing in turn a record from the data set (see Appendix A1 for more detail). For each individual data set, we choose the calibration target season to match the cumulative seasonal sensitivity of the included proxy records.

Hemispheric Temperature from Climate Model Data
We produce hemispheric timeseries of climate model output using three different methods. First, we use the weighted average of all model grid points within the latitudinal band and within the target season of the proxy reconstructions (CESM full band). We additionally generate a timeseries from the model data sampled at the grid points closest to the proxy locations within the target season of the proxy reconstruction (CESM proxy location). The resulting timeseries closely resembles the full band and uses seasonally homogeneous data. Last, we generate a pseudo proxy reconstruction (PPR), mimicking the heterogeneous seasonal composition of the proxy data and its reconstruction method. For each proxy record, we extract monthly climate model data at the grid point closest to its location, limited to the period at which proxy data is available. We use the monthly average within each proxy record's seasonal sensitivity window. Thus, we obtain a timeseries of model data for each proxy record that reproduces its seasonal characteristics. To concentrate on the aspect of seasonal sensitivity, we have not added noise to the pseudo proxy records, making them perfect representations of modeled climate variability. We replicate the proxy reconstruction method with the resulting pseudo proxy data sets for each model ensemble member, producing a total of 13 × 100 PPR's for each proxy target data set (see Appendix A2 for more detail).

Orbital Forcing and Millennial Trend Patterns in Climate Model Data
The change of the orbital parameters over time induces an insolation trend, which varies with month and latitude ( Figure 2). The monthly insolation trend takes the form of a propagating wave, moving forward in the year over the course of the Holocene (top). The signal is strongly zonally dependent, such that at high latitudes maximum and minimum trend are seasonally close together at all periods. During the last millennium, high northern latitudes are expected to experience an increase in insolation during April and May while insolation decreases in July and August. Thus, the magnitude of the trend changes rapidly between May and August. A similar seasonal variation in trends occurs at high southern latitudes. This finding is of particular importance, as proxy data sets are often biased toward both higher latitudes and summer temperature. Thus they directly encompass this highly variable range, as was pointed out by . This indicates that even small differences in seasonal sensitivity could lead to notable differences between temperature reconstructions.
We find a similar sinusoidal trend signal in globally averaged surface temperature of the last millennium produced by all-forced climate simulations (Figures 3a and 3b). This signal is found in all individual climate models used for our analysis ( Figure S3) as well as in CESM simulations including orbital forcing only (Figure 3c). In contrast, other forcings create a negative trend signal, which is approximately constant throughout the year (Figure 3d). The orbital only simulations are the only ones which display a warming trend, which coincides roughly with the period and region of zero trend in the all-forced simulations. Hence the detection of negative millennial temperature trends at high latitudes alone is not sufficient as attribution to orbital forcing as these could be caused by a number of different forcings. Therefore, if we wish to attribute millennial trends to orbital forcing it is the relative seasonal pattern, which is the most important pattern.
Furthermore, we observe that the temperature lags behind the insolation by around one to two months. This can be explained by the thermal inertia due to the high heat capacity of the ocean mixed layer, causing a delayed response of the surface temperature to insolation (Prescott & Collins, 1951). This mechanism also explains the larger amplitude of the seasonal trend cycle of land only data. We also note that there is a slight asymmetry in the lag between the minimum and maximum peak trends. A number of seasonally and geographically dependent mechanisms could cause an asymmetry like this (Donohoe et al., 2020). For example, the larger depth of the mixed layer in winter could lead to an increased heat capacity and thus a slower temperature increase in early spring compared to late summer (Bathen, 1972). Other mechanisms could include sea ice cover, which could induce positive or negative feedbacks depending on the season and the sign of the insolation trend (Fischer & Jungclaus, 2011).
The lag of the temperature response behind the forcing has important implications for proxy record trends, which often respond to temperature changes during a narrow seasonal window. In the case of mid-latitudinal tree-ring data, negative temperature trends would only influence proxies, which are sensitive to late summer temperatures (solid box, Figure 3). For higher latitude data, the lag would particularly impact summer proxies (dashed box). This seasonally lagged response to the forcing trend is a key finding of the study and could have important implications, particularly for previous studies, which have compared insolation trends instantaneously with trends found in proxy records (Esper, 2002;Kaufman et al., 2009;Klippel et al., 2020).

Millennial Trends in Proxy Reconstructions
We calculate the millennial trends up to AD 1850 in our proxy reconstruction ensembles, varying the start year of the fit to account for changing proxy availability and compared the results with data from CESM simulations (Figures 4a and 4c). In the N-TREND data, we observe a weaker trend for early summer and a stronger trend for late summer data, a pattern which is replicated by the model data. We conclude that the trend pattern in the latter could have been caused by orbital forcing (Figure 4b), as no other forcings LÜCKE ET AL.

(a) (b) (c) (d) (e)
produce a seasonal difference (see Section 4). The complete proxy data set shows a stronger trend than the seasonal subsets, which is not replicated by the equivalent model data. This could be a result of the data availability or reconstruction method, as the PPR does not conserve the trend difference between late summer and all summer, as observed in the model data. The PAGES Arctic 2k* data are slightly less conclusive. We observe a difference in trend between the complete data set and the summer only data. However, there is no clear orbital forcing pattern detectable in the model data for JJA versus annual temperature. While we expect the trend to be weaker in summer due to the positive orbital forcing, no clear signal can be found when comparing JJA and annual data (Figure 4d) and the model trends overlap widely. There is, however, a significant difference in trend in the proxy reconstruction. While the trend of the summer reconstruction is well in line with the model results, the trend in the complete data set, calibrated to annual data, produces a much stronger trend than found in the model. This could potentially be a result of data quality or availability, remaining seasonal sensitivity in annually calibrated data or the reconstruction method (compare the slightly inflated trends in the PPR).

Discussion and Conclusion
Our results from comparing patterns in seasonally restricted proxy reconstructions are promising for attributing seasonal trends to orbital forcing. However, both magnitude and relative difference of the trends found in model data and proxy reconstructions are subject to various uncertainties, which may have influ-LÜCKE ET AL.

(a) (b) (c) (d) (e)
enced our results. In the case of the models this includes, for example, model errors and uncertainties in the forced response and internal variability. However, all climate models agree well on the magnitude of the cooling trend as well as the sinusoidal seasonal trend pattern ( Figures S6 and S7). Small differences exist regarding the phase shift between insolation and temperature trend, but the models roughly agree on a lag of one month. This robustness of the model results provides a strong theoretical basis for a clear seasonal fingerprint of orbital forcing over the last millennium, which we should expect to see in proxy temperature records over this period.
In the case of the proxy results, we cannot exclude the possibility that the trend pattern in the proxy data is caused by methodological or spectral proxy biases. Methodological biases could arise from the calibration, for example, weaker correlation with instrumental data could suppress the trend (which could be an issue in particular for N-TREND early summer data), but also from detrending methods. It is known that certain detrending methods do not preserve low-frequency variability (Melvin et al., 2012), which could induce a flat bias in the reconstruction. This could in particular influence N-TREND late summer (includes two such records: IDA, KOL, see Table S1), but also early summer (also includes IDA). On the other hand, trends could also be overestimated following detrending -a known issue for at least one record (QUEw) included in N-TREND all and late summer (Anchukaitis et al. (2017)). Furthermore, spectral biases inherent to the proxy data could arise from biological memory processes and overestimate the trend (Lücke et al., 2019). We also cannot exclude the possibility that spectral differences exist between ring-width and latewood density as suggested by Esper et al. (2012), and that our results simply express the characteristics associated with LÜCKE ET AL.
10.1029/2020GL088776 8 of 13 (d) (c) these datatypes. Limited data availability of the seasonal data sets and potential data quality issues, as well as uncertainties in the seasonal sensitivity could further reduce confidence.
The interpretation of the results is increasingly complicated for inhomogeneous proxy data sets. We notice that for both the N-TREND and especially for the PAGES data set, the trend of the reconstruction including all data differs more from the model than for the seasonally homogeneous data subsets. This applies in particular to the PAGES data set, which includes multiproxy data whose seasonal sensitivity vary strongly. While it is thought that ice-cores and marine/lake sediments may be better than tree rings at capturing low-frequency variability , the summer reconstruction (tree-rings and lake sediments only) seems to be more in line with the model results. This could reflect how different seasonal proxy characteristics complicate the interpretation of results from multiproxy data, and accordingly decrease our confidence in the results.
We conclude that orbital forcing has likely played a key role in long-term seasonal forcing during the last millennium. While other forcings cause negative millennial temperature trends throughout the year, orbital forcing is the only one inducing a characteristic seasonal fingerprint that may lead to differential millennial cooling across seasons. Thus, negative trends in regional high latitude summer proxy data alone cannot be attributed to orbital forcings as other forcings may also cause negative trends during that time. The trends detected in modeled surface temperature lag behind insolation trends by around a month. Since the period of tree-ring formation encompasses both months of positive and negative orbital trends, trends in tree-ring data will likely depend on their specific seasonal sensitivity. It is therefore critical when interpreting longterm trends in proxy reconstructions to consider the specific seasonal sensitivity of the individual proxy records as well as the lagged relationship between the forcing and the temperature response.
Our analysis shows that climate models and proxy data show similar seasonal trend patterns, which can, in the case of the model data, be attributed to orbital forcing. However, a lack of clearly seasonally resolved, high-quality proxy data makes it impossible to reliably detect orbital forcing in last millennium proxy reconstructions at present, even though differential trends in proxy data suggest an orbital influence.
Our findings could resolve previous data discrepancies, including the difference in long-term trends found in ring width and density data (Esper et al., 2012;Fuentes et al., 2018;Klippel et al., 2020). It also provides an explanation for differing low-frequency variability between proxy reconstructions and model simulations (e.g. the difference between MCA and LIA). We strongly advise against the mixing of proxy data with unknown or differing seasonal sensitivities when it comes to reconstructing long-term trends, or focus on multidecadal variability for such reconstructions (e.g. Neukom et al. (2019)).

Conflict of Interest
The authors declare no conflict of interests.

Data Availability Statement
All model data are publicly available: PMIP/CMIP model output is available via the Earth System Grid at: http://pcmdi9.llnl.gov/. The CESM data can be downloaded from: http://www.cesm.ucar.edu/projects/community-projects/LME/. The CSIRO Mk3L data from: https://www.ncdc.noaa.gov/paleo-search/ study/16337. Both proxy data sets are publically available. The N-TREND data set can be downloaded at: https://ntrenddendro.wordpress.com/tr-data/. The PAGES Arctic 2k v1.1.1 can be downloaded from: https://www.ncdc.noaa.gov/paleo-search/study/16973. All data generated or analyzed during this study will be publically available on the University of Edinburgh DataShare server. Acknowledgments L. Lücke. was supported by a studentship from the Natural Environment Research Council (NERC) E3 Doctoral training partnership (grant number NE/ L002558/1). A. P. Schurer and G.Hegerl. were supported by NERC under the Belmont forum, Grant PacMedy (NE/ P006752/1). The authors acknowledge the World Climate Research Program's Working Group on Coupled Modeling, which is responsible for CMIP, and thank all the climate modeling groups for producing and making available their model output. The authors acknowledge the Northern Hemisphere Tree-Ring Network Development (N-TREND) and the Past Global Changes (PAGES) project for providing publicly available data.