Complex spatio-temporal structure of the Holocene Thermal Maximum

Inconsistencies between Holocene climate reconstructions and numerical model simulations question the robustness of climate models and proxy temperature records. Climate reconstructions suggest an early-middle Holocene Thermal Maximum (HTM) followed by gradual cooling, whereas climate models indicate continuous warming. This discrepancy either implies seasonal biases in proxy-based climate reconstructions, or that the climate model sensitivity to forcings and feedbacks needs to be reevaluated. Here, we analyze a global database of Holocene paleotemperature records to investigate the spatiotemporal structure of the HTM. Continental proxy records at mid and high latitudes of the Northern Hemisphere portray a “classic” HTM (8–4 ka). In contrast, marine proxy records from the same latitudes reveal an earlier HTM (11–7ka), while a clear temperature anomaly is missing in the tropics. The results indicate a heterogeneous response to climate forcing and highlight the lack of globally synchronous HTM.

Inconsistencies between Holocene climate reconstructions and numerical model simulations question the robustness of climate models and proxy temperature records. Climate reconstructions suggest an early-middle Holocene Thermal Maximum (HTM) followed by gradual cooling, whereas climate models indicate continuous warming. This discrepancy either implies seasonal biases in proxy-based climate reconstructions, or that the climate model sensitivity to forcings and feedbacks needs to be reevaluated. Here, we analyze a global database of Holocene paleotemperature records to investigate the spatiotemporal structure of the HTM. Continental proxy records at mid and high latitudes of the Northern Hemisphere portray a "classic" HTM (8-4 ka). In contrast, marine proxy records from the same latitudes reveal an earlier HTM (11-7ka), while a clear temperature anomaly is missing in the tropics. The results indicate a heterogeneous response to climate forcing and highlight the lack of globally synchronous HTM. Natural climate variability results from multiple forcings and feedbacks with heterogenous spatiotemporal manifestations. Greenhouse gases, volcanic radiative forcing, and solar irradiance apply rather homogeneously across the Earth's surface, while insolation varies both latitudinally and seasonally. In addition, the climate system response may be amplified or dampened by feedbacks inherent to changes in physiography, albedo, and by variations in oceanic and/or atmospheric circulation that (re)distribute heat across the Earth's surface. Our understanding of climate processes is limited by the rather short temporal span and heterogenous spatial coverage of instrumental records. Evidence of past climate variability gleaned through the testimony of geological archives thus offers a unique opportunity to contextualize ongoing changes and to assess climate model performance on timescales going beyond the decadal climate variability recorded in the instrumental period.
The temperature at the Earth's surface responds directly to global radiative forcing and thus provides fundamental insights into the state of the climate system. Over the past decades, quantitative indicators of past temperature (hereafter called "proxies") based on different types of archives have been used to reconstruct climate variability over a range of timescales. The improvement of both spatial coverage and temporal resolution of temperature proxy records led to the development of regional and global temperature reconstructions, which have allowed the scientific community to highlight the unprecedented nature of anthropogenic climate change across the common era 1,2 and the Holocene 3-6 . Global temperature reconstructions consistently depict a Holocene Thermal Maximum (HTM) typically ranging between 10 and 5 ka 4,5 with a maximal probability centered around 6.45 ka 4 . The HTM was followed by global cooling until the end of the nineteenth century CE, interrupted by rapid and sustained warming characterizing the industrial era towards the present. Yet, the cooling trend inferred from proxy records, often attributed to declining high northern latitude insolation, cannot be resolved in numerical simulations 7 . Indeed, in climate models, the simulated global mean temperature is predominantly driven by the ice-sheet extent and atmospheric greenhouse gas concentrations, which in synergy impose continuous warming over the course of the Holocene 7 .
This discrepancy between proxy data and model simulations, commonly referred to as "The Holocene Temperature Conundrum" 7 , casts doubt on the conceptual framework underlying temperature proxy interpretation and on climate model skill. For instance, it has been suggested that temperature reconstructions may be seasonally biased 7,8 and/or that the global mean value is skewed because of the overrepresentation of northern North Atlantic sea-surface temperature (SST) records [5][6][7] . However, model-data inconsistencies may equally well result from geographically divergent trends due to sea-ice dynamics 9 , polar amplification 10 , insufficient model resolution 11 , and boundary conditions used in numerical simulations 12 . Although the HTM has been intensively studied from a global perspective 3-7 , its spatio-temporal characteristics have received relatively little attention, even though the local and regional trends differ markedly from the globally averaged reconstructions 3,13 .
In this study, we seek to document the spatiotemporal expression of the HTM in the marine and continental realms to shed light on the forcings and feedbacks underpinning the evolution of Holocene climate 9 .

Temperature 12k database analyses
In order to investigate the spatiotemporal expression of the HTM, we analyze the Temperature 12k compilation 3 , which includes 1319 globally distributed paleotemperature records with an average sampling and age control resolution of 164 and 1000 years, respectively. We initially select records that cover at least 80% of the Holocene. This subset includes 233 (184 annual; 32 summer and 17 winter) marine records and 470 (159 annual; 184 summer and 127 winter) continental records (see details in "Methods"). The terrestrial subset is dominated by pollen-based temperature reconstructions from the Northern Hemisphere, with virtually all records (91%) distributed between 40°N and 70°N (Fig. 1a). The SST subset is more widely distributed, although somewhat concentrated towards continental margins and the Northern Hemisphere (Fig. 1b, see details in "Methods").
The Temperature 12k compilation provides age and temperature ensembles based on dated intervals and proxy calibration (100 members each) for 207 of the selected marine records and for 462 continental records. We use these to assess the robustness of our inferences against reconstruction and chronological uncertainties. We generate 10,000 time series based on the age and temperature ensembles for each record, identify the age of the highest temperature, and calculate the temperature anomaly compared to the mean Holocene value (see details in "Methods"). We subsequently analyze the temporal and spatial distribution of the selected local HTM age ensembles using one-and two-dimensional kernel probability density function (Figs. 1-7, Gaussian kernel, with a bandwidth of 500 years and 10°in latitude, see details in "Methods").

Continental HTM
Continental temperature reconstructions between 40°N and 70°N most often reach their maximum values between 8 and 4 ka, with a slight dependence on proxy type (Figs. 1, 5, and 7). The median HTM age based on pollen records (~5.7 ka) lags continental temperature reconstructions based on chironomids (~7.3 ka) by~1600 years (Fig. 7). Notwithstanding proxy type, the integral of the probability density functions of continental HTM ages between 4 and 8 ka amounts to 40.4%, while representing only 33.3% of the Holocene time interval. In other words, 60% of the records reached maximum temperatures either before or after the 8-4 ka interval, highlighting a dispersion of local HTM ages compared to global reconstructions. The HTM occurs at around 7 ka across North America, but rather later (5 ka) across western Europe, explaining the broad bimodal distribution characterizing the continental records (Fig. 6b). Many of the records are considered to reflect mean annual temperature, while others are assumed to reflect seasonal conditions. However, the probability density function (PDF) of HTM age calculated for winter, summer and annual temperature reconstructions are virtually indistinguishable.

Marine HTM
Considered globally, marine records do not display a pronounced clustering of HTM ages (Fig. 1). However, the global average eclipses marked spatial variability in HTM timing. Whereas the tropics conform with the global marine reconstruction (30°S-30°N), time series from the mid-high latitudes of the Northern Hemisphere clearly precede the global average (9.7 ka), thus leading the NH continental temperature optimum by~3500 years ( Figs. 1 and 5). Although records south of 30°S remain sparse, they also suggest an early HTM (12-10 ka), hinting at hemispheric symmetry (Fig. 6). As for terrestrial records, the HTM age distributions based on annual and seasonal proxies are  Spatiotemporal structure of the Holocene Thermal Maximum (HTM) derived from marine and continental temperature records. Dots represent median HTM age and median HTM amplitude (anomaly compared to mean Holocene temperatures) calculated from the age and temperature ensembles for each continental (A) and marine record (B). Gray shadings and contours depict the mean of 10,000 kernel probability density function maps (PDF; Gaussian kernel with a bandwidth of 10°in latitude and 500 years) calculated using the temperature and age model ensembles (see also Fig. 5 for contour label). Left and bottom panels correspond to median kernel probability density calculated against latitude (10°bandwidth) and age (500 years bandwidth), respectively for different subsets, with dashed straight lines corresponding to a homogenous distribution (color shading shows 10-90 percentiles).
undistinguishable, yet the number of winter temperature time series remains insufficient (n = 15) to derive statistically robust conclusions (Fig. 6). The HTM anomaly mostly remains below 6°C and shows a clear meridional pattern with higher temperature amplitudes in the mid latitudes of both hemispheres (Figs. 1, 2, and 7).

Latitudinal pattern in the expression of the marine Holocene Thermal Optimum
The spatiotemporal structure of the HTM in the marine realm indicates that Holocene temperature trends vary by latitude, with early (11-7 ka) and high amplitude HTM manifestations at mid to high latitudes of the Northern Hemisphere, and muted expression in the tropics, in accordance with previous studies 3,5,13 . The boundary between these climate regimes is located around 30°N (Fig. 1b).
The latitudinal pattern in the timing of marine HTM is most likely steered by orbital forcing, because it is the largest contributor to radiative forcing, especially during the early Holocene (Figs. 3, 8-10). However, insolation varies with both latitude and season. The response of the climate system as recorded in proxy timeseries could hence be regionally and temporally variable.   During the Early Holocene mean annual insolation was significantly higher than at present, by up to 5 W/m 2 near the poles and 2.5 W/m 2 at 65°(Figs. 3 and 10) but was lower at mid to low latitude (−1 W/m 2 at the equator). It is therefore conceivable that the zonal trend portrayed by the marine temperature records reflects changes in mean annual insolation. However, forcing with almost similar amplitude in the tropics during late Holocene (2 W/m 2 ) due to increasing greenhouse gas concentrations and insolation, was not associated with a clear temperature maximum (Fig. 3). This makes changes in the mean annual forcing alone an unlikely control on the HTM expression. Instead, latitude-dependent mechanisms likely amplified radiative forcing north of 60°N which impacted climate southward to 30°N as a result of polar amplification feedbacks due to sea ice albedo 14 . However, changes in seasonal insolation might offer a more direct explanation of the observed early marine HTM at mid latitudes.
The sensitivity of local mean annual temperatures to seasonal insolation is nonlinear and varies spatially, with most of the tropical to subpolar Earth surface, between 60°S and 60°N, predominantly influenced by summer insolation 15 . During the Holocene, the maximum summer insolation in the Northern Hemisphere occurred between 10 and 7 ka with anomalies of more than 25 W/m 2 compared to modern (Fig. 9). This is consistent with the inferred early HTM timing in extratropical regions of the Northern Hemisphere (Fig. 3).
The similarity between summer, winter and annual HTM in the marine records evidenced here, despite contrasted seasonal insolation trends, may shed light on the climatic feedbacks affecting local HTM manifestation. While some numerical simulations indicate strong dependency of monthly temperature on monthly insolation, implying divergent seasonal trends over the course of the Holocene 10,16 , other models indicate year-long effects of summer insolation on annual temperature in the ocean 10 at high-to mid-latitudes. The latter display strong polar amplification, annual temperature dependency on summer sea ice loss and better agreement with annual 10 and seasonal proxy data (this study), offering a plausible solution to the HTM conundrum. However, it relies on seasonality attributed to proxies which remains imperfectly resolved and debated 8,17,18 . The PDF of the arguably sparser Southern Hemisphere records shows a dual HTM (early and late, Fig. 6) with low confidence. If the interhemispheric symmetry of an early Holocene marine HTM is real despite the asymmetrical seasonal forcing, it must either reflect contrasted sensitivity to seasonal insolation at high latitudes 15 (potentially due to sea ice dynamics and/or continental configuration), indicate the dominance of high latitude annual insolation forcing over summer insolation on Holocene climate (see above) or suggests the strong global influence of glacier retreat.
The early HTM peak in the Southern Hemisphere probably also relates to changes in the Atlantic Meridional Oceanic Circulation (AMOC) which affected interhemispheric heat transport 13 through the bipolar seesaw 19 across the deglaciation, leading to early interglacial conditions (12 ka) in the Southern Hemisphere 20 . The late HTM peak in the SH (Fig. 6) may result from proxy or sites more sensitive to the austral summer insolation that peaked late in the Holocene (Fig. 9). Clearly, a better spatiotemporal coverage of paleotemperature records from the Southern Hemisphere is required to confidently evaluate the relative importance of the aforementioned processes on   Holocene climate, in particular the impact of annual and seasonal insolation changes.

Delayed continental timing of the Holocene Thermal Optimum
Our analysis demonstrates that the HTM occurred earlier in the ocean than on land at the mid latitudes of the Northern Hemisphere, suggesting that different mechanisms and feedbacks drive temperature changes in oceanic and continental realms, or that contrasting biases affect marine and continental temperature proxies. Continental HTM occurrences between 30 and 70°N peak between 4 and 8 ka, approximately 3.5 kyr later than in the ocean at the same latitudes. Assuming that both marine and continental temperature proxies are not seasonally biased, the comparatively early HTM in the ocean could simply relate to a higher sensitivity of marine temperatures to summer insolation at mid latitudes 15 . In general, the temperature response to seasonal insolation at high latitudes differs in the ocean and on land because of the high thermal inertia of seawater and limited sensitivity to winter insolation due to vertical mixing 15 , not to mention the minimum winter temperature threshold of −1.8°C that corresponds to the seawater freezing temperature. Hence, the oceancontinent difference could result from a systematic seasonal bias towards summer in the high latitude marine records. However, because of seasonal variability in the abundance of the proxies, both foraminifera-18 and alkenone-based temperature reconstructions 21 are more likely to reflect annual mean conditions at mid-latitudes, whereas pollen-based temperatures are more likely to be biased toward the flowering season (summer) 22 . Because both summer and annual insolation in the Northern Hemisphere peak early in the Holocene, seasonal bias in the proxy records is unlikely to explain the ocean-land difference in the HTM timing.
Marine and continental climate would have also responded differently to climatic feedbacks such as glacier fluctuations, snow cover and vegetation changes 23 . The Eurasian ice sheet disappeared at around 8 ka 24 , while the final deglaciation of the Laurentide ice sheet occurred later at 6.7 ka 25 . The terrestrial HTM coincided with the stabilization of sea level (Fig. 3) when ice sheets receded to their modern limits, suggesting that dwindling continental ice sheets may have delayed warming on land. Although, the HTM occurred slightly earlier close to the Laurentide ice sheet than in Western Europe (Fig. 6), highlighting the complex and spatially heterogenous climatic response to dwindling ice sheets 12,26 .
The generally younger HTM age on land may also partly relate to a specific bias affecting pollen-based temperature estimates. Amongst the continental temperature records, pollen-based records are The category named "all" includes every records while the"N30°N" category includes all records north of 30°N. The mean (white dot) and median (colored horizontal bar) are displayed for each subset.
dominant, and display a generally younger HTM age distribution than the second most abundant continental reconstructions based on chironomid assemblages (Fig. 7). Vegetation changes may not only be modulated by temperature, but also by the hydrological regime, soil development, species migration modes, fire history and anthropogenic pressure among others. For instance, glacial erosion affects the soil formation, particularly at high latitude 27 , leaving a long-lasting imprint of glacial climate on interglacial vegetation 28 . As such, the northward migration and expansion of vegetation during the deglaciation possibly led to a transient disequilibrium between the local vegetation composition, pollen assemblages and climate during the early Holocene. Moreover, changes in vegetation cover at mid latitudes during the deglaciation and the Holocene might have impacted the continental climate through changes in albedo 23,29 and may have contributed to the time lag with respect to the ocean.

Latitudinal pattern in the amplitude of the marine HTM
The marine HTM anomaly displays a latitudinal pattern with higher amplitudes centered around 45°S, 45°N and 65°N (Fig. 2). This pattern bears similarity with the pattern in modern SST gradients and the position of oceanic frontal zones (Fig. 2). Today, steep SST gradients occur at the poleward boundary of the oceanic subtropical gyres, around 45°of latitude both in Southern and Northern Hemispheres, and at 65°N, north of the subpolar North Atlantic gyre, close to the gateway between the North Atlantic and the Arctic Ocean (Fig. 2). Thus, the large Holocene temperature anomalies recorded at similar latitudes as modern frontal zones, may reflect enhanced local temperature sensitivity to latitudinal migrations of oceanic fronts and suggests that surface ocean dynamics played a significant role in the recording of the spatiotemporal evolution of Holocene temperatures. High-amplitude variations in SSTs near western boundary currents and frontal zones have been previously recognized 30 . Indeed, the inclusion of records from frontal zones in global reconstructions has been challenged due to their high sensitivity to shifts in surface currents 6,8 . We provide an objective analysis of these features which indicate that a high amplitude HTM is not restricted to western ocean basins, but also affects eastern flanks (Fig. 2). A southward shift of the Gulf Stream and Kuroshio may explain high-amplitude temperature trends near western boundary currents 30 , although an equatorward contraction of the subtropical oceanic gyres is a better candidate in eastern oceanic basins where the spatial temperature gradient is low. Moreover, reduction of the Atlantic Meridional Overturning (AMOC) strength from mid to late Holocene 31 would have also contributed to cooling 32 . Irrespective of the exact cause of spatial discrepancies, our analysis highlights the importance of small-scale ocean dynamics for the evolution of Holocene seawater temperature. Importantly, numerical simulations from coupled model experiments of mid-Holocene climate exhibit a large dispersion in simulated temperature at frontal zones 33 , likely due to insufficient spatial resolution. This disagreement among simulations corroborates the critical role of oceanic circulation, and the relevance of exploring ocean front dynamics further for simulating long term climate change and addressing the HTM conundrum in more detail.

Complex Holocene temperature patterns
Our analyses reveal a dynamic and spatially variable evolution of Holocene temperature that contrasts with the notion of a globally synchronous mid Holocene thermal optimum. The observed spatiotemporal structure of the HTM resulted from geographically divergent temperature trends likely related to differential seasonal and spatial responses to climate forcing and feedbacks. The spatiotemporal dynamics should be considered in data-model comparison for the mid Holocene time slice at 6 ka 34 , in particular, the delayed continental HTM on land compared to the ocean, and the latitudinal contrast in the timing and amplitude of the marine Holocene Thermal Maximum we identified here. Proxy records suggest similar annual and seasonal temperature trends both on land and in the ocean across the Holocene. Seasonal trends in the ocean better conforms with numerical simulations displaying a high degree of Arctic amplification due to sea ice loss, which might reconcile proxy-based reconstructions and numerical simulations 10 . Importantly, the comparison between global mean temperature reconstructions and numerical simulations, at the crux of the HTM conundrum, has limited mechanistic implications. Instead, the spatiotemporal structure of the temperature trends over the Holocene bears more valuable information on the relative importance of the forcing mechanisms and feedbacks acting on climate. The high amplitude temperature variations in frontal zones highlight the importance of oceanic circulation on Holocene climate variability and suggest that model-data (dis)agreement partly hinges on resolving frontal dynamics in climate simulations, in addition to the degree of polar amplification. Gaining confidence in the regional aspects of future climate projections 35 will therefore partly rely on a better representation of oceanic frontal zones in numeric simulations, but also on improved proxy coverage on land outside the 30°−70°N latitudinal band and in the southern Hemisphere ocean, as well as on development of seasonal proxy records.

Data description
We analyze the Temp12k compilation described in ref. 3, to identify age and amplitude of the local HTM and assess uncertainties in reconstructed local HTM characteristics.
In order to avoid bias due to uncomplete coverage of individual records, we first identify records with sufficient coverage by binning individual records into 500-year bins based on their original age model. We choose 500 years because of the minimum 400-year resolution criterion for the dataset selection 3 and because the broadness of the bins mitigates the influence of chronological uncertainties. The records used in our analyses have at least one value in 80% of the twenty-four 500 year long bins between 0 and 12 ka, thus 20 of the 24 bins. Selected records include 233 (184 annual; 32 summer and 17 winter) marine records and 470 (159 annual; 184 summer and 127 winter) continental records (Fig. 4). The Temp12k compilation provides age and temperature ensembles (100 realizations each) for 207 of the previously selected marine records and for 462 continental records. Age ensembles were obtained using Bayesian age-depth modeling approach implemented in Bacon 3,36 . Temperature ensembles for pollen data are derived from the Modern Analogue Technique applied by 6 while temperature ensemble records from marine sediment are based on Bayesian procedures [37][38][39][40] or on multiple generations of analytical and calibration methods (see details in ref. 3). When temperature uncertainties ensemble were not available, we added noise from a standard Gaussian distribution scaled to the temperature uncertainty attributed to each record by the Temp12k consortium to the raw reconstruction (uncertainty on temperature ranges from ±1 to ±3°C depending on the proxy type according to 3 ) and note that these estimates are rather on the conservative side.
In order to characterize the local HTM patterns and to quantify uncertainties associated to ages and temperature reconstruction, we combine all age and reconstruction ensemble members to derive 10,000 realizations of each temperature time series.

Local HTM analyses
We then identify the maximum temperature for each realization, and obtain its age and amplitude with respect to the Holocene mean. We then use the spread among the 10,000 realizations for each record to quantify the uncertainty in characteristics of the structure of the local HTM.
To analyze the temporal and spatial distribution of the timing of the HTM we use a two-dimensional kernel probability density function (Fig. 5, Gaussian kernel, with a bandwidth of 500 years and 10°in latitude). We calculate the kernel density maps 10,000 using only one realization per records for each iteration (mean shown in Fig. 1). The mean and standard deviation of the resulting two-dimensional PDFs allow to control the robustness of the results and to account for uncertainties in both age and temperature reconstruction. The standard deviation of the PDF remains small relatively to mean PDF values in highest density domain (below 20 and 30% of local PDF for continental and marine records respectively, Fig. 5), suggesting that spatiotemporal patterns discussed in the main text are robust. We also tested other kernel functions (Epanechnikov, box (uniform), triangle) but found no significant differences at the considered bandwidth, suggesting that uncertainty in the data has a larger effect than the choice of kernel. We therefore used a Gaussian kernel with a bandwidth of 500 years in all analyses.
To determine the effect of seasonality and location we repeated our analyses on subsets of the data using one-dimensional kernel probability density function (Gaussian kernel, with a bandwidth of 500 years, Fig. 6). These analyses highlight the robustness of the difference in the timing of the HTM between the continent and the ocean at mid Northern hemisphere latitudes. They also show the effects of seasonality and location in the continental records discussed in the main text.
We assess the effect of proxy sensors on the distribution of HTM ages and amplitude using kernel density-based violin plots (Fig. 7). Despite the large spread among the proxies (which partly reflects regional differences) the early timing of the high latitude marine temperature maxima appears in these plots. Although mean and median HTM age rather occurs around the mid Holocene for both continental and marine environments, marine HTM display a skewed distribution towards the late Holocene, and the mode of the marine HTM often occurs a few thousands of years before the mean and median, and does not result from sensor-specific biases.
Almost all proxy types display few anomalously high amplitudes, which are probably the result of the effect of extraneous variables on temperature estimates. We however have not excluded those from our analyses as their number is too low to significantly affect our results.

Climate forcing
We show the forcings that apply relatively homogenously over the earth surface (volcanic radiative forcing 41 , decadal solar irradiance changes 42 , and greenhouse gases atmospheric concentrations forcing including pCO 2 43 , pCH 4  term sum used in Fig. 3. Note that pCO 2 , pCH 4 and pN 2 O radiative forcing (w/m 2 ) were calculated relatively to preindustrial forcing accounting for the overlap in absorption band between N 2 O and CH 4 45 . Seasonal and monthly insolation anomalies compared to modern insolation for the last 12 ka are displayed in Fig. 9, and mean annual insolation is shown in Fig. 10.

Code availability
Computer codes are available upon request from O.C.