# Past terrestrial hydroclimate sensitivity controlled by Earth system feedbacks

Ran Feng, Tripti Bhattacharya, Bette Otto-Bliesner, Esther Brady, Alan Haywood, Julia Tindall, Stephen Hunter, Ayako Abe-Ouchi, Wing-Le Chan, Masa Kageyama, et al.

► To cite this version:

Ran Feng, Tripti Bhattacharya, Bette Otto-Bliesner, Esther Brady, Alan Haywood, et al. Past terrestrial hydroclimate sensitivity controlled by Earth system feedbacks. Nature Communications, 2022, 13 (1), pp.1306. 10.1038/s41467-022-28814-7. hal-03620581

HAL Id: hal-03620581

https://hal.science/hal-03620581

Submitted on 26 Mar 2022

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>
<figure><img src="image_2.png" /></figure>

ARTICLE

https://doi.org/10.1038/s41467-022-28814-7

OPEN

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

Check for updates

# Past terrestrial hydroclimate sensitivity controlled by Earth system feedbacks

Ran Feng<fnref n="1" />, Tripti Bhattacharya², Bette L. Otto-Bliesner³, Esther C. Brady³, Alan M. Haywood<fnref n="4" />, Julia C. Tindall<fnref n="4" />, Stephen J. Hunter<fnref n="4" />, Ayako Abe-Ouchi⁵, Wing-Le Chan⁵, Masa Kageyama⁶, Camille Contoux⁶, Chuncheng Guo⁷, Xiangyu Li⁸, Gerrit Lohmann⁹,¹⁰, Christian Stepanek⁹, Ning Tan¹¹, Qiong Zhang¹², Zhongshi Zhang⁸, Zixuan Han¹³, Charles J. R. Williams¹⁴, Daniel J. Lunt¹⁴, Harry J. Dowsett¹⁵, Deepak Chandan¹⁶ & W. Richard Peltier¹⁶

Despite tectonic conditions and atmospheric CO₂ levels (pCO₂) similar to those of present-day, geological reconstructions from the mid-Pliocene (3.3-3.0 Ma) document high lake levels in the Sahel and mesic conditions in subtropical Eurasia, suggesting drastic reorganizations of subtropical terrestrial hydroclimate during this interval. Here, using a compilation of proxy data and multi-model paleoclimate simulations, we show that the mid-Pliocene hydroclimate state is not driven by direct CO₂ radiative forcing but by a loss of northern high-latitude ice sheets and continental greening. These ice sheet and vegetation changes are long-term Earth system feedbacks to elevated pCO₂. Further, the moist conditions in the Sahel and subtropical Eurasia during the mid-Pliocene are a product of enhanced tropospheric humidity and a stationary wave response to the surface warming pattern, which varies strongly with land cover changes. These findings highlight the potential for amplified terrestrial hydroclimate responses over long timescales to a sustained CO₂ forcing.

<footnote n="1">¹Department of Geosciences, College of Liberal Arts and Sciences, University of Connecticut, Mansfield, CT, USA. ²Department of Earth and Environmental Sciences, Syracuse University, Syracuse, NY, USA. ³Climate and Global Dynamics Laboratory, National Center for Atmospheric Research, Boulder, CO, USA.</footnote>

<footnote n="4">⁴School of Earth and Environment, University of Leeds, Woodhouse Lane, Leeds, West Yorkshire LS29JT, UK. ⁵Atmosphere and Ocean Research Institute, University of Tokyo, Kashiwa, Japan. ⁶Laboratoire des Sciences du Climat et de l'Environnement, LSCE/IPSL, CEA-CNRS-UVSQ, Université Paris-Saclay, F-91191 Gif-sur-Yvette, France. ⁷NORCE Norwegian Research Centre, Bjerknes Centre for Climate Research, 5007 Bergen, Norway. ⁸Department of Atmospheric Science, School of Environnemental Studies, China University of Geoscience, Wuhan 430074, China. ⁹Alfred Wegener Institute-Helmholtz Centre for Polar and Marine Research, Bremerhaven, Germany. ¹⁰University of Bremen, Bremen, Germany. ¹¹Key Laboratory of Cenozoic Geology and Environment, Institute of Geology and Geophysics, Chinese Academy of Sciences, Beijing 100029, China. ¹²Department of Physical Geography and Bolin Centre for Climate Research, Stockholm University, Stockholm, Sweden. ¹³College of Oceanography, Hohai University, Nanjing, China. ¹⁴School of Geographical Sciences and Cabot Institute, University of Bristol, University Road, Bristol BS8 1SS, UK. ¹⁵Florence Bascom Geoscience Center, U. S. Geological Survey, Reston, VA, USA. ¹⁶Department of Physics, University of Toronto, Toronto, CA, Canada. email: ran.feng@uconn.edu</footnote>
Geologic evidence suggests dramatic reorganizations of subtropical climate during past greenhouse climate intervals, including the mid-Piacenzian Warm Period (3.3–3.0 Ma, commonly referred to as the mid-Pliocene). Multiple proxies of hydroclimate indicate large, deep lakes and reduced dust flux across North Africa during the mid-Pliocene¹,², and more mesic vegetation in South and East Asia³,⁴. Additional sedimentary and paleobotanical data also points to wetter subtropical Eurasia conditions prior to the intensification of northern hemisphere glaciation (Supplementary Table 1). These changes imply large increases in precipitation (P) minus evaporation (E), referred to as P–E. All are associated with pCO₂ levels of approximately 400 ppm⁵,⁶, similar to today’s level. Elevated pCO₂ can lead to moderate increases in precipitation across these regions as a result of tropospheric moistening⁷, enhanced land–sea thermal contrast⁸, and enhanced inland moisture advection⁹,¹⁰. However, evaporation also increases due to surface warming. As a result, predicted changes in terrestrial water balance (P–E), and associated changes in soil moisture and runoff remain equivocal across subtropical continents¹¹.

Modeled hydroclimate responses to CO₂ increases broadly follow the “wet-gets-wetter, dry-gets-drier” paradigm⁷ with strong modulations from land warming patterns, the land-sea warming contrast¹², tropical SST patterns¹³, and feedbacks from soil moisture and CO₂ fertilization effects on leaf phenology¹⁴,¹⁵. In simulations featuring middle-of-the-road future warming scenarios, predicted changes in P–E by the end of the 21st century (2081–2100) are minimal across the subtropical Sahel and East Asia relative to both the historical period (1986–2005) (Fig. 1a and d) and preindustrial (PI) (Supplementary Fig. 1). These future scenarios are often thought to be comparable to the mid-Pliocene climate¹⁶,¹⁷(Fig. 1a). Sources for the disparity between the hydroclimate state recorded by mid-Pliocene proxies and future simulations are unknown. One potential source might be the transient nature of future climate change as opposed to the quasi-equilibrium nature of mid-Pliocene climate. Yet, even equilibrium simulations with only the mid-Pliocene CO₂ forcing fail to produce moist subtropical terrestrial conditions¹⁸. Our simulation also confirms this result (Fig. 1).

From the perspective of atmospheric dynamics, two leading hypotheses have been proposed to explain the wetter subtropical continents during past warm climates. Both hypotheses were proposed when many older-generation models at lower spatial resolutions were not able to simulate wetter subtropical continents. One hypothesis emphasizes the hydroclimate impact of an El Niño-like Pacific mean state. SST records from the tropical Pacific document greater warming across the eastern equatorial Pacific than the western Pacific warm pool¹⁹⁻²³, resulting in an El-Niño-like pattern of SST anomalies. An El Niño-like Pacific SST pattern may strengthen and shift the subtropical jet equatorward, enhancing the transient eddy-driven moisture convergence and ascent²⁴,²⁵. The other hypothesis focuses on the role of a sluggish Hadley Circulation that results from a relaxed meridional SST gradient²⁶⁻²⁸. A weaker Hadley Circulation may lead to weakened zonal mean moisture divergence from the subtropics and, in turn, reduced aridity¹⁸,²⁹,³⁰. Nonetheless, with advancements in model development and boundary conditions, newer earth system model simulations show substantial differences in simulated past climate states²³,³¹,³², particularly in mid-Pliocene SST patterns³³, polar warmth³⁴, Atlantic Overturning Circulation³⁵, and precipitation³⁶.

Here, using atmosphere-ocean coupled global climate model (GCM) simulations from the most recent Pliocene Model Intercomparison Project Phase II (PlioMIP2)³³,³⁷,³⁸, we demonstrate that most current generation GCMs can reproduce the pattern of mid-Pliocene hydroclimate of the subtropical Sahel and East Asia

suggested by proxies without any paleoclimate-specific changes to model parameterizations. We focus on these regions given the large, coherent signal across the majority of PlioMIP2 simulations and the convergence between proxy data and model simulations (Fig. 1b, c). We further develop several new simulations using the Community Earth System Model version 2³⁵,³⁹ to explore the extent to which simulated mid-Pliocene hydroclimate changes can be generated by changes in CO₂, tectonics, or vegetation and ice sheets. In contrast to previous hypotheses, wetter mid-Pliocene conditions in the subtropical Sahel and East Asia are driven by tropospheric moistening and changes to stationary wave dynamics in response to surface warming patterns due to vegetation and ice sheet changes. Both changes are part of the long-term Earth system feedbacks to a sustained CO₂ forcing. Consequently, model skill at simulating Pliocene hydroclimate states strongly scales with earth system sensitivities (ESS) of individual models instead of the equilibrium climate sensitivities (ECS).

## Results

**δ(P–E) pattern in models and proxies.** The last 100 years of simulations by 13 PlioMIP2 GCMs (Supplementary Table 2) were averaged to produce the ensemble mean. Mid-Pliocene changes in P–E, referred to as δ(P–E), and other climate variables are calculated with respect to preindustrial (PI) values. A robust moistening signal that is larger than the intermodel variability is found across the Sahel and subtropical Eurasia (Fig. 1b). This pattern is most pronounced during the boreal summer months (June to September) (Fig. 1c), with little or opposite changes during the boreal winter months (December to March) (Fig. S2). The spatial continuity of positive δ(P–E) from North Africa to subtropical East Asia is not a visual coincidence: models that show a large precipitation increase in the Sahel also tend to show a large precipitation increase in subtropical East Asia (Fig. 1d), indicating similar processes driving hydroclimate changes in both regions. We confirm this with moisture budget diagnostics (see below).

To compare modeled patterns of hydroclimate change to available geologic data, we compiled proxy indicators of mid-Pliocene terrestrial hydroclimate, drawing on existing compilations¹⁸,⁴⁰,⁴¹ as well as our own literature search. We identify a total of 64 proxy records that include sedimentological indicators, palynological, floral or faunal, offshore marine records, and stable isotope analyses of organic and inorganic materials (Table S1 and Fig. S3). These records are interpreted as qualitative or semi-quantitative indicators of mean hydroclimate state given the lack of orbital-scale variabilities documented in these records (see Methods). We quantify the extent to which proxies and models produce the same patterns of wetter, drier, or unchanged mid-Pliocene hydroclimate using a metric known as Gwet’s AC designed for categorical data (see Methods). To account for the unknown sensitivity of proxy hydroclimate indicators to Pliocene P–E changes, model values of δ(P–E) are expressed as % changes from PI values at the proxy sites.

The annual pattern of multi-model mean mid-Pliocene δ(P–E) shows statistically significant agreement with proxy indicators of hydroclimate across the subtropical Sahel and East Asia for a wide range of δ(P–E) thresholds (Fig. 2a). This agreement is strongly driven by the simulated summer δ(P–E) pattern (Fig. 2c). Additionally, a subset of experiments featuring dynamic phenology and terrestrial carbon cycle also show a strong coupling between the positive δ(P–E) and enhanced net-primarily productivity across both regions, consistent with paleoecological reconstructions (Fig. S5). Yet, different models show varying skills at capturing proxy hydroclimate changes (Fig. 2a). For instance,
<figure><img src="image_4.png" /><figcaption>**Fig. 1** Hydroclimate changes in mid-Pliocene proxy records and simulations of mid-Pliocene and future climate. **a** Changes in precipitation minus evaporation (δ(P-E), mm day⁻¹) between 2081 to 2100 and 1986 to 2005 following Shared Socioeconomical Pathway (SSP) 2-4.5 simulated by models participating in Climate Model Intercomparison Project 6. **b, c** Mid-Pliocene annual and boreal summer (June to September) mean δ(P-E) of the PlioMIP2 ensemble relative to the preindustrial simulated by the same models. Proxy data displayed as filled circles. **d** Correspondence between the subtropical Sahel (10°-20°N, 10°W-25°E) and East Asia (20°-30°N, 80°E-100°E) δ(P-E) simulated by PlioMIP2 experiments and a subset of SSP2-4.5 experiments with the same models (model names are marked with asterisks). **e-h** δ(P-E) in response to full Pliocene climate forcing conditions (Fall), CO₂ forcing alone (FCO₂), changes in geography and topography (Fgeotop), and changes in vegetation and icesheet (Fvegice) simulated by Community Earth System Model version 2. The SSP2-4.5 ensemble includes models of BCC-CSM2-MR, CESM2, CESM2-WACCM, CanESM5, CNRM-CM6-1, CNRM-ESM2-1, EC-Earth3.3, GISS-E2-1-G, HadGEM3-GC31-LL, IPSL-CM6A-LR, MIROC6, MIROC-ES2L, MRI-ESM2-0, NESM3, and UKESM1-0-LL. Area significant against the multi-model spread is hatched in (a)-(c) identified by Welch's t-test (*p* < 0.1).</figcaption></figure>
<figure><img src="image_5.png" /></figure>

**Fig. 2** Degree of agreement between mid-Pliocene proxy hydroclimate indicators and simulations. Proxy-model fit is assessed using a measure of categorical agreement between two datasets called Gwet's AC (unitless). For a given % threshold of simulated changes in precipitation minus evaporation ($\delta(P-E)$) relative to the preindustrial, higher (lower) values indicate that proxies and models agree (disagree) that given locations are wet, dry, or neutral. The area within the dashed line indicates statistically significant agreements. **a** Gwet's AC agreement between individual models, Multimodel mean (MMM), and proxies at different thresholds. CMIP6 models are identified with asterisks. **b** Agreement between proxies and simulated $\delta(P-E)$ in response to single or combined climate forcings ($F_{\text{all}}$) from changes of CO$_2$ ($F_{\text{CO2}}$), tectonics ($F_{\text{geotop}}$), and vegetation and ice sheets ($F_{\text{vegice}}$). **c** Gwet's AC agreement between annual and boreal summer averages at proxy sites for PlioMIP2 ensemble mean.

IPSL-CM6, and EC-Earth3.3, two models with the highest level of agreement between mid-Pliocene proxies and simulations, show expansive inland wetter conditions across North Africa, Mediterranean, and subtropical East Asia (Fig. S6). In contrast, NorESM-L and GISS-E2-1-G, two of the models with the lowest proxy-model agreement, show highly mixed or muted $\delta(P-E)$ across this region. Similar to the multi-model mean result, the agreement between individual models and proxy hydroclimate indicators is also strongly driven by the simulated summer $\delta(P-E)$ (Fig. S4).

**CO$_2$ or boundary conditions in driving positive $\delta(P-E)$.** Single forcing experiments are commonly used to quantify climate responses to individual forcing agent$^{42}$. In our case, three sets of simulations are constructed with CESM2$^{39}$ to quantify contributions to $\delta(P-E)$ from a range of mid-Pliocene forcings: a 400 ppm CO$_2$ ($F_{\text{CO2}}$), changes in biome distribution and ice sheets ($F_{\text{vegice}}$), and changes in geography and topography ($F_{\text{geotop}}$) (Methods). The separation of $F_{\text{vegice}}$ and $F_{\text{CO2}}$ is designed to separate the influences of vegetation and ice sheet changes from the direct influences of CO$_2$ changes. The former represents Earth system feedbacks, which are not typically considered when evaluating equilibrium climate responses to CO$_2$ forcing$^{43}$. This experimental scheme assumes decorrelation between climate responses to $F_{\text{CO2}}$ and $F_{\text{geotop}}$, which may not be the case for paleoclimate conditions because various feedbacks have been shown to depend on the background climate warmth at high CO$_2$ levels$^{44,45}$. However, our simulations support this decorrelation under moderate $F_{\text{CO2}}$ and $F_{\text{geotop}}$. Simulated responses to $F_{\text{geotop}}$ are consistent with modern or mid-Pliocene pCO$_2$, and simulated responses to $F_{\text{CO2}}$ are consistent with modern or mid-Pliocene topography and geography (Fig. S7). For $F_{\text{vegice}}$ land surface and vegetation changes considered here are a consequence of mid-Pliocene $F_{\text{CO2}}$ and $F_{\text{geotop}}$$^{43}$. Therefore, $F_{\text{vegice}}$ is separated from the combined $F_{\text{CO2}}$ and $F_{\text{geotop}}$.

$F_{\text{vegice}}$ explains ~78% (2.2 mm day$^{-1}$) of the regional mean $\delta(P-E)$ induced by $F_{\text{all}}$ (2.8 mm day$^{-1}$) across the subtropical Sahel and East Asia, while contributions from $F_{\text{CO2}}$ and $F_{\text{geotop}}$ are small (Fig. 1f, g). Furthermore, only $F_{\text{vegice}}$ produces a similar level of proxy-model agreement in $\delta(P-E)$ compared to $F_{\text{all}}$ (Fig. 2b). Both $F_{\text{CO2}}$ and combined $F_{\text{geotop}}$ and $F_{\text{CO2}}$ produce low values of Gwet's AC (Fig. 2b). The proxy-model agreement seen in the full forcing ($F_{\text{all}}$) experiment is only reproducible in the experiment forced with $F_{\text{vegice}}$.

Compared to $F_{\text{all}}$, $F_{\text{vegice}}$ accounts for ~50% of the increase in the global mean surface temperature (2.7 °C) (Fig. 3) and 58% of the increase (0.45 g/kg) in the global mean tropospheric (100 to 1000 hPa) specific humidity. In the Northern subtropics (between 20°N and 30°N), $F_{\text{vegice}}$ explains an even greater fraction of tropospheric moistening (61%, or 0.56 g/kg). In contrast, $F_{\text{CO2}}$ accounts for 45% of the global mean surface warming (2.4 °C) (Fig. S11), 31% (0.24 g/kg) of the global mean tropospheric moistening, and only 26% of the tropospheric moistening in the northern subtropics. $F_{\text{geotop}}$ has much smaller influences on temperature and moisture responses globally. The influence of $F_{\text{geotop}}$ is slightly elevated in the northern subtropics accounting for 13% of the increase (0.11 g/kg) in the tropospheric humidity. Warming due to $F_{\text{vegice}}$ is attributable to both lowered surface albedo and enhanced evapotranspiration. Areas where boreal forest shifts and expands northward and where mid-latitude deserts become vegetated (Fig. 3a) feature substantially lowered surface albedo (Fig. 3b). Expanded boreal forests also show large increases in upward latent heat flux, suggesting enhanced water vapor feedback$^{46}$ (Fig. 3c). A similar increase in latent heat flux also occurs in the northern subtropics where $\delta(P-E)$ is positive.

**The dynamical linkage between climate forcing conditions and mid-Pliocene hydroclimate.** In order to identify the dynamical linkage between climate forcing conditions and hydroclimate in the subtropical Sahel and East Asia, we adapt previously published moisture budget diagnostics (MBD)$^{47}$ to decompose simulated June to September $\delta(P-E)$ by individual models (Methods) into changes in the seasonal cycle of tropospheric Sahel and East Asia, while contributions from $F_{\text{CO2}}$ and $F_{\text{geotop}}$ are small (Fig. 1f, g). Furthermore, only $F_{\text{vegice}}$ produces a similar level of proxy-model agreement in $\delta(P-E)$ compared to $F_{\text{all}}$ (Fig. 2b). Both $F_{\text{CO2}}$ and combined $F_{\text{geotop}}$ and $F_{\text{CO2}}$ produce low values of Gwet's AC (Fig. 2b). The proxy-model agreement seen in the full forcing ($F_{\text{all}}$) experiment is only reproducible in the experiment forced with $F_{\text{vegice}}$.
<figure><img src="image_6.png" /><figcaption>**Fig. 3** Climate responses to the full mid-Pliocene climate forcing ($F_{\text{all}}$) and forcing from vegetation and ice sheet changes ($F_{\text{vegice}}$). **a** vegetation and ice sheet distribution: the plant functional type (PFT) with the highest percentage in a grid cell is shown. Changes of **b** surface albedo ($\delta a, \%$) and **c** surface upward latent heat flux ($\delta \text{LHF}, \text{W m}^{-2}$) induced by $F_{\text{vegice}}$. Label abbreviations in **a** are, NE needleleaf evergreen, BE broadleaf evergreen, BD broadleaf deciduous, Temp temperate, Bore boreal, Deci deciduous, Trop tropical. **d–f** Simulated responses of surface temperature ($\Delta T, ^{\circ}\text{C}$) (color shaded), and wave number 1 of the 600 hPa stationary wave ($\psi_{\alpha}$ contour, $\text{m}^2 \text{s}^{-1}$) and corresponding rotational winds (vectors, $\text{m s}^{-1}$) to **d** all forcings ($F_{\text{all}}$) in the PlioMIP2 ensemble, **e** $F_{\text{all}}$ in CESM2, and **f** $F_{\text{vegice}}$ in CESM2. The diagnosed thermal wind directions based on the $\Delta T$ pattern are also shown.</figcaption></figure>

moisture content ($\delta(P-E)_t$), changes in the zonal mean circulation dynamics ($\delta(P-E)_{[V]}$) and stationary wave dynamics ($\delta(P-E)_{Q^*}$), changes in the tropospheric moisture content ($\delta(P-E)_Q$), changes in the interactions between moisture and circulation dynamics ($\delta(P-E)_{VQ}$), and a residual term ($\delta(P-E)_{\text{RES}}$) (Fig. 4 and Fig. S8). The stationary wave response is quantified as the temporal mean departure from the zonal mean following the classic circulation decomposition$^{48}$. The residual term combines the effects of variability of transient eddies and changes of topography (see Methods).

As revealed by MBD, positive $\delta(P-E)$ in the subtropical Sahel and East Asia primarily arises from changes in stationary wave dynamics ($\delta(P-E)_{V^*}$) and increased tropospheric moisture content ($\delta(P-E)_Q$) (Fig. 4b, c). Contributions from both $\delta(P-E)_t$ and $\delta(P-E)_{VQ}$ are insignificant (Fig. S8). Intermodel spread in simulated $\delta(P-E)$ strongly scales with the spread in $\delta(P-E)_{V^*}$ and $\delta(P-E)_Q$, with little dependency on other terms (Fig. S10). Moreover, the MBD results are consistent between the PlioMIP2 ensemble mean, and CESM2 experiments with $F_{\text{all}}$ and $F_{\text{vegice}}$ (Fig. 4e, f, and Fig. S8). These results suggest a strong dynamical linkage between $\delta(P-E)_{V^*}$, $\delta(P-E)_Q$, and $\delta(P-E)$ across PlioMIP2 models and that this dynamical linkage can mostly be reproduced with $F_{\text{vegice}}$.

We further decompose the stationary wave response into components associated with different wave numbers through Fourier decomposition of stream function anomalies ($\psi_a$) at 600 hPa. $\psi_a$ is calculated as departures from the zonal mean and PI (Fig. 3d–f). As shown by the Fourier decomposition, zonal wavenumber 1 of $\psi_a$ accounts for 67% of the total wave energy between 0 and 90°N. This wave component features cyclones separately centered over western Europe and North Africa, and anticyclones centered over the North Pacific (Fig. 3d–f). Following the southern edge of the cyclical wave centers, a moisture transport corridor emerges: westerly winds bring moisture to North Africa from the tropical Atlantic; southerly and southerly winds bring moisture from the tropical Indian and Pacific Oceans towards the Indian subcontinent and East Asia. Furthermore, the rotational winds of this wave component are connected to the surface warming pattern via the thermal wind

relationship. This pattern is broadly reproduced in the simulation forced with $F_{\text{vegice}}$ (Fig. 3f).

Positive $\delta(P-E)$ in both regions also results from increased tropospheric moisture content ($\delta(P-E)_Q$). Under PI conditions, low-level winds converge towards North Africa and subtropical East Asia, with the diverging flow in the adjacent regions during the boreal summer (Fig. 4c). This circulation is a key feature of the regional summer monsoon$^{49}$. Even without changing this circulation, elevated tropospheric moisture content can result in greater moisture convergence and positive $\delta(P-E)_Q$ in the subtropical Sahel and East Asia, and greater moisture divergence and negative $\delta(P-E)_Q$ in the adjacent regions. This $\delta(P-E)_Q$ pattern is a known signature of thermodynamic response to elevated CO$_2$, i.e., the wet-gets-wetter, dry-gets-drier paradigm$^{7,12,50}$. As shown in $F_{\text{vegice}}$, a similar thermodynamic response can be induced through land cover changes.

## Discussion

In our simulations, mechanisms causing the positive $\delta(P-E)$ in the subtropical Sahel and East Asia are different from previous studies. The MBD reveals that both $\delta(P-E)_{[V]}$ and $\delta(P-E)_{\text{RES}}$ are small (Fig. 4a, d). A small $\delta(P-E)_{[V]}$ suggests little contribution from the zonal mean Hadley Circulation change. If Hadley Circulation change were to drive boreal summer $\delta(P-E)$, coherent $\delta(P-E)$ across longitudes would be expected$^{18}$. This is not the case for PlioMIP2 simulations (Fig. 1b, and Fig. S6). Most PlioMIP2 simulations display a clear latitudinal offset in positive $\delta(P-E)$ between East Asia and West Pacific (Fig. S6) and between Sahel and Central America.

A small $\delta(P-E)_{\text{RES}}$ suggests a small net influence from changes of transient eddies and topography on $\delta(P-E)$ (Fig. 4d). Tropical SST warming during the El Niño events influences extratropical precipitation by strengthening and shifting the positions of storm tracks and the subtropical jet equatorward$^{25}$. These changes are not observed in mid-Pliocene simulations of the northern extratropics (Fig. S9). In this area, storm activity weakens as shown by changes in the 850 hPa eddy kinetic energy. The subtropical jet, measured as the
<figure><img src="image_7.png" /><figcaption>**Fig. 4** Moisture budget diagnostics of mid-Pliocene boreal summer changes in precipitation minus evaporation (δ(P-E), mm day⁻¹). **a-d** Contributions to δ(P-E) from zonal mean circulation (δ(P-E)ₜ), changes of stationary wave dynamics (δ(P-E)ₜ), tropospheric moistening (δ(P-E)Q), and a residual term(δ(P-E)RES). **e, f** Regional mean (red boxes in (a)-(d)) moisture budget for the subtropical Sahel and East Asia simulated by PlioMIP2 ensemble, and CESM2 with Fₐll and Fₐvegice. Error bars show 1 standard deviation of model spread.</figcaption></figure>

200 hPa zonal wind speed, also weakens across Eurasia but remains in a similar position to PI. These changes are consistent with a reduction in the poleward temperature gradient⁵¹ (Fig. S9a, b). Moreover, the pattern of precipitation variability induced by El Niño SSTs⁵² has a classic signature of negative δ(P-E) in North Africa paired with positive δ(P-E) in Southeast Asia. This signature is not observed in mid-Pliocene simulations: δ(P-E) is broadly consistent across both regions (Fig. 1b and S6). Models suggest that El Niño-like SSTs, therefore, do not play a significant role in driving mid-Pliocene δ(P-E) in these regions.

However, δ(P-E)RES is more noticeable in the North Pacific and subtropical western North America. The impact of El Niño-like SSTs could be important for δ(P-E) in both regions¹³. Past hydroclimate changes in western North America are shown to be sensitive to changes in eddy moisture transport associated with varying conditions of atmospheric rivers⁵³⁻⁵⁵. δ(P-E) changes in this region are likely driven by processes different from the subtropical Sahel and East Asia.

The stationary wave response identified here is distinct from the response of the ITCZ shift. Both the Hadley Circulation change and ITCZ shift highlight the importance of changing the zonal mean energy budget on circulation dynamics¹⁴,¹⁸,⁵⁶⁻⁶⁰. Instead, the stationary wave response in PlioMIP2 experiments is generated by a zonally heterogeneous pattern of energy perturbation and warming. Both are induced by continental greening. Circulation changes in the lower troposphere can be explained by this surface warming pattern through the diagnostic thermal winds, which closely follow the zonal wavenumber 1 of ψₐ (Fig. 3). This “pattern effect” has often been overlooked when examining past hydroclimate changes.

Why is Fₐvegice more effective at altering terrestrial hydroclimate compared to FₐCO₂ in the subtropical Sahel and East Asia? Land cover changes are known to be key to generating surface temperature and circulation responses across North Africa and subtropical East Asia⁶¹⁻⁶³ and may even alter the strength of the Atlantic meridional overturning circulation⁶⁴,⁶⁵. In North Africa, expansion of desert results in enhanced surface albedo, surface cooling, strengthened diabatic subsidence, and dust emission⁶⁶,⁶⁷. These responses may further suppress moist convection and perpetuate desert⁶⁸⁻⁷¹. A positive vegetation-precipitation feedback results in multiple equilibrium states of North African vegetation over the late Quaternary⁶². Different from those mechanisms, in our simulations, atmospheric circulation responses to Fₐvegice facilitate moisture transport towards the subtropical Sahel and East Asia via a stationary wave response. This response is closely tied to the large-scale warming pattern generated by Fₐvegice, which does not occur in FₐCO₂ (Fig. S11). Additionally, this response also differs from the regional wave response driven by diabatic heating from the South Asian monsoon moist convection as suggested by the previous studies⁶³. The latter features contrasting wave center and hydroclimate responses between South Asia and the western Sahel, distinct from our findings of a continental-scale wave pattern with consistent responses between these two regions.

Different hydroclimate responses to FₐCO₂ and Fₐvegice may also reflect different soil moisture and plant physiological responses. Increasing CO₂ may favor the reduction of soil moisture and partitioning of surface heat flux toward sensible heat, leading to enhanced surface warming. This in turn lowers relative humidity above the surface and diminishes continental cloud cover,
contributing to negative $δ(P-E)$ on land. A similar response of surface heat flux partitioning and moisture deficit can be caused by the reduction of leaf transpiration in response to CO₂ fertilization. Soil moisture feedback and CO₂ fertilization drive much of the predicted future subtropical terrestrial hydroclimate change¹⁵, and may also contribute to the muted $δ(P-E)$ response to FCO₂ in our experiments despite a modest increase in precipitation. In contrast, Fvegice features no direct physiological effect of CO₂ and produces large increases in latent heat flux relative to sensible heat flux across continents (Fig. S12), creating favorable conditions for a more humid troposphere. Humidification of the subtropical Sahel and East Asia, therefore, reflect a synergy of dynamic and thermodynamic responses to Fvegice. Reduced ice sheet cover, expansive northern high-latitude boreal forests, and the vegetated Sahel and Central Asia have all been recorded during other Cenozoic warm intervals⁷²⁻⁷⁴. These land cover changes were likely instrumental in driving global mean temperature and the terrestrial hydrological cycle throughout the Cenozoic.

Our results suggest an alternate view of the role of vegetation changes in modulating continental hydroclimate. Changes in regional circulation in the form of stationary wave responses lead to strengthened low-level winds that import moisture into the subtropical Sahel and East Asia from the tropical Atlantic and Indian Oceans, respectively. This inland moisture influx is further amplified by enhanced tropospheric humidity. This mechanism is distinct from previous mechanisms that highlight the role of the zonal mean Hadley Circulation or local ecosystem responses to land surface changes. Although complete feedbacks from vegetation and ice sheets are not prognostically simulated in most mid-Pliocene simulations (Supplementary Table 2), the presented simulations highlight radiative perturbations from proxy-constrained changes in vegetation and ice sheets are key to generating hydroclimate changes in the northern subtropics. In comparison, prescribed mid-Pliocene topography and land-sea distribution have little effect (Fig. 1g).

A key implication of our findings is that the mid-Pliocene continental hydroclimate is more appropriately viewed as part of the Earth system feedbacks instead of an immediate response to FCO₂, and hence, requires appropriate boundary conditions of vegetation and ice sheet to simulate. This inference is supported by the strong relationship between a model's ability to capture the mid-Pliocene hydroclimate state and its simulated global mean warming (Fig. 5a). The latter reflects model diversity in Earth System Sensitivity (ESS), which incorporates surface temperature responses to long-term biome range shift and ice sheet changes in addition to CO₂ changes. In contrast, the relationship between a model's ability to capture the mid-Pliocene hydroclimate state and its equilibrium climate sensitivity (ECS) is weak (Fig. 5b). ECS of an individual model is estimated from the doubling CO₂ experiment with the same model and resolution and without biome range and ice sheet changes³³. Published studies have mostly focused on the impact of Earth system feedback on temperature responses⁴³,⁷⁵. Here, we demonstrate that these feedbacks are key for understanding paleohydrological responses. Inconsistencies in interpretations of proxy hydroclimate records¹⁸,⁷⁶ may be resolved by considering Fvegice.

Lastly, our results offer a resolution to the apparent discrepancies between the projections of hydroclimate changes in the Sahel and subtropical East Asia following the middle-of-the-road scenarios and strong geologic evidence for moist Pliocene climate at similar levels of CO₂. While the near-future is dominated by the short-term response to CO₂ radiative forcing and internal variability, Pliocene hydroclimate reflects long-term adjustments of the Earth system that incorporates responses from vegetation and ice sheets, which occur on timescales longer than most future climate projections. Feedbacks of vegetation and ice sheets to CO₂ increase are known to amplify the response of

equilibrium surface temperature to radiative forcing⁴³,⁷⁵. We highlight that these relatively slow Earth system feedbacks are also critical for understanding Earth's hydroclimate responses to varying CO₂. Therefore, changes in vegetation and ice sheet distributions should be carefully considered when simulating past and future climates.

## Methods

**Proxy-model comparison.** Our proxy compilation builds on previous efforts to compile Pliocene hydroclimate records. These include the compilations created by refs. ¹⁸,⁴⁰,⁴¹. The records included in these sources span a variety of proxy types, from sedimentological indicators of lacustrine environments, palynological indicators of vegetation composition, faunal remains, and stable isotope records from organic and inorganic materials. We added to these compilations by including new records of terrestrial hydroclimate dating to the Pliocene archived on the NOAA Paleoclimatology Database (https://www.ncdc.noaa.gov/data-access/paleoclimatology-data) and Pangaea Database (https://www.pangaea.de/).

To identify the average hydroclimate signal during the mid-Pliocene, we filtered available records based on the precision of their age models. Specifically, proxy records were required to include at least two age control points, with one age control point after the mid-Pliocene and a clearly identifiable age control point prior to mid-Piacenzian. For many records, this “basal” tie point was often in the early Pliocene or Miocene. Only including records with multiple age control points reduces the likelihood that samples inferred to be from the mid-Pliocene actually come from earlier in the Pliocene or the early Pleistocene. For most proxy records, age control derives from a combination of magnetostratigraphy (e.g., the Gauss–Matuyama boundary), radiometric dating, or biostratigraphic information. In some cases, identification of independently-dated tephra layers or correlation with the benthic oxygen isotope stack provides age control. These filtering criteria allowed us to retain a compilation of 62 records. Of the 62 records included, 30 come from paleobiological indicators like faunal remains or pollen, while the other 34 records are drawn from interpretations of sedimentary sequences or stable isotopic analyses of organic and inorganic materials (SI). A supplementary excel file (Supplementary Data 1) with details on the proxy records, including methods, chronology, and original citation, is included with this manuscript. Because of the broad, continental-scale coherence of the model signal, proxy-model synthesis results are similar when using the compilations found in other sources e.g., ref. ⁴¹.

We rely on the author's original interpretation about whether the record reflects, on average, wetter or drier conditions, or no change in hydroclimate during the mid-Pliocene compared to late Quaternary/modern conditions. Many records are discontinuous, and cannot provide quantitative comparisons between mid-Pliocene climate and late quaternary or PI conditions, and low-resolution records do not resolve climate cycles within the mid-Pliocene. Our proxy data, therefore, represent comparisons of the change between average mid-Pliocene conditions and late Holocene or modern conditions. In contrast, PlioMIP2 modeling experiments are designed to target a particular orbital interval during the Pliocene (MIS KM5c) with present-day orbital forcings³⁷. Hence, the impact of orbital forcing is not included by PlioMIP2 experiments by design. The influence of orbital forcing on hydroclimate conditions is not recorded by most proxy records due to limited temporal resolution and dating accuracy of these records. Moreover, several high-resolution, continuous hydroclimate records from North Africa have shown relatively small hydroclimate variability during the late Pliocene⁷⁷⁻⁷⁹. Keeping in mind this limitation, our analyses focus on qualitatively assessing the agreement between proxies and models. Despite this limitation, our compilation shows similar large-scale hydroclimate features compared to previous Pliocene hydroclimate synthesis efforts, notably, wetter conditions evidenced by widespread lakes, reduced dust flux, and more mesic environments across southern and eastern Asia, North Africa, and parts of Europe. For proxy-model comparison, we include records spanning 0–67°N and 23°W–172°E. Details of broad regional trends are discussed in the SI.

**Statistical analysis.** Because the majority of these records have qualitative or semi-quantitative interpretations, we assess the fit between proxies and models qualitatively. We classify modeled mid-Pliocene $δ(P-E)$ as indicating wetter, drier, or no change relative to PI simulations from the same model, and we classify proxies as indicating wetter, drier, or no change based on the author's original interpretation of mid-Pliocene vs. present-day or late Quaternary conditions. To convert continuous, quantitative model outputs of $δ(P-E)$ into categorical data, we classify model output as showing wetter, drier, or no change at different thresholds of % change in P-E. For instance, choosing a 10% threshold, we would classify models that show more than a 10% increase in P-E as showing wetter conditions for a particular site. We vary the threshold at 1% intervals between 1 and 100% change and separately calculate proxy-model agreement.

The degree of proxy-model agreement was assessed using Gwet's AC statistic, a metric of inter-rater agreement similar to Cohen's kappa (κ). While the latter has been used in previous proxy-model comparison studies, the nature of the Pliocene proxy data, where the vast majority of records show wetter conditions, renders Cohen's κ statistic less robust⁸⁰. Cohen's κ is known to perform poorly as a result of Cohen's paradox whereby the statistic underestimates the true agreement between
<figure><img src="image_8.png" /><figcaption>Fig. 5 Relationship between earth system sensitivity (ESS, °C), equilibrium climate sensitivity (ECS, °C), and ability to reproduce mid-Pliocene proxy hydroclimate state by different models. Proxy data-model agreement measured by Gwet's ACs as a function of a global mean surface warming (δTs, which quantifies ESSs), and b ECSs of these models. δTs and ECSs are from ref. ³³. Boxplots show the spread of Gwet's AC (1st, 2nd, and 3rd quantiles, upper and lower fences, and outliers if any) from using 20-60% of δ(P-E) as thresholds (Method).</figcaption></figure>

raters in cases of skewed distributions of ratings across categories⁸⁰. We, therefore, use a related metric known as Gwet's AC, which is known to be resistant to this paradox. Gwet's AC statistic is similar to Cohen's in that the AC statistic is

$$AC = \frac{P_a - P_{exp}}{100\% - P_{exp}}$$

Where P_a is the actual percentage of agreement between proxies and models (e.g., the percentage of wetter sites that are correctly classified by the model as wetter, the percentage of drier sites correctly classified by the model, etc.), and P_exp reflects the expected percentage of agreement between raters (e.g., proxies and models) by chance alone. The Gwet's AC statistic is identical to Cohen's κ but uses a different formulation for P_exp that is not susceptible to Cohen's paradox. The statistical significance of Gwet's AC statistic was calculated according to the error estimator and methods outlined in ref. ⁸⁰. Results of statistical significance testing are presented in Fig. S3.

To avoid weighting our analysis towards regions with a greater density of proxy records less than 150 km that featured the same sign of the change (e.g., both showing wetter, drier, or no change) from other sites combined into one site. However, in cases with records showing opposite signs of change, we retain both records since in many cases there is not enough information to determine which record is more reliable, and excluding both would decrease the data coverage of our proxy compilation in key regions like East Asia (Fig. S2).

In the absence of a priori evidence that proxies reflect hydroclimate in a particular season, patterns of the proxy-model agreement are assessed using annually averaged P-E. We independently assess the fit between proxies and models for winter (December–February) and summer (June–September) separately (Fig. S4). The fit is much higher across all models and the multi-model mean for June–September rainfall. Proxy-model agreement for winter rainfall alone shows very low values of Gwet's AC statistic, suggesting that winter rainfall changes do not explain a significant component of the pattern seen in the proxy record. Furthermore, intermodel spread in Gwet's AC values on annually averaged rainfall show a strong correlation with intermodel Gwet's AC results calculated for June–September P-E (Fig. 2) at nearly all threshold values of P-E. This suggests that the summer seasonal signal, which is also the largest signal across model simulations (Fig. 1), drives the overall pattern of proxy-model agreement across models.

**PlioMIP2 ensemble.** We use a suite of model simulations conducted as part of the 2nd Pliocene Model Intercomparison Project (PlioMIP2)³⁵,⁸¹⁻⁸⁷. Boundary conditions are derived from the PRISM4 dataset³⁸. PRISM4 boundary conditions include information on land distributions, topography and bathymetry, vegetation, soils, lakes, and land ice cover. Two experimental protocols were developed, one implementing Pliocene conditions with a modern land/sea mask, and the enhanced experiment that included all boundary conditions³⁷,³⁸ pCO₂ is prescribed at 400 ppm. Other trace gases, orbital parameters, and solar output were set to be identical to the PI control simulation of each model.

Modeling groups were given the option to either prescribe vegetation changes or simulate vegetation changes using a dynamic global vegetation model. For the latter experiment, model simulations were started with PI vegetation and the model was allowed to spin up until a new equilibrium distribution of vegetation was achieved. We provide basic information on the configuration of each of the PlioMIP2 models used in our analyses in Supplementary Table 2. This information is also provided in the previous study³³. We note that only one modeling group in the suite of

simulations we analyzed opted to use a dynamic configuration for vegetation, suggesting that the choice of dynamic vs. static vegetation is not the primary source of spread across model results. We also note that models show varying levels of agreement with the signal in proxies (Fig. 2), suggesting that our results are not trivial (e.g., models with prescribed vegetation reproduce said vegetation) since the spatial pattern of hydroclimate in each model appears to depend on model design.

**Sensitivity experiments using community earth system model version 2 (CESM2).** Despite the overall similarity in geography and topography to present-day, mid-Pliocene boundary conditions feature several changes in ocean gateways, islands, and lake distributions that may influence regional climate³⁸,⁴¹,⁸⁸⁻⁹⁰. PlioMIP2 simulations also feature expanded grassland replacing the subtropical desert of North Africa and Central Asia and afforestation at northern high latitudes as well as deglaciated western Antarctic and most of the Greenland³⁸.

To isolate the mechanisms responsible for Pliocene hydroclimate changes, we carried out several new experiments with CESM2 that decompose simulated responses to the full mid-Pliocene climate forcing (F_all) into responses to forcings from elevated CO₂ (F_CO₂), changes in paleo-geography and -topography (F_geotop), and changes in biome distribution and ice sheets (F_vegice). We also constrain potential state dependency of responses induced by F_CO₂ and F_geotop. Because the forcing from F_vegice only emerges as a result of long-term CO₂ changes and a fixed F_geotop condition, we build this state dependency into the experimental approach. These new experiments separately feature (1) a 400 ppm CO₂ and PI geography and topography (E400); (2) PI vegetation and ice sheets, but otherwise mid-Pliocene CO₂, geography and topography (Eo400), (3) mid-Pliocene geography and topography, but PI CO₂, and vegetation and ice sheets (Eo280). The design and naming convention of these simulations roughly follow the Tier II PlioMIP2 protocol³⁷. All simulations were run with 0.9 × 1.6° resolution for atmosphere and land, and 1° nominal resolution for ocean and sea ice components, resulting in ~100 km resolution of all model components.

These new simulations are carried out for a minimum of 300 model years. Diagnostics of equilibrium by the global mean top of the atmosphere radiation imbalance and surface temperature are shown in Fig. S13. To produce climatology, we average the last 100 years of the model simulation. Eo400, Eo400, Eo280 together with the published full forcing (Eoi400) and PI experiments (E280) allow the decomposition of climate responses (denoted as R) to F_all into the sum of responses to individual forcings: R(F_all) = R(F_vegice) + R(F_geotop) + R(CO₂), for which R(F_all) = R(Eoi400) - R(E280); R(F_geotop) = R(Eo400) - R(E400) or R(F_geotop) = R(Eo280) - R(E280); R(F_CO₂) = R(Eo400) - R(E280) or R(F_CO₂) = R(Eo400) - R(Eo280); R(F_vegice) = R(Eoi400) - R(Eo400). The comparison of R(F_geotop), estimated with mid-Pliocene and PI levels of pCO₂, quantifies the dependency of R(F_geotop) on the background pCO₂. Similarly, the comparison of R(F_CO₂), estimated with mid-Pliocene and PI geography and topography, quantifies the dependency of R(F_CO₂) on the background geography and topography (Fig. S7).

**Development of moisture budget diagnostics.** To diagnose the causes of changes in terrestrial water balance, measured by changes in precipitation (P) minus evaporation (E) and referred to as δ(P-E), we further developed and applied the moisture budget diagnostics⁴⁷ to PlioMIP2 simulations. With only monthly data available, our derivation aims to facilitate the comparison of a pair of experiments and the evaluation of existing hypotheses regarding the cause of mid-Pliocene hydroclimate
change (e.g., identifying whether mid-Pliocene hydroclimate anomalies results from zonal mean changes or stationary wave changes).

Following Eq. (13) in ref. 47 on pressure coordinates, P-E is balanced by changes in the moisture tendency and moisture convergence

$$P - E = - \frac{1}{g\rho_w} \frac{\partial}{\partial t} \int_0^{P_s} qdp - \frac{1}{g\rho_w} \nabla \cdot \int_0^{P_s} uqdp \quad (1)$$

In Eq. (1), g is geopotential acceleration, $\rho_w$ is the density of water, $P_s$ is surface pressure, q is specific humidity, p is pressure, and u is horizontal wind. For a pair of experiments, experiment 1 is the control case (e.g., PI) and experiment 2 is the sensitivity experiment (e.g., Pliocene), the small perturbation method tells us that $q_2 = q_1 + \delta q$; $u_2 = u_1 + \delta u$; $(P - E)_2 = (P - E)_1 + \delta(P - E)$; and $P_{s2} = P_{s1} + \delta P_s$. The experiment (1 or 2) is identified by the subscript. $\delta$ terms are the small perturbations. We can decompose the perturbation to P-E into contributions from the anomalous moisture tendency, changes in convergence due to changes in wind, moisture, and the interaction of these changes, as well as a residual term, which are shown on the right side of Eq. (2) from the left to the right:

$$\delta(P - E) = - \frac{1}{g\rho_w} \frac{\partial}{\partial t} \int_0^{P_{s1}} \delta qdp - \frac{1}{g\rho_w} \nabla \cdot \int_0^{P_{s1}} (q_1 \delta u + u_1 \delta q + \delta u \delta q) dp + \text{res1} \quad (2)$$

$$\text{res1} = - \frac{1}{g\rho_w} \frac{\partial}{\partial t} q_2 \delta P_s - \frac{1}{g\rho_w} \nabla \cdot (u_2 q_2 \delta P_s) \quad (3)$$

Residual term 1 (RES 1, Eq. 3) quantifies influences from changing surface pressure, which is dominated by topographic changes between the mid-Pliocene and PI. Apply Reynold's decomposition, we can separate the total change of a given quantity into its temporal mean (e.g., monthly or annual climatology), denoted via an overbar ($\bar{\cdot}$), and higher frequency temporal fluctuations, denoted via a prime ($\cdot'$). Such that

$$\delta q = \delta \bar{q} + \delta q' \quad \delta u = \delta \bar{u} + \delta u' \quad (4)$$

$$\delta(P - E) = \bar{\delta}(P - E) + \delta(P - E)' \cdot u_1 = \bar{u}_1 + u_1' \cdot q_1 = \bar{q}_1 + q_1' \cdot P_{s1} = \bar{P}_{s1} + P_{s1}' \quad (5)$$

Given that we are interested in understanding the contributions from zonal mean circulation changes compared to other changes, we further separate $\delta u$ into changes in zonal mean, indicated by square brackets ([]) and changes in deviations from the zonal mean (i.e., the stationary wave), which is indicated by an asterisk (*)

$$\delta \bar{u} = [\delta \bar{u}] + \delta u^*$$

Using these expansions, we can rephrase the anomalous moisture budget in Eq. (2) to separate out the contributions due to zonal mean and stationary waves. The updated form of Eq. 2 is

$$\bar{\delta}(P - E) = - \frac{1}{g\rho_w} \frac{\partial}{\partial t} \int_0^{P_{s1}} \bar{\delta} qdp - \frac{1}{g\rho_w} \nabla \cdot \int_0^{P_{s1}} (\bar{u}_1 \delta \bar{q} + \bar{q}_1 [(\delta \bar{u}) + \delta \bar{u}^*]) + \delta \bar{u} \delta \bar{q}) dp + \text{res1} + \text{res2} \quad (6)$$

Residual term 2 (Eq. 7) quantifies combined effects of transient eddies moisture transport and influxes from the surface:

$$\text{res2} = - \frac{1}{g\rho_w} \frac{\partial}{\partial t} \overline{\delta q' P_s'} + \frac{1}{g\rho_w} \nabla \cdot \left[ [(\bar{u}_1 \delta \bar{q} + \bar{u}_1' \delta \bar{q} + q_1' \delta u + \bar{q}_1 \delta u' + \delta u' \delta \bar{q} + \delta u \delta \bar{q})] P_s \right] \\ - \frac{1}{g\rho_w} \nabla \cdot \int_0^{P_{s1}} [(\delta \bar{u} \delta \bar{u} + q_1' \delta u' + u_1' \delta q')] dp \quad (7)$$

Calculating residual terms 1 and 2 require high-frequency outputs of surface pressure, three-dimensional specific humidity, and horizontal winds, which are not available. Also, even at the 6-hourly resolution, the calculated eddy terms are insufficient to close the moisture budget with reanalysis data47. Thereby, RES 1 and RES 2 are not explicitly calculated in our calculations and quantified as a combined residual (i.e., the difference between P-E changes and the sum of other terms in the moisture budget equation). Yet, the decomposition demonstrates physical interpretations of these residuals. Monthly climatologies of surface pressure, precipitation, evaporation, three-dimensional horizontal winds, and specific humidity are used to calculate the remaining terms in Eq. 6. From the left to right, the first term in Eq. 6 ($-\frac{1}{g\rho_w} \frac{\partial}{\partial t} \int_0^{P_{s1}} \delta qdp$) is the moisture tendency term, which quantifies contributions to $\delta(P - E)$ from changes in seasonal cycle. The term $-\frac{1}{g\rho_w} \nabla \cdot \int_0^{P_{s1}} (\bar{u}_1 \delta \bar{q}) dp$ describes contributions from changes in climatological mean tropospheric moisture content; $-\frac{1}{g\rho_w} \nabla \cdot \int_0^{P_{s1}} (\bar{q}_1 [(\delta \bar{u})] dp$ describes

contributions from changes in the zonal mean circulation; $-\frac{1}{g\rho_w} \nabla \cdot \int_0^{P_{s1}} (\delta u \delta \bar{q}) dp$ describes contributions from changes in stationary wave kinetics.

Finally, $-\frac{1}{g\rho_w} \nabla \cdot \int_0^{P_{s1}} (\delta u \delta \bar{q}) dp$ describes contributions from covarying changes in mean moisture content and horizontal circulation.

## Data availability

The PlioMIP2 simulations analyzed in this study that are part of the Climate Model Intercomparison Project 6 (CMIP6) have been deposited to the Earth System Grid Federation with the access link: https://doi.org/10.22033/ESGF/CMIP6.4804 for EC-Earth3, https://doi.org/10.22033/ESGF/CMIP6.5230 for IPSL-CM6A-LR, https://doi.org/10.22033/ESGF/CMIP6.7227 for GISS-E2-1-G, https://doi.org/10.22033/ESGF/CMIP6.7675 for CESM2, and https://doi.org/10.22033/ESGF/CMIP6.12130 for HadGEM. The processed multi-model ensemble of PlioMIP2 and CESM2 data are available at Zenedo with the access link: https://zenodo.org/record/5706370#.YZbc6y1h0eY.

## Code availability

Source code of CESM2³⁹ can be downloaded from https://escomp.github.io/CESM/versions/cesm2.1/html/downloading_cesm.html.

The code for moisture budget analysis⁹¹ is published at Zenedo with the access link: https://zenodo.org/record/5706370#.YZbc6y1h0eY

Received: 3 September 2021; Accepted: 1 February 2022; Published online: 14 March 2022

## References

1. Hailemichael, M., Aronson, J. L., Savin, S., Tevesz, M. J. & Carter, J. G. $\delta$18O in mollusk shells from Pliocene Lake Hadar and modern Ethiopian lakes: Implications for history of the Ethiopian monsoon. *Palaeogeogr. Palaeoclimatol. Palaeoecol.* **186**, 81–99 (2002).

2. Schuster, M. et al. Chad Basin: paleoenvironments of the Sahara since the Late Miocene. *Comptes Rendus Geosci.* **341**, 603–611 (2009).

3. Sanyal, P., Bhattacharya, S. K., Kumar, R., Ghosh, S. K. & Sangode, S. J. Mio-Pliocene monsoonal record from Himalayan foreland basin (Indian Siwalik) and its relation to vegetational change. *Palaeogeogr. Palaeoclimatol. Palaeoecol.* **205**, 23–41 (2004).

4. Wang, H. et al. Asian monsoon rainfall variation during the Pliocene forced by global temperature change. *Nat. Commun.* **10**, 1–8 (2019).

5. Haywood, A. M., Dowsett, H. J. & Dolan, A. M. Integrating geological archives and climate models for the mid-Pliocene warm period. *Nat. Commun.* **7**, 1–14 (2016).

6. La Vega, De, E., Chalk, T. B., Wilson, P. A., Bysani, R. P. & Foster, G. L. Atmospheric CO 2 during the mid-Piacenzian warm period and the M2 glaciation. *Sci. Rep.* **10**, 1–8 (2020).

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

8. He, J. & Soden, B. J. A re-examination of the projected subtropical precipitation decline. *Nat. Clim. Change* **7**, 53–57 (2017).

9. Chou, C. & Neelin, J. D. Mechanisms limiting the northward extent of the northern summer monsoons over North America, Asia, and Africa. *J. Clim.* **16**, 406–425 (2003).

10. Bhattacharya, T., Tierney, J. E., Addison, J. A. & Murray, J. W. Ice-sheet modulation of deglacial North American monsoon intensification. *Nat. Geosci.* **11**, 848–852 (2018).

11. Scheff, J. & Frierson, D. M. Robust future precipitation declines in CMIP5 largely reflect the poleward expansion of model subtropical dry zones. *Geophys. Res. Lett.* **39**, L18704 https://doi.org/10.1029/2012GL052910 (2012).

12. Byrne, M. P. & O'Gorman, P. A. The response of precipitation minus evapotranspiration to climate warming: why the 'wet-get-wetter, dry-get-drier' seaping does not hold over land. *J. Clim.* **28**, 8078–8092 (2015).

13. Seager, R. & Vecchi, G. A. Greenhouse warming and the 21st century hydroclimate of southwestern North America. *Proc. Natl. Acad. Sci. USA* **107**, 21277–21282 (2010).

14. Swann, A. L., Fung, I. Y. & Chiang, J. C. Mid-latitude afforestation shifts general circulation and tropical precipitation. *Proc. Natl. Acad. Sci. USA* **109**, 712–716 (2012).

15. Berg, A. et al. Land-atmosphere feedbacks amplify aridity increase over land under global warming. *Nat. Clim. Change* **6**, 869–874 (2016).

16. Burke, K. D. et al. Pliocene and Eocene provide best analogs for near-future climates. *Proc. Natl. Acad. Sci. USA* **115**, 13288–13293 (2018).

17. Tierney, J. E. et al. Past climates inform our future. *Science* **370**, eay3701 (2020).

18. Burls, N. J. & Fedorov, A. V. Wetter subtropics in a warmer world: contrasting past and future hydrological cycles. *Proc. Natl. Acad. Sci. USA* **114**, 12888–12893 (2017).

19. Dowsett, H. J. & Robinson, M. M. Mid-Pliocene equatorial Pacific sea surface temperature reconstruction: a multi-proxy perspective. *Philos. Trans. R. Soc. A* **367**, 109–125 (2009).

20. Fedorov, A. V. et al. Patterns and mechanisms of early Pliocene warmth. *Nature* **496**, 43–49 (2013).

21. Zhang, Y. G., Pagani, M. & Liu, Z. A 12-million-year temperature history of the tropical pacific ocean. *Science* **344**, 84–87 (2014).
22. Ford, H. L., Ravelo, A. C. & Hovan, S. A deep Eastern Equatorial Pacific thermocline during the early Pliocene warm period. *Earth Planet. Sci. Lett.* **355**, 152–161 (2012).

23. Tierney, J. E., Haywood, A. M., Feng, R., Bhattacharya, T. & Otto-Bliesner, B. L. Pliocene warmth consistent with greenhouse gas forcing. *Geophys. Res. Lett.* **46**, 9136–9144 (2019).

24. Goldner, A., Herold, N. & Huber, M. The challenge of simulating the warmth of the mid-Miocene climatic optimum in CESM1. *Climate* **10**, 536 (2014).

25. Seager, R. et al. Mechanisms of ENSO-forcing of hemispherically symmetric precipitation variability. *Q. J. R. Met. Soc.* **131**, 1501–1527 (2005).

26. Brierley, C. M. et al. Greatly expanded tropical warm pool and weakened Hadley circulation in the early Pliocene. *Science* **323**, 1714–1718 (2009).

27. Fedorov, A. V., Burls, N. J., Lawrence, K. T. & Peterson, L. C. Tightly linked zonal and meridional sea surface temperature gradients over the past five million years. *Nat. Geosci.* **8**, 975–980 (2015).

28. Sun, Y., Ramstein, G., Contoux, C. & Zhou, T. A comparative study of large-scale atmospheric circulation in the context of a future scenario (RCP4.5) and past warmth (mid-Pliocene). *Climate* **9**, 1613–1627 (2013).

29. Rind, D. Latitudinal temperature gradients and climate change. *J. Geophys. Res.* **103**, 5943–5971 (1998).

30. Feng, R., Poulsen, C. J. & Werner, M. Tropical circulation intensification and tectonic extension recorded by Neogene terrestrial 8180 records of the western United States. *Geology* **44**, 971–974 (2016).

31. Zhu, J., Poulsen, C. J. & Tierney, J. E. Simulation of Eocene extreme warmth and high climate sensitivity through cloud feedbacks. *Sci. Adv.* **5**, eaax1874 (2019).

32. Feng, R. et al. Contributions of aerosol-cloud interactions to mid-Piacenzian seasonally ice-free Arctic Ocean. *Geophys. Res. Lett.* **46**, 9920–9929 (2019).

33. Haywood, A. M. et al. A return to large-scale features of Pliocene climate: the Pliocene Model Intercomparison Project Phase 2. *Clim. Past* **16**, 2095–2123 (2020).

34. de Nooijer, W. et al. Evaluation of Arctic warming in mid-Pliocene climate simulations. *Climate* **16**, 2325–2341 (2020).

35. Feng, R., Otto-Bliesner, B. L., Brady, E. C. & Rosenbloom, N. Increased Climate Response and Earth System Sensitivity from CCSM4 to CESM2 in mid-Pliocene simulations. *J. Adv. Model. Earth Syst.* https://doi.org/10.1029/2019MS002033 (2020).

36. Berntell, E. et al. Mid-Pliocene West African Monsoon Rainfall as simulated in the PlioMIP2 ensemble. *Clim. Past* **17**, 1777–1794 (2021).

37. Haywood, A. et al. The Pliocene Model Intercomparison project (PlioMIP) phase 2: scientific objectives and experimental design. *Climate* **12**, 663–675 (2016).

38. Dowsett, H. et al. The PRISM4 (mid-Piacenzian) paleoenvironmental reconstruction. *Climate* **12**, 1519–1538 (2016).

39. Danabasoglu, G. et al. The Community Earth System Model version 2 (eCSM2). *J. Adv. Model. Earth Syst.* **12**, e2019MS001916 (2020).

40. Salzmann, U., Haywood, A. M., Lunt, D. J., Valdes, P. J. & Hill, D. J. A new global biome reconstruction and data-model comparison for the Middle Pliocene. *Glob. Ecol. Biogeogr.* **17**, 432–447 (2008).

41. Pound, M. J., Tindall, J. & Pickering, S. J. Late Pliocene lakes and soils: a global data set for the analysis of climate feedbacks in a warmer world. *Clim. Past* **10**, 167–180 (2014).

42. Broccoli, A. J. & Manabe, S. The influence of continental ice, atmospheric CO₂, and land albedo on the climate of the last glacial maximum. *Clim. Dyn.* **1**, 87–99 (1987).

43. Rohling, E. J. et al. Making sense of palaeoclimate sensitivity. *Nature* **491**, 683–691 (2012).

44. Colman, R. & McAvaney, B. Climate feedbacks under a very broad range of forcing. *Geophys. Res. Lett.* **36**, L01702 https://doi.org/10.1029/2008GL03626 (2009).

45. Caballero, R. & Huber, M. State-dependent climate sensitivity in past warm climates and its implications for future climate projections. *Proc. Natl. Acad. Sci. USA* **110**, 14162–14167 (2013).

46. Swann, A. L., Fung, I. Y., Levis, S., Bonan, G. B. & Doney, S. C. Changes in Arctic vegetation amplify high-latitude warming through the greenhouse effect. *Proc. Natl. Acad. Sci. USA* **107**, 1295–1300 (2010).

47. Seager, R. & Henderson, N. Diagnostic computation of moisture budgets in the ERA-Interim reanalysis with reference to analysis of CMIP3-archived atmospheric model data. *J. Clim.* **26**, 7876–7901 (2013).

48. Peixoto, J. P. & Oort, A. H. Various decompositions of the circulation. *Phys. Clim.* **61**–64 (American Institute of Physics, New York, NY, 1992).

49. Trenberth, K. E., Stepaniak, D. P. & Caron, J. M. The global monsoon as seen through the divergent atmospheric circulation. *J. Clim.* **13**, 3969–3993 (2000).

50. Biasutti, M. Rainfall trends in the African Sahel: characteristics, processes, and causes. *Wiley Interdiscip. Rev. 10*, e591 (2019).

51. Coumou, D., Lehmann, J. & Beckmann, J. The weakening summer circulation in the Northern Hemisphere mid-latitudes. *Science* **348**, 324–327 (2015).

52. Dai, A. & Wigley, T. Global patterns of ENSO-induced precipitation. *Geophys. Res. Lett.* **27**, 1283–1286 (2000).

53. Skinner, C. B., Lora, J. M., Payne, A. E. & Poulsen, C. J. Atmospheric river changes shaped mid-latitude hydroclimate since the mid-Holocene. *Earth Planet. Sci. Lett.* **541**, 116293 (2020).

54. Shields, C. A., Kiehl, J. T., Rush, W., Rothstein, M. & Snyder, M. A. Atmospheric rivers in high-resolution simulations of the Paleocene Eocene Thermal Maximum (PETM). *Palaeogeogr. Palaeoclimatol. Palaeoecol.* **567**, 110293 (2021).

55. Menemenlis, S., Lora, J. M., Löfverström, M. & Chandan, D. Influence of stationary waves on mid-Pliocene atmospheric rivers and hydroclimate. *Glob. Planet. Change* **204**, 103557 (2021).

56. Broccoli, A. J., Dahl, K. A. & Stouffer, R. J. Response of the ITCZ to Northern Hemisphere cooling. *Geophys. Res. Lett.* **33**, L01702 https://doi.org/10.1029/2005GL02454 (2006).

57. Frierson, D. M. et al. Contribution of ocean overturning circulation to tropical rainfall peak in the Northern Hemisphere. *Nat. Geosci.* **6**, 940 (2013).

58. Schneider, T., Bischoff, T. & Haug, G. H. Migrations and dynamics of the intertropical convergence zone. *Nature* **513**, 45–53 (2014).

59. Biasutti, M. et al. Global energetics and local physics as drivers of past, present and future monsoons. *Nat. Geosci.* **11**, 392–400 (2018).

60. Tabor, C., Otto Bliesner, B. & Liu, Z. Speleothems of South American and Asian Monsoons influenced by a Green Sahara. *Geophys. Res. Lett.* https://doi.org/10.1029/2020GL089695 (2020).

61. Charney, J. G. Dynamics of deserts and drought in the Sahel. *Q. J. R. Met. Soc.* **101**, 193–202 (1975).

62. Claussen, M. et al. Simulation of an abrupt change in Saharan vegetation in the mid-Holocene. *Geophys. Res. Lett.* **26**, 2037–2040 (1999).

63. Rodwell, M. J. & Hoskins, B. J. Monsoons and the dynamics of deserts. *Q. J. R. Met. Soc.* **122**, 1385–1404 (1996).

64. Alkama, R., Kageyama, M. & Ramstein, G. A sensitivity study to global desertification in cold and warm climates: results from the IPSL OAGCM model. *Clim. Dyn.* **38**, 1629–1647 (2012).

65. Armstrong, E., Valdes, P., House, J. & Singarayer, J. Investigating the feedbacks between CO 2, vegetation and the AMOC in a coupled climate model. *Clim. Dyn.* **53**, 2485–2500 (2019).

66. Pausata, F. S., Messori, G. & Zhang, Q. Impacts of dust reduction on the northward expansion of the African monsoon during the Green Sahara period. *Earth Planet. Sci. Lett.* **434**, 298–307 (2016).

67. O'ishi, R. et al. PMIP4/CMIP6 last interglacial simulations using three different versions of MIROC: importance of vegetation. *Climate* **17**, 21–36 (2021).

68. Rosenfeld, D., Rudich, Y. & Lahav, R. Desert dust suppressing precipitation: a possible desertification feedback loop. *Proc. Natl. Acad. Sci. USA* **98**, 5975–5980 (2001).

69. Yoshioka, M. et al. Impact of desert dust radiative forcing on Sahel precipitation: relative importance of dust compared to sea surface temperature variations, vegetation changes, and greenhouse gas warming. *J. Clim.* **20**, 1445–1467 (2007).

70. Min, Q.-L. et al. Evidence of mineral dust altering cloud microphysics and precipitation. *Atmos. Chem. Phys.* **9**, 3223–3231 (2009).

71. Zhao, C., Liu, X., Ruby Leung, L. & Hagos, S. Radiative impact of mineral dust on monsoon precipitation variability over West Africa. *Atmos. Chem. Phys.* **11**, 1879–1893 (2011).

72. Pound, M. J., Haywood, A. M., Salzmann, U. & Riding, J. B. Global vegetation dynamics and latitudinal temperature gradients during the Mid to Late Miocene (15.97–5.33 Ma). *Earth Sci. Rev.* **112**, 1–22 (2012).

73. Utescher, T. & Mosbrugger, V. Eocene vegetation patterns reconstructed from plant diversity—a global perspective. *Palaeogeogr. Palaeoclimatol. Palaeoecol.* **247**, 243–271 (2007).

74. Caves, J. K. et al. The Neogene de-greening of Central Asia. *Geology* **44**, 887–890 (2016).

75. Lunt, D. J. et al. Earth system sensitivity inferred from Pliocene modelling and data. *Nat. Geosci.* **3**, 60–64 (2009).

76. Sun, J. Quantifying East Asian summer monsoon dynamics in the ECP4.5 scenario with reference to the mid-Piacenzian Warm Period. *Geophys. Res. Lett.* **45**, 12–523–12–533 (2018).

77. Feakins, S. J., Demenocal, P. B. & Eglinton, T. I. Biomarker records of late Neogene changes in northeast African vegetation. *Geology* **33**, 977–980 (2005).

78. Liddy, H. M., Feakins, S. J. & Tierney, J. E. Cooling and drying in northeast Africa across the Pliocene. *Earth Planet. Sci. Lett.* **449**, 430–438 (2016).

79. Demenocal, P. B. Plio-pleistocene African climate. *Science* **270**, 53–59 (1995).

80. G wet, K. L. Testing the difference of correlated agreement coefficients for statistical significance. *Educ. Psychol. Meas.* **76**, 609–637 (2016).

81. Hunter, S. J., Haywood, A. M., Dolan, A. M. & Tindall, J. C. The HadCM3 contribution to PlioMIP phase 2. *Climate* **15**, 1691–1713 (2019).

82. Stepanek, C., Samakinwa, E., Knorr, G. & Lohmann, G. Contribution of the coupled atmosphere-ocean-sea ice-vegetation model COSMOS to the PlioMIP2. *Climate* **16**, 2275–2323 (2020).

83. Zhang, Q. et al. Simulating the mid-Holocene, Last Interglacial and mid-Pliocene climate with EC-Earth3-LR. *Geosci. Model Dev.* **14**, 1147–1169 (2021).
84. Chan, W.-L. & Abe-Ouchi, A. Pliocene Model Intercomparison Project (PlioMIP2) simulations using the Model for Interdisciplinary Research on Climate (MIROC4m). *Clim. Past* **16**, 1523–1545 https://doi.org/10.5194/cp-16-1523-2020 (2020).

85. Tan, N. et al. Modeling a modern-like pCO 2 warm period (Marine Isotope Stage KM5c) with two versions of an Institut Pierre Simon Laplace atmosphere-ocean coupled general circulation model. *Climate* **16**, 1–16 (2020).

86. Li, X., Guo, C., Zhang, Z., Otterå, O. H. & Zhang, R. PlioMIP2 simulations with NorESM-L and NorESM1-F. *Climate* **16**, 183–197 (2020).

87. Williams, C. J. et al. Simulation of the mid-Pliocene Warm Period using HadGEM3: experimental design and results from model-model and model-data comparison. *Clim. Past* **17**, 2139–2163 https://doi.org/10.5194/cp-17-2139-2021 (2021).

88. Brierley, C. M. & Fedorov, A. V. Comparing the impacts of Miocene–Pliocene changes in inter-ocean gateways on climate: Central American Seaway, Bering Strait, and Indonesia. *Earth Planet. Sci. Lett.* **444**, 116–130 (2016).

89. Otto-Bliesner, B. L. et al. Amplified North Atlantic warming in the late Pliocene by changes in Arctic gateways. *Geophys. Res. Lett.* **44**, 957–964 (2017).

90. Hill, D. J. The non-analogue nature of Pliocene temperature gradients. *Earth Planet. Sci. Lett.* **425**, 232–241 (2015).

91. Feng, R. et al. Past terrestrial hydroclimate sensitivity controlled by Earth system feedbacks. *Zenodo* https://zenodo.org/record/5706370#.YZbc6y1h0eY (2021).

## Acknowledgements

The authors would like to thank all modeling groups who provided PMIP4 outputs for this analysis, WCRP, CMIP panel, PCMDI, ESGF infrastructures for sharing data, WCRP, and CLIVAR for supporting the PMIP project. R.F., T.B., B.L.O., and E.C.B acknowledge support from U.S. National Science Foundation grant numbers 1814029 and 1903650 (R.F.), 1903148 and 2103015 (T.B.) and 1852977 (B.L.O. and E.C.B.). X.L. and N.T. acknowledge support from the National Science Foundation of China grant numbers 42005042 (X.L.) and 41888101 (N.T.). D.L. acknowledges support from NERC (Natural Environment Research Council), SWEET Large Grant number NE/P01903X/1. C.C. acknowledges support from France ANR HADOc grant number ANR-17-CE31-0010. Q.Z. acknowledges support from Swedish Research Council (Vetenskapsrådet) grant numbers. 2013-06476 and 2017-04232. W.L.C. and A.A.O. acknowledge support from Japanese JSPS Kakenhi grant 17H06104 and NEXT Kakenhi grant 17H06323. H.D. acknowledges support from USGS Paleoclimate Research and Development Program. C.S. and G.L. acknowledge funding via the Alfred Wegener Institute's research programme PACES2. C.S. received funding via the Helmholtz Climate Initiative REKLIM. The PRISM4 reconstruction and boundary conditions used in the presented simulations were funded by the U.S. Geological Survey Climate and Land Use Change Research and Development Program. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government. The CESM2 simulations are performed with high-performance computing support from Cheyenne (doi:10.5065/D6RX99HX) provided by NCAR's Computational and Information Systems Laboratory, sponsored by the National Science Foundation. The IPSL-CM6A-LR simulation was run on the Très Grande Infrastructure de Calcul (TGCC) at Commissariat à l'Energie Atomique (gencpmip6 project) under the allocations 2016-A0030107732, 2017-R0040110492 and 2018-R0040110492 (project gencpmip6) provided by

GENCI (Grand Equipement National de Calcul Intensif). The model simulations with EC-Earth3 and the data analysis were performed using resources provided by ECMWF's computing and the Swedish National Infrastructure for Computing (SNIC) at the National Supercomputer Centre (NSC), which is partially funded by the Swedish Research Council through grant agreement no. 2018-05973. A.A.O. and W.L.C. acknowledge JAMSTEC for use of the Earth Simulator supercomputer. COSMOS simulations have been conducted at the Computing and Data Centre of the Alfred Wegener Institute - Helmholtz Centre for Polar and Marine Research on a NEC SX-ACE high-performance vector computer.

## Author contributions

R.F., T.B. B.L.O.-B., and A.M.H. conceptualized this study. R.F., T.B., M.K., and D.L. devised the methods for analysis. R.F. and T.B. led the investigation of the question, compiled data, and implemented the methods. R.F. and T.B. wrote the original draft. H.D. contributed data to create boundary conditions of PlioMIP2 simulations. E.C.B., A.H., J.T., S.H., A.A.-O., W.-L.C., M.K., C.C., C.G., X.L., G.L., C.S., N.T., Q.Z., Z.Z., Z.H., C.J.R.W., D.J.L., D.C., and W.R.P. contributed simulation data to this study, and reviewed and edited various versions of this paper.

## Competing interests

The authors declare no competing interests.

## Additional information

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41467-022-28814-7.

**Correspondence** and requests for materials should be addressed to Ran Feng.

**Peer review information** *Nature Communications* thanks Chris Brierley and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. Peer reviewer reports are available.

**Reprints and permission information** is available at http://www.nature.com/reprints

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

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

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2022