Top-of-atmosphere radiation fields over the last millennium reconstructed from proxies

TL;DR

This study tackles the challenge of quantifying Earth's energy budget variability over the last millennium by reconstructing gridded, seasonal top-of-atmosphere (TOA) radiation fields and upper-ocean heat content (OHC) from proxies, with clear separation of shortwave and longwave components. The authors implement an online data-assimilation framework that blends proxy information with climate-model dynamics using linear inverse models (LIMs), an ensemble Kalman filter (EnKF), and proxy system models (PSMs), trained on CMIP6 last-millennium simulations to produce three model-prior reconstructions and a 1200-member multi-model ensemble. Key findings reveal a millennium-scale cooling trend with intermittent energy gains, a southwestward shift of Indo-Pacific convection and sea-ice growth, and a volcanic-forcing–driven OHC loss followed by 5- to 10-year recovery; the recent TOA energy imbalance after 1980 is unprecedented in the pre-industrial era. The work provides a spatially and seasonally resolved TOA-radiation and EEI record that disentangles SW/LW processes, offering a robust benchmark for understanding millennial climate variability and evaluating climate-model performance, while highlighting limitations related to proxy sparsity and historical-period partitioning.

Abstract

Earth's energy imbalance at the top of the atmosphere is a key climate system metric, but its variability is poorly constrained by the short observational record and large uncertainty in coupled climate models. While existing ocean heat content reconstructions offer a longer perspective, they cannot separate the contributions of shortwave and longwave radiation, obscuring the underlying processes. We extend the energy budget record by reconstructing the top-of-atmosphere radiation and related surface variables over the last millennium (850--2000 CE). Our method combines proxy data and model dynamics using seasonal, gridded data assimilation, exploiting the covariance of radiation with proxies sensitive to surface temperatures. The method validates well in a pseudoproxy experiment and against instrumental observations. We find that a last-millennium cooling trend coincides with heat loss that gradually slowed, although there are intermittent multidecadal periods of energy gain. The cooling trend is associated with a southwestward shift of Indo--Pacific convection and growth of sea ice, with seasonal sea ice trends following orbital-scale changes in polar insolation. We also find that the upper-ocean heat content following large volcanic eruptions does not begin to recover until 5--10 years later, suggesting the initiation of the Little Ice Age by decadally-paced eruptions in the early 1100s and late 1200s. Indeed, the latter period marks the decade of largest energy loss over the last millennium. Our reconstruction reveals that the energy imbalance for all 20-yr periods after 1980 is unprecedented in the pre-industrial period.

Figures (10)

• Figure 1: Annual-mean correlation between the SST at a point in the Niño 3.4 region (yellow pentagon; e.g., location of a coral proxy) and the TOA (a) reflected SW and (b) outgoing LW radiation globally. The correlation is from last-millennium simulations over 850--1850, averaged over four models.
• Figure 2: Comparison of the pseudoproxy imperfect-model reconstruction with the truth simulation. (a--d) Global-mean time series and their correlation coefficients r over 850--1850. The reconstruction (black) and the truth (yellow) are anomalies with a 20-yr low-pass filter. The truth anomalies are relative to 1961--1990, and the reconstruction has been shifted to match the truth in the 850--1850 mean. DA is the mean of the three single-model reconstructions with the MPI, CESM, and MRI model priors. Shading denotes the 5th--95th percentile range. (e--h) Spatial correlation of annual anomalies over 850--1850. Numbers in the top right corners represent global-mean values.
• Figure 3: Annual TOA radiation anomalies from our reconstruction, the DEEP-C combined product, and ERBE satellite measurements. The anomalies are relative to 1985--1999. Shading denotes the 5th--95th percentile range. (a--c) Global-mean time series and their correlation coefficients r. The p-values indicate if r is different from zero. (d--f) Spatial correlation with DEEP-C. Numbers in the top right corners represent global-mean values.
• Figure 4: Annual upper-ocean heat content from our reconstruction and the Cheng2017 instrumental dataset. The linear trend over 1960--2000 has been removed to emphasize interannual to decadal variations as in Church2005. Stratospheric aerosol optical depth (SAOD) at 550 nm is from the CMIP7 historical volcanic forcing dataset Aubry2025.
• Figure 5: Global-mean time series over the last millennium. DA refers to the combined ensemble of the multi-model reconstructions. Shading denotes the 5th--95th percentile range. Bold lines have a 20-yr low-pass filter, thin gray lines are annual values. (a) SAT anomalies relative to 1961--1990, comparing to the multi-method ensemble and median from PAGESConsortium2019 and the HadCRUT5 instrumental dataset Morice2021. Carets indicate volcanic eruptions, scaled by their volcanic stratospheric sulfur injection, from the eVolv2k Toohey2017 and the CMIP7 historical volcanic forcing datasets (Aubry2025; small eruptions with a VSSI below 0.5 Tg S removed). Those marked in orange are composited in Fig. \ref{['fig:composite_volcano']}. (b) Absolute energy imbalance. The DA EEI anomalies are made absolute by shifting to a zero-mean during 1600--1900, chosen to coincide with a period of relatively constant temperature. Colored lines show the rate of OHC change from the Gebbie2019 and Wu2025a reconstructions, expressed as equivalent EEI. The early 1100s and late 1200s volcanic clusters are annotated. (c) RSR and OLR. Atmosphere-only simulations (amip-hist; Zhou2016a), averaged over seven models, are shown in dashed lines as anomalies relative to 1961--1990. The DA anomalies in (c) are shifted to match the simulation 1871--1900 mean.