# Globally coherent water cycle response to temperature changes during the past two millennia

Bronwen Konecky, Guillaume Leduc, Marie-Alexandrine Sicre, Emilie Pauline Dassié

► To cite this version:

Bronwen Konecky, Guillaume Leduc, Marie-Alexandrine Sicre, Emilie Pauline Dassié. Globally coherent water cycle response to temperature changes during the past two millennia. Nature Geoscience, 2023, 10.1038/s41561-023-01291-3. hal-04301540

HAL Id: hal-04301540

https://hal.science/hal-04301540v1

Submitted on 23 Nov 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Title: Globally coherent water cycle response to temperature changes during the past two millennia

**Authors:** Bronwen L. Konecky¹*, Nicholas P. McKay², Georgina M. Falster¹,³, Samantha L. Stevenson⁴, Matt J. Fischer⁵, Alyssa R. Atwood⁶, Diane M. Thompson⁷, Matthew D. Jones⁸, Jonathan J. Tyler⁹, Kristine L. DeLong¹⁰, Belen Martrat¹¹, Elizabeth K. Thomas¹², Jessica L. Conroy¹³, Sylvia G. Dee¹⁴, Lukas Jonkers¹⁵, Olga V. Churakova (Sidorova)¹⁶,¹⁷, Zoltán Kern¹⁸,¹⁹, Thomas Opel²⁰, Trevor J. Porter²¹, Hussein R. Sayani²², Grzegorz Skrzypek²³, and Iso2k Project Members†

## Affiliations:

¹Department of Earth and Planetary Sciences, Washington University, Saint Louis, Missouri, 63108, USA

²School of Earth and Sustainability, Northern Arizona University, Flagstaff, AZ, 86011, USA

³Research School of Earth Sciences and ARC Centre of Excellence for Climate Extremes, Australian National University, Canberra, ACT, 2601, Australia

⁴Bren School of Environmental Science & Management, University of California, Santa Barbara, Santa Barbara, CA, 93106, USA

⁵NSTLI Environment, ANSTO, Sydney, NSW, 2234, Australia

⁶Department of Earth, Ocean, and Atmospheric Sciences, Florida State University, Tallahassee, Florida, 32306, USA

⁷Department of Geosciences, University of Arizona, Tucson, Arizona, 85719, USA

⁸School of Geography, University of Nottingham, Nottingham, NG7 2RD, UK

⁹Department of Earth Sciences, The University of Adelaide, Adelaide, South Australia, 5005, Australia

¹⁰Department of Geography and Anthropology, Coastal Studies Institute, Louisiana State University, Baton Rouge, LA, 70803, USA

¹¹Department of Environmental Chemistry, Institute of Environmental Assessment and Water Research (IDAEA-CSIC), Barcelona, Barcelona, 08034, Spain

¹²Department of Geology, University at Buffalo, Buffalo, NY, 14260, USA

¹³Department of Geology, Department of Plant Biology, University of Illinois at Urbana-Champaign, Urbana, IL, 61822, USA

¹⁴Department of Earth, Environmental, and Planetary Sciences, Rice University, Houston, Texas, 77005, USA

¹⁵MARUM Center for Marine Environmental Sciences, Bremen University, Bremen, 28359, Germany

¹⁶Institute of Ecology and Geography, Siberian Federal University, Krasnoyarsk, 660041, Russian Federation

¹⁷Department of Forest Dynamics, Swiss Federal Institute for Forest, Snow and Landscape Research WSL, Birmensdorf, CH-8903, Switzerland

¹⁸Institute for Geological and Geochemical Research, Research Centre for Astronomy and Earth Sciences, Eötvös Loránd Research Network, Budapest, H-1112, Hungary

¹⁹CSFK, MTA Centre of Excellence, Budapest, H-1121, Hungary

²⁰Polar Terrestrial Environmental Systems, Alfred Wegener Institute Helmholtz Centre for Polar and Marine Research, Potsdam, 14473, Germany

²¹Department of Geography, Geomatics and Environment, University of Toronto - Mississauga, Mississauga, Ontario, L5L1C6, Canada

²²School of Earth and Atmospheric Science, Georgia Institute of Technology, Atlanta, GA, 30332, USA

²³West Australian Biogeochemistry Centre, School of Biological Sciences, The University of Western Australia, Perth, Western Australia, 6009, Australia
† Iso2k Project Members includes all above named authors and: Nerilie J. Abram (Australian National University, Australia), Kerstin Braun (Arizona State University, USA), Matthieu Carré (Sorbonne Universités (UPMC)-CNRS-IRD-MNHN, France), Olivier Cartapanis (Aix-Marseille University, France), Laia Comas-Bru (University of Reading, United Kingdom), Mark A. Curran (Australian Antarctic Division, Australia), Emilie P. Dassié (University of Bordeaux, France), Michael Deininger (Johannes Gutenberg University Mainz, Germany), Dmitry V. Divine (Norwegian Polar Institute, Norway), Alessandro Incarbona (Palermo University, Italy), Darrell S. Kaufman (Northern Arizona University, USA), Nikita Kaushal (University of Oxford, United Kingdom), Robert M. Klaebe (The University of Adelaide, Australia), Hannah R. Kolus (Northern Arizona University, USA), Guillaume Leduc (Aix Marseille University, CNRS, IRD, INRAE, Coll France, CEREGE, France), Shreyas R. Managave (Indian Institute of Science Education and Research, India), P. Graham Mortyn (Universitat Autonoma de Barcelona (UAB), Spain), Andrew D. Moy (Australian Antarctic Division, Australia), Anais J. Orsi (Laboratoire des Sciences du Climat et de L'Environnement, France), Judson W. Partin (University of Texas at Austin, USA), Heidi A. Roop (University of Minnesota, USA), Marie-Alexandrine Sicre (CNRS and Sorbonne Université, France), Lucien von Gunten (PAGES International Project Office, Switzerland), and Kei Yoshimura (The University of Tokyo, Japan)

* Corresponding author. Email: bkonecky@wustl.edu
The response of the global water cycle to changes in global surface temperature remains an outstanding question in future climate projections and past climate reconstructions. The stable hydrogen and oxygen isotope compositions of precipitation ($\delta_{precip}$), meteoric water ($\delta_{MW}$), and seawater ($\delta_{SW}$) integrate processes from microphysical to global scales and thus are uniquely positioned to track global hydroclimate variations. Here, we evaluate global hydroclimate during the past 2,000 years using a globally distributed compilation of proxies for $\delta_{precip}$, $\delta_{MW}$, and $\delta_{SW}$. We show that global mean surface temperature exerted a coherent influence on global $\delta_{precip}$ and $\delta_{MW}$ throughout the past two millennia, driven by global ocean evaporation and condensation processes, with lower values during the Little Ice Age (1450-1850) and higher values after the onset of anthropogenic warming (~1850).

The Pacific Walker Circulation is a predominant source of regional variability, particularly since 1850. Our results demonstrate rapid adjustments in global precipitation and atmospheric circulation patterns – within decades – as the planet warms and cools.

Recent global syntheses of paleoclimate “proxy” data have constrained global mean surface temperature (GMST) changes during the past 2000 years (i.e., the Common Era, CE), providing critical context for anthropogenic warming¹,². Yet despite the importance of water resources to society, contemporaneous variations in the global water cycle — including precipitation, evapotranspiration, atmospheric circulation, and modes of climate variability that affect these processes — remain underconstrained³,⁴.

Stable hydrogen and oxygen isotope ratios ($\delta^{2}\text{H}$ and $\delta^{18}\text{O}$) in environmental waters are well positioned to provide a global picture of hydroclimate. Evaporation, condensation, freezing, and
other phase changes in the water cycle differentially impact (fractionate) heavy versus light isotopes, causing the δ²H and δ¹⁸O of precipitation (δprecip), precipitation-derived meteoric waters such as lake and soil water (δMW), and seawater (δSW) to integrate and record hydrological processes on timescales from minutes to millions of years⁵⁻¹². Variations in δprecip, δMW, and δSW are subsequently incorporated into diverse geologic materials including speleothem and coral carbonate, glacial ice, and tree cellulose. By synthesizing such data from a variety of sources, it is therefore possible to infer changes in a powerful suite of hydroclimatic variables: δprecip, which reflects atmospheric factors such as condensation temperature, precipitation amount, rainout history, and moisture source⁵; δMW from lake, soil, and groundwaters, which reflects variations in δprecip and surface water evaporation; and δSW, which reflects δprecip, seawater evaporation, and mixing¹³,¹⁴.

We analyzed proxies for δprecip, δMW, and δSW from the recently published PAGES Iso2k database, which contains 759 globally-distributed paleoclimate records from coral, tree, ice, speleothem, lake, and marine sites¹⁵. The database includes extensive metadata designed to facilitate cross-archive comparison, including interpretations from the original publications and supplemented with information from a team of over 50 archive experts (see ref. ¹⁵ for details on database design). Metadata fields include original climatic interpretations, proxy system transformations, and the multiple environmental drivers of the isotopic composition of the measured material in each record (hereafter “isotope interpretation”). Each record is further classified into one of three primary isotope interpretation groups: 1) δprecip; 2) Effective Moisture (EM), i.e., the balance between precipitation and evaporation, with higher EM reflecting higher precipitation relative to evaporation; or 3) the *in situ* temperature of the environmental medium
during the formation of the proxy sensor or archive¹⁵ (we note most records in the temperature category are marine carbonates whose δ¹⁸O primarily reflects seawater temperature, with only minor influence of δsw¹⁶). Regional and global analyses were performed on these three isotope interpretation groups to distinguish patterns in different reservoirs of the water cycle (precipitation, surface water, seawater), without a priori assumptions about the climatic drivers of each record's variability (e.g., upstream monsoon intensity, regional air temperature), which are more subject to change as records are re-interpreted over time (Methods).

For each group we created composite records of global δ¹⁸O anomalies relative to the 0–2000 mean (Δ¹⁸O), including δ²H records scaled to δ¹⁸O-equivalent variance, using a dynamic compositing method that was previously employed to reconstruct paleotemperature in a manner that robustly handles proxy timeseries of different lengths, resolutions, and coverage periods (ref. ¹⁷ and Methods). In addition to calculating composites, we performed Principal Component Analysis (PCA) on a subset of higher-resolution (i.e., ≤30 year bins) individual proxy records with >85% temporal data coverage during the Last Millennium (LM; 850–1850) to reveal dominant spatiotemporal modes of variability (Methods). We compared these results to an ensemble of three full-forcing LM experiments with the water isotope-enabled Community Earth System Model (iCESM)¹⁸⁻²⁰.

The three Δ¹⁸O composites (Fig. 1) display similar patterns, with notable differences in the magnitude of centennial-scale variability. During the first millennium, composite Δ¹⁸O of δprecip-driven records (hereafter, composite Δ¹⁸Oδprecip) was relatively stable, whereas composite Δ¹⁸O of EM- and temperature-driven records (composite Δ¹⁸OEM and Δ¹⁸OTemp) increased and
decreased, respectively. During the LM, all three composites show a monotonic trend from ~800 to ~1700, and a reversal of that trend since the 19th century. These patterns broadly echo the temporal evolution of GMST during the LM¹, but with different magnitudes depending on the primary environmental interpretation of δ¹⁸O. From 1000 to 1850, a global cooling of 0.25 (±0.1) °C ¹ corresponds to a change of -0.27 (± 0.00021), +0.02 (± 0.00019), and +0.09‰ (± 0.0002) in composite Δ¹⁸Oδprecip, Δ¹⁸OEM, and Δ¹⁸Otemp, respectively (Fig. 1A-C). Subsequent warming of 0.65 (± 0.06) °C from 1850 to 2000 corresponds to a change of +0.56 (± 0.00012), +0.62 (± 0.00017), and -0.16‰ (± 0.00023), respectively. Uncertainties in these estimates are based on differences between the respective time periods across the full ensemble (Methods).

Composite Δ¹⁸OEM displays pronounced centennial-scale variability, with distinct positive excursions from 300–500, 700–900, and 1800–2000. Variability in composite Δ¹⁸OEM is at least twice the magnitude of composite Δ¹⁸Otemp or Δ¹⁸Oδprecip (Fig. 1), likely due to the strong influence of surface-water evaporation on lake and seawater δ¹⁸O, which amplifies the δprecip signal relative to noise⁵,⁷.

The first principal component (PC1) of each category of records is dominated by a monotonic trend over the LM, with smaller centennial-scale fluctuations (Fig. 2A). Similar to the composites, the gross trend in each PC1 corresponds to a decrease in GMST. Site loadings on each PC1, however, differ by region and by the primary environmental driver of δ¹⁸O (Fig. 2B-D). For example, positive trends and negative PC1 loadings are evident in δ¹⁸Otemp at almost all extra-tropical locations (Fig. 2D), consistent with the impact of ocean cooling on the δ¹⁸O of
167 marine carbonates. Contrastingly, trends are insignificant in the Indo-Pacific Warm Pool where
168 $\delta^{18}$O$_{sw}$ influence likely confounds the temperature signal$^{21}$.

169
170 **Influence of temperature on global $\delta_{MW}$**
171 Together, the $\Delta^{18}$O composites and PC1 suggest that GMST exerts a first-order control on
172 temporal changes in global $\delta_{MW}$ during the CE. The relationship of composite $\Delta^{18}$O$_{\delta precip}$ with
173 temperature is 0.68 ± 0.20 ‰/°C for the full CE and 0.78 ± 0.19 ‰/°C from 850-2000, based on
174 regression of the 30-year binned values (Methods) (Fig. 1A). A positive relationship between
175 GMST and $\delta_{precip}$ may be expected from high latitude ice cores$^{22-24}$, but positive relationships in
176 composite $\Delta^{18}$O$_{\delta precip}$ persist even when such records are excluded (Extended Data Fig. 1), and
177 also occur at mid and low latitudes, especially after 1850 (Extended Data Fig. 2). In addition to
178 the composites, the positive relationship between $\delta_{precip}$ and GMST is evident in PC1 of the $\delta_{precip}$
179 records spanning the entire CE (Extended Data Fig. 3A-B) and in iCESM simulations (Fig.
180 3A,D; Extended Data Fig. 4). In iCESM, the regression slope between GMST and global, 30-
181 year smoothed, mean annual $\delta^{18}$O$_{precip}$ from 850-2000 is 0.25‰/°C, and 0.48‰/°C when
182 calculated using Iso2k site locations (which, for $\delta_{precip}$ records, excludes nearly all ocean grid
183 cells and regions where site-level regression slopes are close to 0; Extended Data Fig. 4). The
184 iCESM slopes are shallower than in the Iso2k data, but the global slope is consistent with earlier
185 GCM estimates of 0.3‰/°C during the last deglaciation$^{25}$ and simplified model estimates of
186 zonal-mean temporal slopes of ~0.1-0.4‰/°C for most latitudes$^{26}$. The discrepancy between the
187 Iso2k and iCESM temporal slopes may reflect model biases in extratropical moisture transport
188 and positive $\delta^{18}$O$_{precip}$ biases at the high latitudes$^{18}$. However, observationally, the global
189 temporal slope between GMST and $\delta^{18}$O$_{precip}$ has not been quantified observationally– unlike
spatially-derived slopes between site-level air temperature and δ¹⁸Oₚrecip²³,²⁷,²⁸— because modern δ¹⁸Oₚrecip measurements in most regions are either absent or too short and discontinuous (most records <10 years) to do so beyond the scale of a few years²⁹,³⁰. Despite differences in magnitude, both iCESM and Iso2k data demonstrate that a positive relationship between GMST and δMW is a persistent feature of the global water cycle for the CE. The 30-year-binned Iso2k temporal slope of 0.68 ± 0.20 ‰/°C should be considered a benchmark to be tested as longer observations and reanalyses become available.

Ice core and marine sediment studies have long recognized the importance of air temperature in driving high latitude δMW, and of global ice volume in driving global δSW, on glacial-interglacial timescales and across the Cenozoic²⁴,²⁶. However, the nature of this relationship has remained uncertain on the shorter (decadal to centennial) timescales most relevant for understanding modern climate change and its impacts on water resources. Our results provide the first observational evidence that GMST drives temporal changes in δMW and δSW, and therefore changes in the hydrological cycle, on such timescales. Mechanisms other than global ice volume, which has not changed substantially during the CE, are therefore required to explain this relationship.

Although the Iso2k analyses are the first to document it, stable isotope theory and experimental studies provide ample foundation for a positive imprint of GMST on global δMW, which is what our spatially-distributed (albeit not spatially continuous) proxy network approximates. At the global scale, a relationship between δMW and GMST will integrate all processes that relate local-scale δMW with GMST, while balancing out regional distributions of heavy vs. light
213 isotopologues throughout the water cycle. Variations in global δ<sub>MW</sub> therefore reflect variations in
214 the isotopic composition of the global oceanic and atmospheric reservoirs. Atmospheric δ<sub>vapor</sub> is
215 ultimately governed by the isotopic composition of water evaporated from the oceans.

216 Equilibrium fractionation is greater at lower temperatures<sup>27</sup>, so in a cooler world, the liquid-to-
217 vapor difference in δ<sup>18</sup>O is higher than in warmer conditions (i.e., lower δ<sub>vapor</sub> in the saturated
218 atmospheric layer above the ocean surface). Kinetic fractionation further decreases atmospheric
219 δ<sub>vapor</sub> as newly evaporated vapor diffuses and mixes into the undersaturated free atmosphere,
220 with stronger fractionation when the lower troposphere is less humid<sup>28,29</sup>. A globally cooler and
221 drier troposphere should therefore decrease δ<sub>vapor</sub> to a greater extent than a warmer and more
222 humid troposphere. Global δ<sub>vapor</sub> is further modified by temperature-dependent fractionation as
223 precipitation forms and condensation preferentially removes <sup>18</sup>O and decreases δ<sub>vapor</sub>, with
224 greater isotopic discrimination at colder condensation temperatures<sup>27,30</sup>.

225
226 iCESM simulations indicate multiple mechanisms play a role (Fig. 3). Higher GMST is
227 associated with higher δ<sup>18</sup>O of water vapor flux from evaporating seawater into the saturated
228 boundary layer (0.14‰/°C), following the Craig-Gordon model<sup>13,31</sup>, slightly higher relative
229 humidity with respect to sea surface temperature (0.77‰/°C) (Methods), and ultimately, higher
230 atmospheric δ<sup>18</sup>O<sub>vapor</sub> (0.30‰/°C). Global δ<sub>precip</sub> therefore reflects not only the δ<sub>vapor</sub> determined
231 during ocean evaporation and mixing into the free atmosphere, but also the subsequent depletion
232 of <sup>18</sup>O in atmospheric vapor in the atmosphere through condensation, which via equilibrium
233 fractionation is stronger in a cooler climate. Stronger depletion of heavy isotopes in precipitation
234 also occurs with more vigorous circulation and shorter residence time of atmospheric moisture<sup>32</sup>.
235 For precipitation integrated over the timescales of most proxy systems (weeks to months), δ<sub>precip</sub>
236 could become even lower when global precipitation rates are high relative to the amount of
237 precipitable water in the atmosphere, and higher when atmospheric humidity increases to a
238 greater extent than precipitation rate, as in the 20th century (Fig. 3; ref. 33).

239
240 Theory and simplified model experiments suggest that a cooler global atmosphere could also
241 favor stronger latitudinal gradients in δprecip<sup>5</sup>, as lower temperatures increase distillation along the
242 water's path from the subtropics to the poles<sup>23,34</sup>. A cooler atmosphere may also shift evaporative
243 source regions equatorward, driving a greater fraction of high latitude precipitation to be sourced
244 remotely<sup>35</sup> which strongly affects δprecip at those latitudes<sup>36</sup>. The latitudinal gradient in Iso2k
245 δprecip records between 40° to 70° N and S is approximately -0.48‰ per degree latitude, in
246 agreement with observed modern slopes of ~ -0.3 to -0.6‰ per degree latitude<sup>23,34</sup>. However, we
247 observe no difference in this gradient between the Little Ice Age and the 20th century — the
248 globally coldest and warmest intervals of the CE, respectively<sup>1,37</sup> (Extended Data Fig. 5). Hence,
249 either equator-to-pole Rayleigh distillation is counterbalanced by changes in spatial patterns of
250 global evaporative recharge and moisture transport relative to precipitation<sup>35,36</sup>, or the spread in
251 Iso2k δprecip records is too large to resolve changing gradients during the LM when temperature
252 changes are relatively small. As spatial gradients are averaged out at the global scale, further
253 experiments with both simplified models and isotope-enabled GCMs are needed to quantitatively
254 decompose the relative importance of ocean evaporation, condensation, and global precipitation
255 and evaporation patterns on the global-scale GMST- δprecip relationship.

256
257 Patterns of regional variability
Despite its influence on global mean δ<sub>MW</sub>, GMST explains neither the spatial patterns nor the shorter-term variability in δ<sub>MW</sub> and δ<sub>SW</sub> (Fig. 2, Extended Data Fig. 6). Instead, site-level PC1 loadings, LM trends, and centennial anomalies in δ<sup>18</sup>O<sub>precip</sub> and δ<sup>18</sup>O<sub>EM</sub> are spatially heterogeneous (Figs. 1-2, Extended Data Fig. 7). Although non-climatic processes and noise likely contribute to this variability, this result also indicates that regional water balance, atmospheric circulation, and precipitation characteristics dominate regional δ<sub>MW</sub> not just in the late 20th/early 21st centuries<sup>38,39</sup> but also on multi-decadal to millennium-long timescales. In the Arctic and the tropical Andes, negative trends in δ<sup>18</sup>O<sub>precip</sub> and positive PC1 loadings agree with documented Little Ice Age climate changes: high-latitude cooling, which lowered δ<sub>precip</sub> in the Arctic<sup>40</sup>, and intensification of the South American summer monsoon, which lowered δ<sub>precip</sub> in the Andes<sup>41</sup>. In eastern China and the Maritime Continent, positive PC1 loadings and negative trends in δ<sup>18</sup>O<sub>precip</sub> and δ<sup>18</sup>O<sub>EM</sub> over the LM indicate increasing effective moisture, with depletion of heavy isotopes at some sites enhanced by intensifying convection and moisture transport within the East Asian and Australasian summer monsoons<sup>42,43</sup>. In northwestern North America and Central America, negative PC1 loadings and positive trends in δ<sup>18</sup>O<sub>precip</sub> and δ<sup>18</sup>O<sub>EM</sub> indicate increasing δ<sup>18</sup>O<sub>MW</sub> and declining EM over the LM, while southwestern North America shows the opposite pattern, consistent with other LM proxy reconstructions<sup>44,45</sup>. δ<sup>18</sup>O<sub>EM</sub> from lakes exhibit the clearest trends in western North America, likely because δ<sub>precip</sub> in the western North American interior is especially sensitive to changes in water balance<sup>39</sup> and enhanced evaporation during dry periods then amplifies increases in lake water δ<sup>18</sup>O. This explains why coherent patterns emerge in western North American isotopic proxy records despite substantial heterogeneity in compilations that blend isotopic and non-isotopic records<sup>45</sup>.
Opposing precipitation anomalies between southwestern North America and Central Asia versus the Maritime Continent and South Asia during the LM have previously been attributed to variability in the El Niño-Southern Oscillation<sup>46</sup>, the interannual component of the east-west atmospheric overturning circulation over the tropical Pacific known as the Pacific Walker Circulation (PWC)<sup>47</sup>. An underlying influence of the PWC may also explain some of the regional coherency in Iso2k δ<sub>MW</sub>. In modern precipitation, strengthening and weakening of the PWC drives opposing precipitation and δ<sup>18</sup>O<sub>precip</sub> anomalies in Asia and the Americas due to changes in precipitation amount, moisture source, and transport pathways<sup>38</sup>. For example, rerouting of the jet stream during El Niño events alters Pacific Ocean moisture trajectories toward southwestern rather than northwestern North America<sup>48,49</sup>. To the extent that these relationships persist beyond interannual timescales, changes in moisture transport and EM provide a plausible mechanism for multidecadal to centennial patterns in Asian and western North American δ<sup>18</sup>O<sub>EM</sub>.

Evidence for the PWC's influence on Iso2k records is stronger on sub-decadal timescales during the historical period (1850-2005), when instrumental observations are available for direct comparison. We correlated PC1 of 3-year binned Iso2k EM records with observed sea level pressure (SLP) and found that the associated pattern mimics the observed global expression of the PWC (Fig. 4a). Similarly, in iCESM experiments, PC1 of soil water δ<sup>18</sup>O (δ<sup>18</sup>O<sub>soil</sub>) at Iso2k EM sites displays an SLP pattern that resembles iCESM's PWC (Fig. 4b). There is vigorous debate surrounding the extent to which the PWC fluctuates on multidecadal and longer timescales, either due to internal variability in the climate system or to external forcing<sup>50-52</sup>. Our results suggest that the PWC was the predominant influence on interannual δ<sup>18</sup>O<sub>precip</sub> not only
from 1982 to 2015³⁸, but at least since 1850, and likely on multi-decadal timescales prior to
1850.

Our results reveal a remarkably consistent relationship between GMST and δ<sub>MW</sub> throughout the
CE, despite relatively constant ice-ocean boundary conditions. Between 1850-2000, global δ<sub>precip</sub>
increased by at least 0.56‰, indicative of a warmer, more humid troposphere. Global δ<sub>MW</sub> and
regional δ<sub>MW</sub> appear to adjust to changing temperature and atmospheric circulation patterns,
respectively, within decades – similar to the timescale of the forcing itself. Expanded δ<sub>MW</sub>
obervational networks are critical for detecting and attributing shifts in rainfall, drought, and
circulation as the planet continues to warm.

Acknowledgments

Iso2k is a contribution to Phase 3 and 4 of the PAGES 2k Network. PAGES received support
from the Swiss Academy of Sciences, the US National Science Foundation, and the Chinese
Academy of Sciences. We thank three anonymous reviewers for helpful comments which
improved this manuscript.

Support for this work includes NSF-AGS grant 1805141, NSF-AGS PRF 1433408, and a David
and Lucile Packard Foundation Fellowship in Science and Engineering to Bronwen L. Konecky.
Samantha L. Stevenson was supported by NSF-AGS grant 1805143. Georgina M. Falster was
supported by NSF-AGS grant 1805141 to B. Konecky and by the Australian Research Council
Centre for Excellence for Climate Extremes (CE170100023). Jessica Conroy was supported by
NSF-CAREER-1847791. Diane M. Thompson was supported by NSF-CAREER-1945479, NSF-
1931242, and NSF-2002460. Lukas Jonkers was funded through PalMod, the German paleoclimate modeling initiative. PalMod is part of the Research for Sustainable Development initiative funded by the German Federal Ministry of Education and Research (BMBF). Olga Churakova (Sidorova) was supported by RSF project 21-17-00006. Kristine DeLong was supported by South Central Climate Adaptation Science Center Cooperative Agreement G19AC00086. Belen Martrat (BM) was funded by the RYC-2013-14073 programme and the LINKA20102 and CEX2018-000794-S projects. Jonathan Tyler was supported by the Australian Research Council Discovery Project DP190102782. Nerilie J. Abram was supported by the Australian Research Council through a Future Fellowship (FT160100029), Special Research Initiative for the Australian Centre for Excellence in Antarctic Science (SR200100008) and the Centre of Excellence for Climate Extremes (CE170100023).

Author contributions

Analyses presented in the main text and extended data were conceived and performed by BLK, NPM, GMF, SS, MJF, AA, DMT, MDJ, JJT, EKT, JLC, SGD, and LJ.

Results were analyzed and interpreted by BLK, NPM, GMF, SS, MJF, AA, DMT, MDJ, KLD, JJT, BM, and EKT, with input from all authors.

The manuscript was written mainly by BLK, NPM, GMF, SS, MJF, AA, DMT, MDJ, KLD, JJT, BM, EKT, JLC, SGD, LJ, and HRS, with additional contributions from OVC, ZK, TO, TJP, and GS.
349

350 All Iso2k Project Members created the Iso2k database and edited the manuscript.

351

352 BLK directed the project, led the overall design of the study, and led the writing of the manuscript.

354

355 Competing interests:

356

357 The authors declare no competing interests.

358

359
360 **Figure legends/captions**

361 **Fig. 1.** Global composite $\delta^{18}$O anomalies. **a-c,** 30-year binned proxy $\Delta^{18}$O (black line, ensemble median; dark shading, 25$^{th}$–75$^{th}$ percentiles; light shading, 2.5$^{th}$–97.5$^{th}$ percentiles) for each group of Iso2k records$^{15}$, anomalies relative to 2000-year mean. Gray shading depicts ensemble 2.5$^{th}$–97.5$^{th}$ percentile GMST anomaly relative to 1961-1990 mean$^{1}$. Note y-axis in (**c**) is reversed to orient values upwards for warmer periods. **d-f,** Records contributing to each composite. Symbol shape, archive type; size, record duration; shading, correlation between that record and the corresponding composite; bold outline, $p<0.05$. Correlations are Pearson's r, two-sided, with no adjustment for multiple correlation. Maps created in R, using coastlines from Natural Earth.

369

370

371 **Fig. 2.** First principal components and trends in Iso2k records over the Last Millennium. **a,** PC1 of 30-year binned Iso2k records$^{15}$ (850 to 1840), by interpretation group (colored lines), compared to GMST anomaly (gray shading) as in Fig. 1. **B-d,** Spatial loadings on PC1 (symbol shading) for each group. Correlations are Pearson's r, two-sided, with no adjustment for multiple correlation. Inner symbols, archive type; outer symbol shape, slope of significant ($p<0.05$) linear trend in the $\delta^{18}$O of each individual record; labels, variance explained by PC1 and number of records used. Maps created in R, using coastlines from Natural Earth.

372

373

374

375

376

377

378

379 **Fig. 3.** Global hydroclimate and isotopic anomalies in iCESM$^{15}$. **A,** Ensemble mean, amount-weighted mean annual $\Delta^{18}$O$_{precip}$ anomaly relative to 1961 to 1990 mean (thin gray line) and 30-year running mean (thick black line) for all grid cells, and for only grid cells containing Iso2k site locations (blue lines). Note right Y axis in **a** is scaled 2x that of the left. **B-f,** as in **a,** but for

380

381

382
383 **(b)** global column-averaged $\Delta^{18}$O$_{vapor}$, **(c)** $\Delta^{18}$O of evaporative flux from the global oceans, **(d)** GMST, **(e)** relative humidity over the low- and mid-latitude oceans (60°N-60°S), with respect to saturation vapor pressure at sea surface temperature$^{13,28,29}$, and **(f)** global precipitation rate.

386
387 **Fig. 4. Isotopic fingerprint of the Pacific Walker Circulation in Iso2k¹⁵ and iCESM. a,**
388 Shading is the correlation between observed SLP and PC1 of 3-year binned Iso2k EM records
389 (1850 to 2005). Symbol colors denote site loading on PC1. Unfilled contours denote correlation
390 between SLP and PWC index $\Delta$SLp. **b,** Correlation between SLP and ensemble-mean PC1 of 3-
391 year binned 0-10 cm $\delta^{18}$O$_{soil}$ in iCESM experiments (1850 to 2005), using grid cells of Iso2k EM
392 sites (symbols). **c,** As in **b**, but for 30-year binned data from 850-1850. **d-f,** Time series of PC1
393 (black line) and $\Delta$SLP (gray line) corresponding to panels **a-c**. Maps created in MATLAB using
394 m_map for coastlines.

395
396 **Methods**
397 **Disaggregation of analyses by isotope interpretation**
398 All global analyses presented in this paper are disaggregated by isotope interpretation rather than
399 by the original authors' climatic interpretation. This approach avoids building into our synthesis
400 *a priori* assumptions about the specific climatic variable driving $\delta^{18}$O and $\delta^{2}$H variability in each
401 record, the basis for which is not always explained in the original publications. In addition, many
402 records contain multiple climatic interpretations (e.g., speleothem $\delta^{18}$O being driven mainly by
403 monsoonal rainfall amount, but amplified by the ratio of summer vs. winter moisture source
404 changes$^{1-3}$), and the relative importance of each variable was impossible to objectively evaluate
405 for every time scale, region, and proxy system in the database. Finally, climatic interpretations of
406 isotope-based proxy records are continually evolving as new information emerges from
407 environmental monitoring and modeling studies (e.g., moisture source and transport
408 characteristics driving Greenland ice core δ¹⁸O variability⁴⁻⁸; seawater δ¹⁸O driving coral δ¹⁸O in
409 some regions⁹). The isotope interpretation groupings are less subject to interpretation because
410 they chiefly represent isotope systematics and the physical pools of environmental waters that
411 each proxy sensor imbibes.

412
413 We separated δprecip and EM records because EM proxies integrate information about evaporation
414 that is not expected from pure δprecip proxies, which are mainly driven by factors influencing
415 condensation (e.g. air temperature, degree of rainout, import of moisture from different
416 sources)¹⁰,¹¹.

417
418
419 **Calculation of Δ¹⁸O composites**

420 All records were extracted from the PAGES Iso2k Database v 1.0.0¹⁵,⁵³. The database contains
421 metadata on the principal determinants of isotopic composition in the measured material (e.g.,
422 the δ¹⁸O of precipitation that forms glacial ice or cave dripwaters) and the record's climatic
423 interpretation (e.g., atmospheric temperature at condensation level, rainout due to monsoon
424 intensity), as interpreted by the original studies' authors and our team of archive experts¹⁵. For
425 each isotope interpretation group, we calculated an ensemble of 100 composite δ¹⁸O time series
426 for the CE. Prior to calculation, we filtered the database to only include the 'primary' time series
427 for each site, and then grouped records according to the primary driver of isotopic variability i.e.,
EM, temperature, or $\delta_{\text{precip}}^{15,53}$) (entitled 'Effective Moisture', 'Temperature', and 'P_isotope' in the 'isotope Interpretation 1 variable Group' metadata field of the Iso2k database).

For records with isotope interpretation 'P_isotope', we also calculated separate composites for a) glacier ice only, and b) excluding glacier ice (Extended Data Fig. 1). This separation was performed to assess whether composite $\Delta^{18}\text{O}_{\delta_{\text{precip}}}$ is overprinted by the large number of ice core records from high latitudes and high elevations, where temperature-driven isotopic fractionation may disproportionately affect $\delta_{\text{MW}}^{12}$. An especially strong temperature-$\delta_{\text{MW}}$ relationship may be unsurprising in the glacier ice $\delta^{18}\text{O}$ system. Glacier ice reflects $\delta_{\text{precip}}$ more directly than other proxy sensors, which reflect pools of meteoric water (e.g. soil water and lake water) that are influenced by precipitation and other secondary processes such as evaporation or aquifer mixing. Glacier ice is found at high latitudes and altitudes, where cold temperatures drive stronger stable isotope fractionation due to Rayleigh distillation and global patterns of precipitation vs. evaporative recharge$^{5-7,13-15}$. Most of the glacier ice records included in the Iso2k database are consistently interpreted as temperature indicators, and were included in the PAGES 2k$^{12,16}$ temperature database used for GMST calculations.

Despite these considerations, the observed global $\delta_{\text{MW}}$-temperature relationship persists even when glacier ice records are removed from the $\delta_{\text{precip}}$ composite (Extended Data Fig. 1), and in low- and mid-latitudes where local temperature effects on $\delta_{\text{precip}}$ are small$^{13}$ (Fig. 2). The overall patterns in PC1 are also similar regardless of whether glacier ice records are included or excluded (Extended Data Figs. 3, 8). Therefore, the inclusion or exclusion of glacier ice records does not substantially affect the composites or PCA, supporting that the strong temperature-$\delta_{\text{MW}}$
relationship in our data is due to the overall influence of temperature on the global pool of
meteoric water and not due to the strong effect of air temperature on high-latitude ice cores.

Records with ten or fewer data values within the CE were excluded from the analysis. δ²H records (n=45 in composite Δ¹⁸Oδprecip, n=12 in composite Δ¹⁸OEM) were divided by eight to scale the magnitude of their variance with that of δ¹⁸O in global meteoric waters⁵⁴. This was done to avoid erroneously high apparent climate variability at δ²H sites simply from the eight-times-higher variability in δ²H relative to δ¹⁸O, which arises from relative differences in equilibrium fractionation factors between liquid and vapor for these isotope ratios⁵⁴. Local slopes were not available for the 12 records from evaporative water bodies in composite Δ¹⁸OEM, so the global scaling of 8 was used.

Records in the Iso2k database have a wide variety of temporal resolution, length, and coverage over the CE¹⁵. To align records to a common interval and resolution, the data were averaged into equal bins of 30 years spanning the CE, which approximates the average resolution of the lower resolution archive types in the database (marine and lake sediments). Records contributing to each bin were mean-centered but not scaled by variance, as described below. To minimize aliasing, the data were binned following a modified nearest-neighbor annual interpolation procedure¹⁷. This approach accounts for the fact that for the non-annually-resolved data, the age of a sample at a given depth typically represents more than one year of accumulation (and up to several years or decades, or even longer for some low-resolution sedimentary records), and may therefore contain climate information that is relevant to more than one bin. The duration of each sample is not consistently recorded in the Iso2k database (or the primary references), so to
estimate sample coverage we calculated the distribution of gaps between adjacent observations in each timeseries, and used nearest-neighbor interpolation to estimate sample values spanning the intervals that are less than the 75th quantile of the distribution of all gaps between adjacent observations (consistent with Kaufman et al.¹⁷). Consequently, samples with resolution <30 years can potentially contribute to the weighted δ¹⁸O calculation for up to two bins, though with less weight given to samples further from their published age. In the case of records with resolution >30 years, this data-spreading step allows observations to impact the composites across multiple bins, consistent with their interpretation as multidecadal averages.

Record lengths also vary widely in the database, so there is no universal time period of common overlap within the CE. We therefore aligned the records using the Dynamic Compositing approach¹⁷, which uses randomly selected portions of each time series to adjust the δ¹⁸O variance, then iteratively adjusts the mean δ¹⁸O of each time series so that the mean of each record is minimally offset from all other records in the composite. This process was repeated for each of 100 ensemble members. Because the mean value of the records is adjusted during compositing, the composite values are now in relative Δ¹⁸O (in ‰), rather than in their original δ¹⁸O (or pre-scaled δ²H) values on the VSMOW-SLAP scale (i.e., in 1000·δ notation where 0‰ refers to standard mean ocean water). Therefore, for convenience, composite Δ¹⁸O was slightly shifted such that the mean of the ensemble median is 0‰. The final composite Δ¹⁸O values (Fig. 1) are therefore in units of ‰ anomalies relative to the 2000-year mean. Dynamic compositing was used in ref. ¹⁷ to produce calibrated reconstructions of paleotemperature. Unlike surface air temperature observations, long and complete timeseries of globally distributed δprecip, δMW, and δSW observations are insufficient (e.g. for δprecip, large spatial coverage gaps and few records
497 longer than 10 years³⁸), and so calibrated reconstructions are not possible until these
498 observational networks are improved.

499
500
501
502 PCA of Iso2k data
503 We used PCA to calculate Empirical Orthogonal Functions (EOFs) from a subset of records in
504 the Iso2k database that met requirements for temporal coverage (described below; Fig. 2). Prior
505 to calculating the EOFs, we filtered the database as described for the Δ¹⁸O composites, and
506 binned the raw data. PCA was performed using records in each of the three “isotope
507 interpretation” categories; for ‘P_isotope’ records, we performed separate PCA for a) glacier
508 ice, and b) ‘not glacier ice’ (i.e., all P_isotope records that are not from glacier ice). PCA was
509 performed on three temporal subsets: 0-1980 CE (30-year bins), 850-1840 (30-year bins), and
510 1850-2005 (3-year bins). Here we describe the methodology for the 850-1840 PCA that is shown
511 in the main text, then outline any adaptations that were made for the other two intervals.

512
513 Prior to PCA, all records were truncated to a fixed time interval of 850-1840. Other time
514 windows were explored and are presented in Extended Data Figs. 3 and 8. The 850-1840 interval
515 was chosen in order to focus on pre-industrial variability, to avoid issues with records or
516 chronologies not extending fully to either the present day or the start of the CE, and to provide
517 comparison with LM model experiments from 850-1850. (Note a 990-year interval was used
518 because it is divisible by the 30-year bin size). Results of PCA on the 0-1980 interval (with 30-
519 year bins) and 1850-2005 (with 3-year bins) are shown in Extended Data Figs. 3 and 8. All data
520 processing for these intervals was as below.

521
522 Records included in the PCA had to meet stricter time coverage and completeness requirements
523 than records included in the composites. Binned records with <85% temporal coverage during
524 the 850-1840 interval were excluded. EOF analysis was carried out on the matrix of remaining,
525 screened time series using the Data Interpolation Empirical Orthogonal Function (DINEOF)
526 method, which accounts for records with a small number of missing data points⁵⁵,⁵⁶. DINEOF
527 performs best when missing values are scattered randomly throughout the input data time series.
528 In the case of Iso2k records, missing data are generally concentrated at the ends of records, and
529 hence records with a large proportion of missing data have a large impact on the PC loadings
530 (though less impact on the overall PC time series). To check for obvious artifacts induced by the
531 data processing, we visually compared the binned and interpolated records with their raw data
532 equivalents. We performed the PCA on the binned and interpolated data matrix using the rda()
533 function of the 'vegan' package in R⁵⁷. We scaled the interpolated records to unit variance
534 because they are from different proxy systems and hence the raw data have widely ranging
535 variance.

536
537 To test whether variability explained by the principal components could be explained
538 stochastically, we tested the magnitude of the eigenvalues against a stochastic null hypothesis,
539 using a block bootstrap method that accounts for decadal persistence (Extended Data Fig. 9). For
540 the 850-1840 analysis, we split all raw data records that contributed to the EOF into ten-year
541 blocks; we chose ten years because this separates interannual variance from interdecadal
persistence. For each time series, we randomly selected an initial block, identified all blocks with similar means, and randomly sampled one block from these. The next selected block was the successor to that randomly sampled block. Blocks with similar means were identified using nearest neighbor-based ranking. For blocks with no successor (i.e., where there is a gap in the time series), the successor block was randomly selected from all other blocks. The process was then repeated 1000 times to produce an ensemble of time series. This method of re-sampling preserves the effects of 1) irregular time spacing, 2) the autocorrelation function (to lags >30 years), and 3) the time series processing steps (e.g., binning) on the correlation structure of the raw (unprocessed) time series. We followed exactly the same steps for the 0-1980 PCA, using a 30-year bin size.

The 1850-2004 PCA was performed using three-year bins, on records with 80% coverage. Block bootstrap was performed using 3-year blocks because 10-year blocks were too coarse for this short time period. As with the above, records were binned and standardized but not detrended.

## PCA with iCESM experiments

To compare the iCESM with Iso2k results, we performed PCA of 30-year binned iCESM surface soil water $\delta^{18}$O ($\delta^{18}$O$_{soil}$) i.e. the model variable most comparable to evaporation-sensitive systems such as lakes. We used data from Iso2k EM site locations, spanning 850-1850, and then regressed the ensemble-mean first PC against model SLP. The binned $\delta^{18}$O$_{soil}$ was calculated for gridpoint time series corresponding to the “effective moisture” sites, taken from each isotope-enabled Last Millennium Ensemble full-forcing ensemble member from the upper 0-10 cm of the soil profile in iCLM, then averaged to produce an ensemble mean. The upper 10 cm were chosen
because in the model, isotopic fractionation is clearest at this level whereas δ¹⁸O₃₀₈ of deeper soil layers rapidly approaches δ¹⁸Oₚₐₖₐ¹⁸. EOF analysis was performed using δ¹⁸O₃₀₈ because the “lake” land cover type in iCLM4 is non-fractionating¹⁸,⁵⁸ and hence no model variable is available for direct comparison with the large number of lake-based Iso2k δ¹⁸Oₚₐₖₐ proxies. Yet to first order, δ¹⁸O₃₀₈ is similarly controlled by both δ¹⁸Oₚₐₖₐ and surface evaporation, and δ¹⁸O₃₀₈ is also an adequate comparison to tree cellulose and speleothem records in the Iso2k EM category.

572

## 573 Correlations with sea level pressure and calculation of PWC index

Observed SLP and ΔSLP for the correlations against Iso2k HP PC1 (Fig. 4) were taken from the HadSLP2r dataset⁵⁹. Following the treatment of the Iso2k data, the HadSLP dataset was binned to 3 years but not detrended or deseasonalized (Fig. 4).

577

The index for the trans-Pacific SLP gradient (ΔSLP) is defined as anomalies (from the monthly climatology) in the difference between area-mean SLP over the central-eastern Pacific Ocean (160°W-180°W, 5°S-5°N) and the western Pacific and eastern Indian Oceans (80°E-160°E, 5°S-5°N)⁵². Positive ΔSLP values represent an increased zonal pressure gradient, and hence stronger PWC (and vice versa). ΔSLP was calculated using HadSLP2r for comparisons with Iso2k (Fig. 4).

584

## 585 Magnitudes of change in composite Δ¹⁸O and GMST

To estimate the magnitude of change in the three composite Δ¹⁸O time series, we subtracted the composite Δ¹⁸O value at 1000 CE from the composite Δ¹⁸O value at 1850 CE, for all 100

587
composite ensemble members. We likewise calculated the difference between 2000 CE and 1850 CE. Similarly for estimating the change in GMST during this interval, we subtracted the temperature anomaly at 1000 (1850) from the same at 1850 (2000), for all 7000 ensemble members. In all cases, we report the mean and standard deviation of the distributions of magnitudes of change.

## Calculation of isotope-temperature relationships in Iso2k data and iCESM experiments

We calculated the relationship between composite $\Delta^{18}O_{\delta_{\text{precip}}}$ (i.e., an approximation of anomalies in the global mean $\delta^{18}O$ of precipitation, given the spatial distribution of 305 $\delta_{\text{precip}}$-sensitive proxy records) and GMST between 850-2000. We binned all ensemble members from the most recent reconstruction of CE GMST1 to match the $\Delta^{18}O_{\delta_{\text{precip}}}$ composite. We then calculated linear regressions for 10000 unique combinations of composite $\Delta^{18}O_{\delta_{\text{precip}}}$ ensemble member (n = 100) on GMST ensemble member (n = 7000), for the two time intervals. We report the mean and standard deviation of the distribution of regression coefficients, for the two time intervals.

For iCESM, isotope/temperature relationships were calculated for 850-2000 (the period of overlap with Iso2k composites) using area-weighted globally averaged, mean annual surface temperature and area-weighted, amount-weighted, mean annual global average $\delta^{18}O_{\text{precip}}$ from the ensemble-mean of three isotope-enabled Last Millennium Ensemble full-forcing simulations. Before determining the regression slopes, we calculated 30-year running means for both surface temperature and $\delta^{18}O_{\text{precip}}$ (Fig. 3). Regression slopes in the main text are reported for GMST vs. two time series: global mean $\delta^{18}O_{\text{precip}}$ (using all grid
611 cells), and mean $\delta^{18}$O$_{precip}$ calculated only for grid cells containing locations of Iso2k
612 primary time series data contributing to the composites between 850-2000 CE.
613
614 Similarly, to estimate the amount of variance in individual Iso2k primary time series explained
615 by proxy estimates of global temperature (Extended Data Fig. 6), we calculated R$^2$ for the 30-
616 year binned data across the interval 1-2000, ignoring bins that contained no observations.
617 Correlations were only calculated if at least six bins overlapped the bins from the GMST
618 reconstruction.

619

620 Calculation of Relative Humidity normalized to SST

621 Relative humidity was calculated from iCESM experiments for all grid points over the oceans
622 from 60°N-60°S, i.e. the portion of the lower troposphere receiving the majority of evaporated
623 water from the surface oceans. RH$_{SST}$ in Fig. 3 was calculated as the relative humidity of the
624 surface-most model layer, normalized to the saturation vapor pressure at the temperature of the
625 surface ocean rather than the air, following the physical principles of the Craig-Gordon model
626 for an evaporating water body$^{13,28,29}$. All global-mean time series in Fig. 3 are area-weighted.

627

628 Calculation of trends

629 Trends (Fig. 2) were calculated as the slope of the linear regression from 850-1840 for a subset
630 of Iso2k records meeting the following criteria: 1) designated ‘primary time series’, 2)
631 containing at least one data point in the first and last 50 years of the time interval; and 3)
632 containing at least 20 data points over the full time interval (i.e., $\ge$50 year average resolution).
For iCESM data, trends were calculated as the slope of the linear regression for monthly data covering the 850-1850 and 1850-2005 time intervals.

## 200-year standardized anomalies

We created standardized anomaly ('z score') maps to aid interpretation of temporal variability in the composite time series (Extended Data Fig. 7). Similarly to the analyses described above, we filtered the database to only include the 'primary' isotope ($\delta^{18}$O or $\delta^{2}$H) time series for each site, and then grouped records according to the primary driver of isotopic variability. We then filtered this subset of data sets to only include records spanning >600 years within the CE. We averaged those records into 200-year bins, and then calculated z-scores for each bin by subtracting the mean of all data points within the CE from the bin average, and then dividing by the standard deviation of all data points within the CE. We performed this analysis using both 'odd' (100, 300, 500 etc.) and 'even' (200, 400, 600) centuries as bin centers and showed only the time periods relevant to the main text in Extended Data Fig. 7.
647

648 **Data Availability Statement**

649 The Iso2k Database⁵³ is available for download at https://doi.org/10.25921/57j8-vs18 and is

650 accessible via the NOAA/WDS Paleo Data landing page at

651 https://www.ncdc.noaa.gov/paleo/study/29593. Composites and principal components datasets

652 generated for this manuscript are available through NOAA/NCEI.

653

654 **Code Availability Statement**

655 Codes to reproduce the results from this manuscript are available via Github

656 https://github.com/nickmckay/iso2kNatureGeoscience2023 and archived on Zenodo.

657

658
EXTENDED DATA FIGURES 1-9

<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

Extended Data Figure 1. Composite $\Delta^{18}\text{O}_{\text{δprecip}}$ calculated with and without glacier ice records. As per Fig. 1, the black line with coloured shading shows the 30-year binned proxy $\delta^{18}\text{O}$ anomaly from Iso2k records¹⁵ (black line, ensemble median; dark shading, first and third quantiles; light shading, 2.5th/97.5th percentiles). (a) Composite of $\Delta^{18}\text{O}_{\text{δprecip}}$ records from all archives other than glacier ice. (b) Composite of only glacier ice records. Records contributing to each bin are mean-centered but not scaled according to that record's variance (see Methods). Anomalies are in ‰ relative to the 2000-year mean. Gray shading depicts the ensemble 2.5 and 97.5 percentile of the 31-year Butterworth-filtered Global Mean Surface Temperature (GMST) anomaly relative to the 1961-1990 mean¹.

659

660

661

662

663

664

665

666

667

668

669

670

671

672

673
<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

674

675 Extended Data Figure 2. Composite Δ¹⁸Oδprecip calculated using Iso2k records¹⁵ falling within 30-degree latitudinal bins. a, 60-90°N (n=76); b, 30-60°N (n=86); c, 0-30°N (n=27); d, 0-30°S (n=39); e, 60-90°S (n=77). As per Fig. 1, black line with coloured shading shows the 30-year binned proxy δ¹⁸O anomaly (black line, ensemble median; dark shading, first and third quantiles; light shading, 2.5th/97.5th percentiles) and gray shading depicts the ensemble 2.5 and 97.5 percentile of the 31-year Butterworth-filtered GMST anomaly relative to the 1961-1990 mean. Black text denotes mean regression slope (± 1 standard deviation) of regional composite

676

677

678

679

680

681
682 $\Delta^{18}$O vs. GMST from 850-2000. Regional composite for 30-60°S not calculated due to
683 insufficient number of records from those latitudes (n=2).

684

<figure><img src="image_9.png" /></figure>

685

686 **Extended Data Figure 3.** First Principal Component (PC1) of Iso2k records¹⁵ during the
687 **full Common Era.** Symbols, lines, and shading are as in Fig. 2 but for the interval 0-1980 (30-
688 year bins), and without trends depicted (i.e. constant shape for outer symbols). $\delta_{\text{precip}}$ PC1
689 explains 19% of the total variance (n = 44). Effective Moisture PC1 explains 19% of the total
690 variance (n = 19). Temperature PC1 explains 25% of the total variance (n = 19). Maps created in
691 R using coastlines from Natural Earth.

692
<figure><img src="image_10.png" /><figcaption>Extended Data Figure 4. iCESM relationship between GMST and grid cell-level δ<sup>18</sup>O<sub>precip</sub>.</figcaption></figure>

Shading depicts the regression coefficient in ‰/°C between mean annual, amount-weighted, 30-year running mean δ<sup>18</sup>O<sub>precip</sub> at every grid cell vs. 30-year running mean, area-weighted GMST for the mean of 3 full-forcing, isotope-enabled Last Millennium Ensemble members. Map created in MATLAB using m_map for coastlines.
<figure><img src="image_11.png" /><figcaption>Extended Data Figure 5. δ¹⁸O vs. latitude in the Iso2k database¹⁵ for the modern era (a-d; 1950-2018) and the Little Ice Age (e-h; 1450-1850). Data are plotted for the three main isotope interpretation categories (a-c, e-g) and for δ¹⁸O of glacier and ground ice (d, h), the archives that preserve precipitation most directly. Black lines indicate polynomial fits to the 0.1, 0.5, and 0.9 quantiles. Black squares on panels (d) and (h) indicate 40 and 70 degrees latitude, the interval for which the mean of the gradient function was calculated (-0.48‰/° latitude for both time periods; see main text). All records are plotted in permil on the VSMOW-SLAP scale for ice core records; PDB or VPDB for all other archives).</figcaption></figure>
715

716

717

718
<figure><img src="image_12.png" /><figcaption>Extended Data Figure 6. Relation between Iso2k¹⁵ primary timeseries and Global Mean Surface Temperature (GMST). (a) Fraction of variance explained (R²) in δₚrecip primary time series by changes in PAGES 2k global mean surface temperature (GMST) over the past 2000 years. Both Iso2k time series and GMST were averaged into 30 year bins before calculating correlations (See Methods). (b) As in (a), but for Effective Moisture primary time series. (c) As in (a), but for Temperature primary time series.</figcaption></figure>
<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /></figure>
731 **Extended Data Figure 7.** Standardized 200-year $\delta^{18}$O anomalies from Effective Moisture and $\delta_{\text{precip}}$ records. **a,** Composite medians for Iso2k$^{15}$ Effective Moisture and $\delta_{\text{precip}}$ driven records, as shown in Fig. 1 of the main text. Red shading denotes intervals of relatively high composite $\Delta^{18}$O$_{EM}$, and blue shading denotes intervals of relatively low composite $\Delta^{18}$O$_{EM}$. **b–f,** Isotopic anomalies in individual EM and $\delta_{\text{precip}}$ records contributing to the composites during those shaded intervals, via standardised anomaly maps for time intervals discussed in the main text (see Methods). These maps only include records with data spanning $\ge$ 600 years.

738 Standardized anomalies at each site are relative to the Common Era mean value for that record.

739 Maps created in R using coastlines from Natural Earth.

740

741

742

743

744

745

746

747
<figure><img src="image_17.png" /></figure>

**Extended Data Figure 8.** First Principal Component (PC1) of Iso2k data¹⁵ for the Historical Period. Symbols, lines, and shading are as per Extended Data Fig. 3, but for the interval 1850-2005, with 3-year bins. $\delta_{\text{precip}}$ PC1 explains 12% of the total variance (n = 109). Effective Moisture PC1 explains 32% of the total variance (n = 29). Temperature PC1 explains 38% of the total variance (n = 27). Maps created in R using coastlines from Natural Earth.
<figure><img src="image_18.png" /><figcaption>Extended Data Figure 9. Significance of Iso2k¹⁵ PCA eigenvalues with respect to a null hypothesis of stochastic forcing with decadal persistence. a, Block bootstrap results from 850-1840 with a 30-year bin width, 10-year block length, and tolerance of up to 15% missing data (i.e., 85% coverage during the time interval; Methods). Yellow symbols depict the eigenvalues (expressed as a percentage of the total variance) of each principal component. Blue bars show the 1-99% confidence intervals of the stochastic null hypothesis (n=1000). Eigenvalues above the 99% confidence interval are significant at the 1% level (one-sided test) and therefore are unlikely to have arisen stochastically, and unlikely to be an artifact of the data processing steps (i.e., binning and interpolation). (b) as in (a) but for 0-1980. (c) as in (a) but for 1850-2000 and with a 10-year bin width and tolerance of up to 10% missing data (90% coverage).</figcaption></figure>
774
775
776 REFERENCES
777
778
779 1. PAGES 2k Consortium. Consistent multidecadal variability in global temperature reconstructions and simulations over the Common Era. *Nature Geoscience* **12**, 643–649 (2019).
780
781
782 2. PAGES 2k Consortium. A global multiproxy database for temperature reconstructions of the Common Era. *Sci Data* **4**, 170088 (2017).
783
784 3. Field, C. B. et al. Part A: Global and sectoral aspects. Contribution of Working Group II to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change. in *Climate Change 2014: Impacts, Adaptation, and Vulnerability* (Cambridge University Press, 2014).
785
786
787 4. Smerdon, J. E. et al. Comparing proxy and model estimates of hydroclimate variability and change over the Common Era. *Proxy Use-Development-Validation/Terrestrial Archives/Centennial-Decadal* (2017) doi:10.5194/cp-2017-37-RC1.
788
789
790 5. Dansgaard, W. Stable isotopes in precipitation. *Tellus* **16**, 436–468 (1964).
791
792
793 6. Rozanski, K. Isotopes in Atmospheric Moisture. in *Isotopes in the Water Cycle* (eds. Aggarwal, P. K., Gat, J. R. & Froehlich, K. F.) 291–302 (Springer, 2005). doi:10.1007/1-4020-3023-1_18.
794
795
796
797
798 7. Bowen, G. J., Cai, Z., Fiorella, R. P. & Putman, A. L. Isotopes in the Water Cycle: Regional- to Global-Scale Patterns and Applications. *Annual Review of Earth and Planetary Sciences* vol. 47 453–479 Preprint at https://doi.org/10.1146/annurev-earth-053018-060220 (2019).
8. Galewsky, J. et al. Stable isotopes in atmospheric water vapor and applications to the
799 hydrologic cycle. *Reviews of Geophysics* **54**, 809–865 (2016).

800 9. Konecky, B. L., Noone, D. C. & Cobb, K. M. The Influence of Competing Hydroclimate

801 Processes on Stable Isotope Ratios in Tropical Rainfall. *Geophys. Res. Lett.* **46**, 1622–1633

802 (2019).

803 10. Jasechko, S. Global isotope hydrogeology—review. *Rev. Geophys.* **57**, 835–965 (2019).

804 11. Sturm, C., Zhang, Q. & Noone, D. An introduction to stable water isotopes in climate models: benefits of forward proxy modelling for paleoclimatology, Clim. Past, 6, 115--129.

806 Preprint at (2010).

807 12. Gat, J. *Isotope Hydrology: A Study of the Water Cycle*. vol. 6 (Imperial College Press,

808 2010).

809 13. Craig, H. & Gordon, L. I. Deuterium and oxygen 18 variations in the ocean and the marine atmosphere. in *Stable Isotopes in Oceanographic Studies and Paleotemperatures* (ed.

811 Tongiorgi, E.) 9–130 (Consiglio Nazionale Delle Ricerche, Laboratorio di Geologia
812 Nucleare - Pisa, 1965).

813 14. Rohling, E. J. PALEOCEANOGRAPHY, PHYSICAL AND CHEMICAL PROXIES |
814 Oxygen Isotope Composition of Seawater. in *Encyclopedia of Quaternary Science (Second
815 Edition)* (eds. Elias, S. A. & Mock, C. J.) 915–922 (Elsevier, 2013). doi:10.1016/B978-0-
816 444-53643-3.00293-4.

817 15. Konecky, B. L. et al. The Iso2k database: a global compilation of paleo-δ¹⁸O and δ²H records to aid understanding of Common Era climate. *Earth Syst. Sci. Data* **12**, 2261–2288 (2020).

820 16. Urey, H. C. Oxygen isotopes in nature and in the laboratory. *Science* **108**, 489–496 (1948).

821 17. Kaufman, D. et al. Holocene global mean surface temperature, a multi-method
reconstruction approach. *Sci Data* **7**, 201 (2020). 822

18. Brady, E. et al. The Connected Isotopic Water Cycle in the Community Earth System Model Version 1. *Journal of Advances in Modeling Earth Systems* vol. 11 2547–2566 Preprint at https://doi.org/10.1029/2019ms001663 (2019). 823 824 825

19. Stevenson, S. et al. Volcanic Eruption Signatures in the Isotope-Enabled Last Millennium Ensemble. *Paleoceanography and Paleoclimatology* **34**, 1534–1552 (2019). 826 827

20. Otto-Bliesner, B. L. et al. Climate Variability and Change since 850 CE: An Ensemble Approach with the Community Earth System Model. *Bull. Am. Meteorol. Soc.* **97**, 735–754 (2016). 828 829 830

21. Thompson, D. M. et al. Identifying hydro-sensitive coral $\delta^{18}$O records for improved high-resolution temperature and salinity reconstructions. *Geophys. Res. Lett.* **49**, (2022). 831 832

22. Dansgaard, W. Stable isotopes in precipitation. *Tellus* **16**, 436–468 (1964). 833

23. Rozanski, K., Araguás-Araguás, L. & Gonfiantini, R. Isotopic patterns in modern global precipitation. in *Climate Change in Continental Isotopic Records* (eds. Stewart, P. K., Lohmann, K. C., McKenzie, J. & Savin, S.) vol. 78 1–36 (American Geophysical Union, 1993). 834 835 836 837

24. Petit, J. R. et al. Climate and atmospheric history of the past 420,000 years from the Vostok ice core, Antarctica. *Nature* **399**, 429–436 (1999). 838 839

25. Guan, J. et al. Understanding the temporal slope of the temperature-water isotope relation during the deglaciation using isoCAM3: The slope equation. *J. Geophys. Res.* **121**, (2016). 840 841

26. Dansgaard, W. Stable isotopes in precipitation. *Tellus* **16**, 436–468 (1964). 842

27. Majoube, M. Fractionnement en oxygene 18 et en deuterium entre l'eau et sa vapeur. *Journal de Chimie Physique* **68**, 1423–1436 (1971). 843 844
28. Craig, H., Gordon, L. I. & Horibe, Y. Isotopic exchange effects in the evaporation of water: 1. Low-temperature experimental results. *J. Geophys. Res. C: Oceans* **68**, 5079–5087 (1963).

29. Horita, J., Rozanski, K. & Cohen, S. Isotope effects in the evaporation of water: a status report of the Craig-Gordon model. *Isotopes Environ. Health Stud.* **44**, 23–49 (2008).

30. Dansgaard, W. Stable isotopes in precipitation. *Tellus* **16**, 436–468 (1964).

31. Nusbaumer, J., Wong, T. E., Bardeen, C. & Noone, D. Evaluating hydrological processes in the Community Atmosphere Model Version 5 (CAM5) using stable isotope ratios of water. *J. Adv. Model. Earth Syst.* **9**, 949–977 (2017).

32. Aggarwal, P. K. et al. Stable isotopes in global precipitation: A unified interpretation based on atmospheric moisture residence time. *Geophys. Res. Lett.* **39**, L11705 (2012).

33. Held, I. M. & Soden, B. J. Robust responses of the hydrological cycle to global warming. *J. Clim.* **19**, 5686–5699 (2006).

34. Gat, J. R., Mook, W. G. & Meijer, H. A. *Environmental Isotopes in the Hydrological Cycle*, 2nd ed. vol. 2 (International Atomic Energy Agency, Vienna, Austria, 2001).

35. Siler, N. et al. The large-scale, long-term coupling of temperature, hydrology, and water isotopes. *J. Clim.* 1–51 (2021) doi:10.1175/jcli-d-20-0563.1

36. Bailey, A., Posmentier, E. & Feng, X. Patterns of evaporation and precipitation drive global isotopic changes in atmospheric moisture. *Geophys. Res. Lett.* **45**, 7093–7101 (2018).

37. Neukom, R., Steiger, N., Gómez-Navarro, J. J., Wang, J. & Werner, J. P. No evidence for globally coherent warm and cold periods over the preindustrial Common Era. *Nature* **571**, 550–554 (2019).

38. Falster, G., Konecky, B., Madhavan, M., Stevenson, S. & Coats, S. Imprint of the Pacific
868 Walker Circulation in Global Precipitation δ18O. *J. Clim.* **34**, 8579–8597 (2021).

869 39. Putman, A. L., Bowen, G. J. & Strong, C. Local and regional modes of hydroclimatic
870 change expressed in modern multidecadal precipitation oxygen isotope trends. *Geophys.
871 Res. Lett.* **48**, (2021).

872 40. Stuiver, M., Grootes, P. M. & Braziunas, T. F. The GISP2 δ18O Climate Record of the Past
873 16,500 Years and the Role of the Sun, Ocean, and Volcanoes. *Quat. Res.* **44**, 341–354
874 (1995).

875 41. Vuille, M. *et al*. A review of the South American Monsoon history as recorded in stable
876 isotopic proxies over the past two millennia. *Clim. Past* **8**, 1309–1321 (2012).

877 42. Zhao, K. *et al*. A high-resolved record of the Asian Summer Monsoon from Dongge Cave,
878 China for the past 1200 years. *Quat. Sci. Rev.* **122**, 250–257 (2015).

879 43. Konecky, B. L. *et al*. Intensification of southwestern Indonesian rainfall over the past
880 millennium. *Geophys. Res. Lett.* **40**, 386–391 (2013).

881 44. Cook, E. R., Seager, R., Cane, M. A. & Stahle, D. W. North American drought:
882 Reconstructions, causes, and consequences. *Earth-Sci. Rev.* **81**, 93–134 (2007).

883 45. Shuman, B. N. *et al*. Placing the Common Era in a Holocene context: millennial to
884 centennial patterns and trends in the hydroclimate of North America over the past 2000
885 years. *Clim. Past* **14**, 665–686 (2018).

886 46. Li, J., Xie, S.-P. & Cook, E. R. El Niño phases embedded in Asian and North American
887 drought reconstructions. *Quat. Sci. Rev.* **85**, 20–34 (2014).

888 47. Bjerknes, J. Atmospheric teleconnections from the equatorial Pacific. *Mon. Weather Rev.*
889 **97**, 163–172 (1969).

890 48. Kong, W. & Chiang, J. C. H. Southward shift of westerlies intensifies the east Asian early
891 summer rainband following El Niño. *Geophys. Res. Lett.* **47**, (2020).

892 49. Seager, R. *et al*. Adjustment of the atmospheric circulation to tropical Pacific SST anomalies: Variability of transient eddy propagation in the Pacific-North America sector.

893
Quart. J. Roy. Meteor. Soc. **136**, 277–296 (2010).

894
Quart. J. Roy. Meteor. Soc. **136**, 277–296 (2010).

895 50. Wittenberg, A. T. Are historical records sufficient to constrain ENSO simulations? *Geophys.*

896
Res. Lett. **36**, L12702 (2009).

897 51. Collins, M. *et al*. The impact of global warming on the tropical Pacific Ocean and El Niño.

898
Nat. Geosci. **3**, 391–397 (2010).

899 52. Vecchi, G. *et al*. Weakening of tropical Pacific atmospheric circulation due to anthropogenic

900
forcing. *Nature* **441**, 73–76 (2006).

901 53. Konecky, B. L. & McKay, N. P. NOAA/WDS Paleoclimatology - The Iso2k Database.

902
NOAA National Centers for Environmental Information. (2020) doi:10.25921/57j8-vs18.

903 54. Craig, H. Isotopic Variations in Meteoric Waters. *Science* **133**, 1702–1703 (1961).

904 55. Beckers, J.-M., Barth, A. & Alvera-Azcárate, A. DINEOF reconstruction of clouded images

905
including error maps–application to the Sea-Surface Temperature around Corsican Island.

906
*Ocean Science* **2**, 183–199 (2006).

907 56. Alvera-Azcárate, A., Barth, A., Sirjacobs, D., Lenartz, F. & Beckers, J.-M. Data

908
Interpolating Empirical Orthogonal Functions (DINEOF): a tool for geophysical data

909
analyses. *Mediterr. Mar. Sci.* **12**, 5–11 (2011).

910 57. Oksanen, J., Kindt, R., Legendre, P. & Others. Vegan: community ecology package. R

911
package version 2.2-1. R Development Core Team, Vienna. Preprint at (2015).

912 58. Wong, T. E., Nusbaumer, J. & Noone, D. C. Evaluation of modeled land-atmosphere

913
exchanges with a comprehensive water isotope fractionation scheme in version 4 of the
914 Community Land Model. *J. Adv. Model. Earth Syst.* **9**, 978–1001 (2017).

915 59. Allan, R. & Ansell, T. A New Globally Complete Monthly Historical Gridded Mean Sea

916 Level Pressure Dataset (HadSLP2): 1850–2004. *J. Clim.* **19**, 5816–5842 (2006).