Past terrestrial hydroclimate driven by Earth System Feedbacks 4

Authors 5 R. Feng1*, T. Bhattacharya2, B. Otto-Bliesner3, E. Brady3, A. Haywood4, J. Tindall4, S. Hunter4, A. AbeOuchi5, W. -L. 6 Chan5, M. Kageyama6, C. Contoux6, C. Guo7, X. Li8, G. Lohmann9, C. Stepanek9, N. Tan10, Q. Zhang11, Z. Zhang8, Z. Han12, 7 J. R. Williams13, D. J. Lunt13, H. Dowsett14 8 9 Affiliations 10 1. Department of Geosciences, College of Liberal Arts and Sciences, University of Connecticut, Connecticut, USA 11 2. Department of Earth and Environmental Sciences, Syracuse University 12 3. Climate and Global Dynamics Laboratory, National Center for Atmospheric Research, Boulder, Colorado, USA 13 4. School of Earth and Environment, University of Leeds, Woodhouse Lane, Leeds, West Yorkshire, LS29JT, UK 14 5. Atmosphere and Ocean Research Institute, University of Tokyo, Kashiwa, Japan. 15 6. Laboratoire des Sciences du Climat et de l’Environnement, LSCE/IPSL, CEA-CNRS-UVSQ, Université Paris-Saclay, F16 91191 Gif-sur-Yvette, France 17 7. NORCE Norwegian Research Centre, Bjerknes Centre for Climate Research, 5007 Bergen, Norway 18 8. Department of Atmospheric Science, School of Environmental studies, China University of Geoscience, Wuhan, 19 430074, China 20 9. Alfred Wegener Institute-Helmholtz Centre for Polar and Marine Research, Bremerhaven, Germany 21 10. Key Laboratory of Cenozoic Geology and Environment, Institute of Geology and Geophysics, Chinese Academy of 22 Sciences, Beijing 100029, China 23 11. Department of Physical Geography and Bolin Centre for Climate Research, Stockholm University, Stockholm, 24 Sweden 25 12. College of Oceanography/Key Laboratory of Marine Hazards Forecasting, Ministry of Natural Resources/Key 26 Laboratory of Ministry of Education for Coastal Disaster and Protection, Hohai University, Nanjing, China 27 13. School of Geographical Sciences and Cabot Institute, University of Bristol, University Road, Bristol, BS8 1SS, UK 28 14. Florence Bascom Geoscience Center, U. S. Geological Survey, Reston, Virginia, USA 29 *Corresponding author. Email: ran.feng@uconn.edu 30 31 None-peer-reviewed preprint submitted to EarthArXiv 32 33 Currently in review 34 35 Submitted to Science Advances – 31 May 2021 36 37 Abstract 38 Geologic evidence suggests drastic reorganizations of subtropical terrestrial hydroclimate during past warm 39 intervals, including the mid-Piacenzian Warm Period (MP, 3.3 to 3.0 Ma). Despite having a similar to present-day 40 atmospheric CO2 level (pCO2), MP featured moist subtropical conditions with high lake levels in Northern Africa, 41 and mesic vegetation and sedimentary facies in subtropical Eurasia. Here, we demonstrate that major loss of the 42 northern high-latitude ice sheets and continental greening, not the pCO2 forcing, are key to generating moist 43 terrestrial conditions in subtropical Sahel and east Asia. In contrast to previous hypotheses, the moist conditions 44 simulated in both regions are a product of enhanced tropospheric humidity and a stationary wave response to the 45 surface warming pattern, both varying strongly in response to land cover changes. These results suggest that past 46 terrestrial hydroclimate states were driven by Earth System Feedbacks, which may outweigh the direct effect of 47 pCO2 forcing. 48 49 Teaser 50


Introduction 57
Geologic evidence suggests dramatic reorganizations of subtropical climate during past greenhouse climate 58 intervals, including the mid-Piacenzian Warm Period (3.3 to 3.0 Ma, MP for short). Multiple proxies of hydroclimate 59 indicate large, deep lakes and reduced dust flux across north Africa(1, 2). There is also evidence for more mesic 60 vegetation in South and East Asia(3, 4). Moreover, proxies generally suggest significant regions of subtropical Eurasia 61 were wetter prior to the onset of northern hemisphere glaciation. These changes, which imply large increases in 62 precipitation minus evaporation (P-E), are associated with pCO2 levels of approximately 400 ppm (5, 6). Moderate 63 increase of precipitation across these regions may be expected due to tropospheric moistening(7), enhanced land-64 sea thermal contrast(8), and enhanced inland moisture advection and ventilation(9, 10). However, evaporation also 65 increases in response to surface warming. As a result, predicted changes in terrestrial water balance (P-E) remain 66 equivocal across broad subtropical continents(11). 67 Model predicted response of terrestrial hydroclimate to CO2 increase broadly follows the "wet-gets-wetter, 68 dry-gets-drier" paradigm(7) with strong modulations from land warming pattern, land-sea warming contrast(12), 69 tropical SST pattern(13), and feedbacks from soil moisture and CO2 fertilization effects on leaf phenology(14, 15). The 70 net effect of these mechanisms potentially results in muted changes in P-E in subtropical Sahel and east Asia seen in 71 simulations featuring middle of the road future warming scenarios, such as the Shared Socioeconomic Pathways 72 (SSP245) (Fig. 1a and d). These future scenarios are often thought to be comparable to the MP climate(16, 17) ( Fig.  73 1a). The disparity between hydroclimate state recorded by MP proxies and future simulations challenges our 74 understanding of terrestrial hydroclimate sensitivity to CO2 radiative forcing. 75 Two leading hypotheses might explain the wetter subtropical continents during past warm climates. One 76 emphasizes the hydroclimate impact of an El Niño-like Pacific mean state. SST records from the tropical Pacific record 77 greater warming across the eastern equatorial Pacific than the western Pacific warm pool (18-22), resulting in an El-78 Niño-like pattern of SST anomalies. An El Niño-like tropical Pacific SST pattern may strengthen and shift the 79 subtropical jet equatorward, enhancing the transient eddy-driven moisture convergence and ascent(23, 24). The 80 other hypothesis focuses on the role of a sluggish Hadley Circulation (HC) that is inferred from a relaxed meridional 81 SST gradient(25-27). A weaker HC may result in weakened zonal mean moisture divergence from the subtropics and, 82 in turn, reduced aridity(28-30). 83 Here, using atmosphere-ocean coupled global climate model (GCM) simulations from the Pliocene Model 84 Intercomparison Project Phase II (PlioMIP2) (31-33), we assess whether the current generation ESMs can reproduce 85 the pattern of MP hydroclimate suggested by proxies. We further develop several new simulations using the 86 Community Earth System Model version 2(34, 35) to explore the extent to which simulated MP hydroclimate changes 87 across Sahel and subtropical Asia East Asia can be generated by changes in CO2 radiative forcing, tectonics, or 88 vegetation and ice sheets. In contrast to previous hypotheses, simulated Pliocene  P-E pattern in models and proxies 96 MP changes in precipitation minus evaporation (δ(P-E)) and other climate variables are calculated with 97 respect to preindustrial (PI) values. The last 100 years of simulations by 13 PlioMIP2 ESMs were averaged to produce 98 the ensemble mean (Table S1). A robust moistening signal that is larger than internal variability of individual models 99 is found across the Sahel and subtropical Eurasia (Fig. 1b). This pattern is the most pronounced during the boreal 100 summer months (June to September) ( Fig. 1c) with little or compensating changes during the winter months 101 (December to March) (Fig. S1). The spatial continuity of positive δ(P-E) from North Africa to subtropical eastern Asia 102 is not a visual coincidence: models that show a small precipitation increase in the Sahel also tend to show a small 103 increase in the southeast Asia (Fig. 1d), suggesting similar processes driving hydroclimate changes in both regions, 104 which is later confirmed by the moisture budget analysis. 105 To compare modeled patterns of hydroclimate change to available geologic data, we compiled proxy 106 indicators of MP terrestrial hydroclimate, drawing on existing compilations(30, 36, 37) as well as our own literature 107 search. We identify a total of 64 proxy records that include sedimentological indicators, palynological, floral or faunal, 108 offshore marine records, and stable isotope analyses of organic and inorganic materials (see Methods) (Fig. S2). These 109 records are most appropriately interpreted as qualitative or semi-quantitative indicators of hydroclimate. We 110 therefore quantify the extent to which proxies and models produce the same patterns of wetter, drier, or unchanged 111 MP hydroclimate changes using a metric known as Gwet's AC designed for categorical data (See Methods). To account 112 for different biases in climatological rainfall across models, model values of P-E are expressed as % changes from pre-113 industrial values. 114 The annual patterns of δ(P-E) revealed by the PlioMIP2 ensemble are consistent with proxy indicators of 115 hydroclimate across Sahel and subtropical Eurasia (Fig. 2), and are strongly driven by the summer δ(P-E) patterns ( Fig.  116  S3). Yet, different models show varying skill at capturing the pattern in the proxy records (Figs. 2a and S4). For  117  instance, IPSL-CM6, and EC-Earth3.3, two models with the highest level of agreement between MP proxies and  118 simulations, show expansive inland wetter conditions across North Africa, the Mediterranean, and subtropical Asia 119 ( Figure S4). In contrast, NorESM-L and GISS-E2-1-G, two of the models with the lowest proxy-model agreement, show 120 highly mixed or muted δ(P-E) across this region (Fig. S4). Despite intermodal spread in the degree of proxy-model 121 agreement, the majority of models, and the multi-model mean, show statistically significant proxy-model agreement 122 at δ(P-E) threshold of 20% above the PI (Fig. S4). While we analyze proxy-model agreement based on annually 123 averaged δ(P-E) in models, intermodal spread in the fit between proxy and model is strongly driven by the summer 124 signal (Fig. 2c) and where mid-latitude deserts becomes vegetated feature substantially lowered surface albedo (Fig. 3c). New 145 territories of boreal forests also show large increases of upward latent heat flux, suggesting enhanced water vapor 146 feedback(40) (Fig. 3b). Similar increase in latent heat flux also occurs in the northern subtropics where δ(P-E) is 147 positive, contributing to enhanced tropospheric humidity and warming.

149
Dynamical linkage between climate forcing conditions and MP hydroclimate 150 In order to identify the dynamical linkage between climate forcing conditions and hydroclimate in subtropical 151 Sahel and east Asia, we adapt the published moisture budget diagnostics (MBD) (41) to decompose simulated June 152 to September δ(P-E) by individual models (Methods) into changes in seasonal cycle of tropospheric moisture content 153 (δ(P-E)t), changes in zonal mean circulation dynamics (δ(P-E)[V]) and stationary wave dynamics (δ(P-E)V*) ( Fig. 4a and  154 4b), changes in tropospheric moisture content (δ(P-E)Q) (Fig. 4c), changes in interactions between moisture and 155 airflow dynamics δ(P-E)VQ, and a residual term (δ(P-E)Resi) (Fig. 4d). Stationary wave is quantified as the temporal 156 mean departure from the zonal mean following the classic circulation decomposition(42). The residual term combines 157 effects of high frequency variability of transient eddies and changes of topography (see Methods).

158
The MBD results are consistent among PlioMIP2 ensemble mean and Fvegice, suggesting common mechanisms 159 for δ(P-E) among these two sets of experiments ( Fig. 4e and f, and However, this is not the case. Positive δ(P-E) displays a clear latitudinal offset between land and ocean, and between 165 subtropical east Asia and North American continents.

166
A small δ(P-E)Resi suggests that the net influence from changes of transient eddies and topography is small 167 (Fig. 4d). In PlioMIP2 experiments, storm track activities and subtropical jet stream in the Northern Hemisphere both 168 weaken as shown by the 850 hPa eddy kinetic energy (EKE) and 200 hPa zonal wind speed ( Fig. S6a and b). These 169 responses are opposite to the expected responses to El Niño-like SSTs (24), and likely driven by other aspects of 170 simulated climate change such as Arctic amplification(43). Moreover, if the El Niño-like SST is key to driving the 171 ensemble mean δ(P-E), negative δ(P-E) in the northern Africa paired with positive δ(P-E) in the southeast Asia are 172 expected as part of the fingerprints of tropospheric teleconnection(44), which is not seen in the simulated MP δ(P-E) 173 (Fig. 1b).

174
Positive δ(P-E) in the subtropical Sahel and subtropical Asia primarily arises from changes in stationary wave 175 dynamics (δ(P-E)V*) (Fig. 2b) and increased zonal mean moisture content (δ(P-E)[Q]). Changes in stationary wave is 176 further decomposed into components associated with different wave numbers through Fourier decomposition of 177 stream function (SF) anomalies at 600 hPa. The SF anomalies are calculated as departures from the zonal mean and 178 reference state (Method) (Fig. 3d -f). As shown by the Fourier decomposition, wave number 1 (wave-1) of SF 179 anomalies accounts for 67% of the total wave energy between 0 -90°N. This wave component features cyclones 180 separately centered in the western Europe and North Africa, and anticyclones centered in the north Pacific ( Fig. 3d -181 f). Following the southern edge of the cyclonic wave centers, a moisture transport corridor emerges: westerly winds 182 bring moisture to North Africa from the tropical Atlantic; southerly and southwesterly winds bring moisture from 183 tropical Indian and Pacific Ocean towards Indian subcontinent and east Asia. Furthermore, the rotational winds of 184 wave-1 show a clear thermal wind relationship with the pattern of surface warming. This pattern of warming can 185 mostly be reproduced with Fvegice. 186 Positive δ(P-E) in both regions also results from increased tropospheric moisture content. Under PI 187 conditions, low level winds converge towards the north Africa and subtropical east Asia, with diverging flow across 188 the nearby ocean regions during the boreal summer (Fig. 3c). This regional circulation is a key feature of regional 189 summer monsoon(45). Even without changing circulation, elevated tropospheric moisture content can result in 190 greater moisture convergence and positive δ(P-E)Q in subtropical Sahel and east Asia, and greater moisture divergence 191 and negative δ(P-E)Q in the northern and southern side of the adjacent regions. This δ(P-E)Q pattern is a known 192 signature of thermodynamic response to elevated CO2, i.e., the wet-gets-wetter, dry-gets-drier paradigm(7, 12, 46). 193 As shown in Fvegice, similar thermodynamic response can be induced through land cover changes besides the CO2 194 change.

196
Discussion 197 Implications to hydrological cycle during Cenozoic warm intervals 198 The stationary wave response identified here is distinct from zonal mean Hadley circulation response or ITCZ 199 shift. These mechanisms highlight the importance of a changing zonal mean energy budget on circulation 200 dynamics(14, 30, 47-51). Instead, the stationary wave response in PlioMIP2 experiments is generated by spatially 201 heterogeneous warming pattern due to perturbations to regional radiation budget induced by continental greening. 202 This "pattern effect" may have been overlooked when examining past hydroclimate changes.

203
Why is Fvegice more effective at altering subtropical terrestrial hydroclimate compared to FCO2? Land cover 204 changes are known to be key to generating authigenic surface temperature and regional hydroclimate responses 205 across north Africa and subtropical east Asia (52-54) and may even alter the strength of Atlantic meridional 206 overturning circulation(55, 56). In North Africa, expansion of desert results in enhanced surface albedo, surface 207 cooling, and strengthened diabatic subsidence and dust emission(57, 58). These responses may further suppress 208 moist convection and perpetuate desert. Positive vegetation-precipitation feedback results in multiple equilibrium 209 states of North African vegetation over the late Quaternary(53). In our simulations, atmospheric circulation responses 210 to Fvegice facilitate moisture transport towards subtropical Sahel and east Asia via stationary wave response. This 211 response is closely tied to the warming pattern generated by Fvegice, which does not occur in FCO2 (Fig. S7).

212
Different hydroclimate responses to FCO2 and Fvegice also reflect the difference between future and MP land 213 surface processes. Increasing CO2 favors reduction of soil moisture and partitioning of surface heat flux towards 214 sensible heat, leading to enhanced surface warming. This in turn lowers relative humidity above the surface, and 215 diminishes continental cloud cover, contributing to negative δ(P-E) on land. A similar response of surface heat flux 216 partitioning and moisture deficit can be caused by the reduction of leaf transpiration in response to CO2 fertilization. 217 Both soil moisture feedback and CO2 fertilization drive predicted future subtropical terrestrial hydroclimate change 218 to a large extent(15). These changes may also contribute to muted δ(P-E) in FCO2 despite a modest increase in 219 precipitation. In contrast, Fvegice features no physiological effect of CO2, and produces large increases in latent over 220 sensible heat flux across continents (Fig. S8), creating favorable conditions for a more humid troposphere. As such, 221 moist subtropical Sahel and east Asia reflect a synergy of dynamic and thermodynamic responses to Fvegice. Reduced 222 ice sheet cover, expansive northern high-latitude boreal forests, and vegetated Sahel and central Asia have been 223 recorded during other Cenozoic warm intervals(59-61). These land cover changes were likely instrumental in driving 224 changes in global mean temperature and terrestrial hydrological cycle throughout Cenozoic.

Past hydroclimate states driven by Earth System Feedbacks 227
Our results suggest an alternate view of the role of vegetation changes in modulating continental 228 hydroclimate. Changes in regional circulation in the form of stationary wave responses lead to strengthened low-level 229 winds that import moisture into the subtropical Sahel and east Asia from tropical Atlantic and Indian Ocean. This 230 process is amplified by enhanced tropospheric moistening due to Fvegice, and distinct from previous mechanisms that 231 highlight the role of the zonal mean HC and authigenic terrestrial responses to land surface changes. Simulated 232 responses of stationary wave and tropospheric humidity clearly dominate δ(P-E) across subtropical Sahel and east 233 Asia among models (Fig. S9). Although full feedbacks from vegetation and ice sheets are not prognostically simulated 234 in most MP simulations (Table S1), these simulations highlight that the radiative perturbations from proxy constrained 235 changes of vegetation and ice sheets are key to generating hydroclimate changes in the Northern subtropics. 236 Moreover, these hydroclimate changes have minimal dependency on the prescribed MP topography and land-sea 237 distribution (Fig. 1g).

238
A key implication of our findings is that MP continental hydroclimate is more appropriately viewed as part 239 of the Earth System feedbacks instead of an immediate response to FCO2. This inference is supported by the strong 240 dependency between Gwet's AC and simulated global mean warming relative to preindustrial among models. The 241 latter reflects model diversity in Earth System Sensitivity, given that all simulations were run with the same CO2, and 242 similar MP vegetation and ice sheet conditions. In contrast, intermodal spread in Gwet's AC shows little dependency 243 on equilibrium climate sensitivity (ECS) among models (ECSs are reported in Ref (33)). 244 Moreover, the proxy data-model convergence suggests that MP terrestrial hydroclimate records do not pose 245 a fundamental challenge to our understanding of the physics of hydroclimate change across subtropical Sahel and 246 eastern Asia. As such, these results offer a resolution to the apparent difference between the projections of 247 hydroclimate changes in the Sahel and subtropical east Asia in the middle-of-the-road scenarios and strong geologic 248 evidence for moist Pliocene climate at similar levels of CO2: while the former is dominated by short-term response to 249 CO2 radiative forcing and internal variability, Pliocene hydroclimate reflects long-term adjustments of the Many records are discontinuous, and cannot provide quantitative comparisons between MP climate and late 282 Quaternary or pre-industrial conditions, and low-resolution records do not resolve climate cycles within the MP. 283 Our proxy data therefore represent comparisons of the change between average MP conditions and late Holocene 284 or modern conditions. In contrast, PlioMIP2 modeling experiments are designed to target a particular orbital 285 interval during the Pliocene (MIS KM5c), and are explicitly compared to pre-industrial control simulations(31). A key 286 assumption behind our proxy model comparison is therefore that the magnitude of Pliocene hydroclimate change is 287 greater than the impact of orbital variations. There is some evidence to support this assumption, since several high-288 resolution, continuous hydroclimate records suggest substantial shift in hydroclimate states from early and mid-289 Pliocene to early-Pleistocene. Documented hydroclimate shifts are well beyond the recorded shorter time-scale 290 hydroclimate variability(65-67). However, this assumption requires further testing with additional, high resolution, 291 well-dated hydroclimate records that can be directly compared to model outputs. Keeping in mind this source of 292 uncertainty, our analyses focus on qualitatively assessing the agreement between proxies and models. Despite 293 these considerations, our compilation shows similar large-scale hydroclimate features compared to previous 294 Pliocene hydroclimate synthesis efforts, notably, wetter conditions evidenced by widespread lakes, reduced dust 295 flux, and more mesic environments across southern and eastern Asia, north Africa, and parts of Europe. For proxy-296 model comparison, we include records spanning 0-67°N and 23°W-172°E. Details of broad regional trends are 297 discussed in the SI. 298 Statistical Analysis 299 Because the majority of these records have qualitative or semi-quantitative interpretations, we assess the fit 300 between proxies and models qualitatively. We classify modeled MP δ(P-E) as indicating wetter, drier, or no change 301 relative to pre-industrial simulations from the same model, and we classify proxies as indicating wetter, drier, or no 302 change based on the author's original interpretation of MP vs. present-day or late Quaternary conditions. To 303 convert continuous, quantitative model outputs of δ(P-E) into categorical data, we classify model output as showing 304 wetter, drier, or no change at different thresholds of % change in P-E. For instance, choosing a 10% threshold, we 305 would classify models that show more than a 10% increase in P-E as showing wetter conditions for a particular site.

306
We vary the threshold at 1% intervals between 1 and 100% change and separately calculate proxy-model 307 agreement. 308 The degree of proxy-model agreement was assessed using Gwet's AC statistic, a metric of inter-rater 309 agreement similar to Cohen's kappa (κ). While the latter has been used in previous proxy-model comparison 310 studies, the nature of the Pliocene proxy data, where the vast majority of records show wetter conditions, renders 311 Cohen's κ statistic less robust(68). Cohen's κ is known to perform poorly as a result of `Cohen's paradox' whereby 312 the statistic underestimates the true agreement between raters in cases of skewed distributions of ratings across 313 categories(68). We therefore use a related metric known as Gwet's AC, which is known to be resistant to this 314 paradox. The Gwet's AC statistic is similar to Cohen's in that the AC statistic is: 315 316 Where Pa is the actual percentage of agreement between proxies and models (e.g. the percentage of wetter sites 319 that are correctly classified by the model as wetter, the percentage of drier sites correctly classified by the model, 320 etc.), and Pexp reflects the expected percentage of agreement between raters (e.g. proxies and models) by chance 321 alone. The Gwet's AC statistic is identical to Cohen's κ but uses a different formulation for Pexp that is not susceptible 322 to Cohen's paradox. Statistical significance of the Gwet's AC statistic was calculated according to the error estimator 323 and methods outlined in ref (68). Results of statistical significance testing are presented in Figure S3.

324
To avoid weighting our analysis towards regions with a greater density of proxy records, records less than 325 150 km that featured the same sign of change (e.g. both showing wetter, drier or no change) from each other were 326 combined into one site. However, in cases with records showing opposite signs of change, we retain both records 327 since in many cases there is not enough information to determine which record is more reliable, and excluding both 328 would decrease the data coverage of our proxy compilation in key regions like East Asia ( Figure S2). 329 In the absence of a priori evidence that proxies reflect hydroclimate in a particular season, patterns of 330 proxy-model agreement are assessed using annually averaged P-E. We independently assess the fit between proxies 331 and models for winter (December -February) and summer (July-September) separately ( Figure S4). The fit is much 332 higher across all models and the multi-model mean for July-September rainfall. Proxy-model agreement for winter 333 rainfall alone show very low values of Gwet's AC statistic, suggesting that winter rainfall changes do not explain a 334 significant component of the pattern seen in proxy record. Furthermore, intermodel spread in Gwet's AC values on 335 annually averaged rainfall show a strong correlation with intermodal Gwet's AC results calculated for July-336 September P-E (Fig. 2) at nearly all threshold values of P-E. This suggests that the summer seasonal signal, which is 337 also the largest signal across model simulations (Fig. 1) trace gases, orbital parameters, and solar output were set to be identical to the pre-industrial control simulation of 347 each model. 348 Modeling groups were given the option to either prescribe vegetation changes or simulate vegetation 349 changes using a dynamic global vegetation model. For the latter experiment, model simulations were started with 350 pre-industrial vegetation and the model was allowed to spin up until a new equilibrium distribution of vegetation 351 was achieved. We provide basic information on configuration of each of the PlioMIP2 models used in our analyses 352 in Table S1, and is provided in the previous study(33). We note that only one modeling group in the suite of 353 simulations we analyzed opted to use a dynamic configuration for vegetation, suggesting that the choice of dynamic 354 vs. static vegetation is not the primary source of spread across model results. We also note that models show 355 varying levels of agreement with the signal in proxies (Fig. 2), suggesting that our results are not trivial (e.g. models 356 with prescribed vegetation reproduce said vegetation) since the spatial pattern of hydroclimate in each model 357 appear to depend on model design.

359
Sensitivity experiments using Community Earth System Model version 2 (CESM2) 360 Despite the overall similarity in geography and topography to present-day, MP boundary conditions feature 361 several changes in ocean gateways, islands, and lake distributions that may influence regional climate(32, 37, 77-79). 362 PlioMIP2 simulations also feature expanded grassland replacing the subtropical desert of North Africa and central 363 Asia and afforestation at northern high latitudes as well as deglaciated western Antarctic and most of the 364 Greenland(32). 365 To isolate the mechanisms responsible for Pliocene hydroclimate changes, we carried out two new 366 experiments with CESM2 that decompose simulated responses to the full MP climate forcing (Fall) into responses to 367 forcings from elevated CO2 (FCO2), changes in paleo-geography and -topography (Fgeotop), and changes in biome 368 distribution and ice sheets (Fvegice). These new experiments separately feature 1) a 400 ppm CO2 and preindustrial 369 boundary conditions (E400); 2) preindustrial vegetation and ice sheets, but otherwise MP CO2 and boundary 370 conditions (Eo400). The design and naming convention of these simulations follow the Tier II PlioMIP2 protocol(31). 371 All simulations are run with 0.9x1.6° resolution for atmosphere and land, and 1° nominal resolution for ocean and ice In Equation (1), is geopotential acceleration, , is the density of water, -is surface pressure, and isspecific 391 humidity, : pressure. For a pair of experiments, experiment 1 is the control case (e.g. preindustrial) and experiment 392 2 is the sensitivity experiment (e.g. Pliocene), the small perturbation method tells us that . = % + ; . = % + 393 ; ( − ) . = ( − ) % + ( − ); -. = -% + -. We can therefore decompose the perturbations to P-E 394 into contributions from the anomalous moisture tendency, changes in convergence due to changes in wind, moisture, 395 and the interaction of these changes, as well as a residual term: Applying Reynold's decomposition, we can separate the total change in a given quantity into its temporal 401 mean (e.g. monthly or annual climatology), denoted via an overbar, and higher frequency temporal fluctuations, 402 denoted via a prime. Such that: Using these expansions, we can rephrase the anomalous moisture budget in equation (2)  Calculating residual terms 1 and 2 require high frequency outputs of surface pressure, three-dimensional specific 426 humidity and horizontal winds, which are not available in PlioMIP2 database. Also, even at 6-hourly resolution, the 427 calculated eddy terms are insufficient to close the moisture budget with reanalysis observational data(41). Thereby, 428 resi 1 and resi 2 are not explicitly calculated in our calculations and quantified as a combined residual (i.e., the 429 difference between P-E changes and the sum of other terms in the moisture budget equation).

431
Monthly climatologies of surface pressure, precipitation, evaporation, and three-dimensional horizontal winds, and 432 specific humidity are used to calculate the remaining terms in equation 6. From the left to right, the first term in 433 Acknowledgments 441 The authors would like to thank all the modelling groups who provided the PMIP3 and PMIP4 outputs for 442 this analysis, WCRP, CMIP panel, PCMDI, ESGF infrastructures for sharing data, WCRP and CLIVAR for 443 supporting the PMIP project. We would also like to thank H. Dowsett for developing PRISM4D datasets.      Proxy recorded Pliocene regional hydroclimate patterns 622 We compiled Pliocene hydroclimate indicators from published records (Table S1). We rely on the of high salinity at 3 Ma towards the end of the MP, which could reflect changes to the regional 635 water budget, but has been interpreted to reflect higher water levels in the Black Sea that slightly        Changes of surface temperature (color shaded), wave number 1 (contour), and rotational winds (vectors) of stationary 717 wave in response to FCO2 simulated by CESM2. Notice that the stationary wave pattern has little dependency on 718 surface temperature changes, but mainly reflects the high-pressure system above Tibet developed during the boreal 719 summer.