Mid-Holocene climate of the Tibetan Plateau and hydroclimate in three major river basins based on high-resolution regional climate simulations

. The Tibetan Plateau (TP) contains the headwaters of major Asian rivers that sustain billions of people and plays an important role in both regional and global climate through thermal and mechanical forcings. Understanding the characteristics and changes to the hydrological regimes on the TP during the mid-Holocene (MH) will help understand the 10 expected future changes. Here, an analysis of the hydroclimates over the headwater regions of three major rivers originating in the TP, namely the Yellow, Yangtze and Brahmaputra rivers is presented, using dynamically downscaled climate simulations constructed using the Weather Research and Forecasting Model (WRF) coupled to the hydrological model WRF-Hydro. Green Sahara (GS) boundary conditions have also been incorporated into the global model so as to capture the remote feedbacks between the Saharan vegetation and the river hydrographs over the TP. Model-data comparisons show that 15 the dynamically downscaled simulations significantly improve the regional climate simulations over the TP in both the modern day and the MH, highlighting the crucial role of downscaling in both present-day and past climates. TP precipitation is also strongly affected by the greening of the Sahara, with a particularly large increase over the southern TP, as well as a delay in the monsoon withdrawal. The simulation results were first validated over the upper basins of the three rivers before the hydrological responses to the MH forcing for the three basins were quantified. Both the upper Yellow and Yangtze rivers 20 exhibit a decline in streamflow during the MH, especially in summer, which is a combined effect of less snowmelt and stronger evapotranspiration. The GS forcing caused a rise in temperature during the MH, as well as larger rainfall but less snowfall and greater evaporative water losses. The Brahmaputra River runoff is simulated to increase in the MH, due to greater net precipitation.


Introduction 25
The Tibetan Plateau (TP), often referred to as the Third Pole, with an average elevation exceeding 4000 m and an area of 2.5 × 10 6 km 2 , constitute the highest and largest plateau in the world. The TP also has the most extensive cryosphere outside the Arctic/Greenland and the Antarctic regions (Yao et al., 2012), including snow cover, glaciers and permafrost. Moreover, the TP is also widely believed to have a huge impact on both the local and global climate through thermal and mechanical forcing (Kutzbach et al., 1993;An et al., 2001;Molnar et al., 2010). Climate over the TP, featuring warm wet summers and 30 cold dry winters, is extremely sensitive to global warming and has undergone rapid and profound changes in recent years (Lehnert et al., 2016;Yao et al., 2019). With an elevation gradually decreasing from the northwest to the southeast, the TP also acts as the "water tower" and the source region of large rivers in Asia that flow down to almost half of humanity, including the Yellow, Yangtze and Brahmaputra (Yarlung Tsangpo) Rivers (Immerzeel et al., 2020).
The Yellow River, which is often referred to as the cradle of the Chinese civilization, rises in the northern part of the TP (Fig.  35 1) and supports nearly one third of the China's population (Huang et al., 2009). It passes across the Loess Plateau and the North China Plain, before emptying into the ocean on the east coast of China. The Yellow River's headwaters are located in the Bayan Har Mountains in the northern TP. Its water resources are vital to the unique but fragile ecosystems in the upper Yellow River Basin (UYEB) and are extremely sensitive to climate change. The longest river in Asia, the Yangtze, originates in the Tanggula Mountains over the eastern part of the TP and flows eastwards to the East China Sea (Fig. 1). Its 40 drainage basin comprises one-fifth of China's land area, which is home to nearly a third of the country's population. The Brahmaputra River, which originates in the southern TP and flows eastwards before turning southward into India in Arunachal Pradesh, is an important trans-boundary river. After entering India, it turns southwestward and eventually enters Bangladesh where it merges with the Padama River in the Ganges Delta. There is dramatic climatic heterogeneity within its drainage area, from the cold dry upper region with an altitude higher than 5000 m, which features glaciers and permanent 45 snow, to the humid and hot downstream area that is only 200 m above sea level. The Brahmaputra River flows across four Asian countries including countries with the two largest populations in the world, China and India. Understanding hydrological processes of these rivers is necessary for water resource management and are of significant economic importance for these countries.
Observations over the TP suggest a warming trend that is twice as large as the global average (IPCC, 2013), with spatial 50 differences in precipitation and glacier changes (Yao et al., 2012). An increasing trend in precipitation has been observed to be occurring in the western and central TP and a smaller precipitation increase has been detected in the eastern TP (Zhou et al., 2019). The TP has generally become wetter and warmer during the past three decades, although the southern part has turned drier (Yang et al., 2014). As an important hydrologic process in mountain regions, the melting of glaciers can greatly modify the characteristics of streamflow, including magnitude, timing, and variability (Kaser et al., 2010;Immerzeel et al., 55 2020). Both temperature and precipitation changes will significantly influence glacier melt (Barnett et al., 2005;Immerzeel et al., 2020), inducing heterogenous responses over different parts of the TP. Higher temperature can significantly enhance glacier melt by removing fresh snow and hence decreasing surface albedo (Barnett et al., 2005). Glacier melt is a particularly important source of streamflow during warm and dry seasons, which will enhance the downstream discharge of the rivers (Kaser et al., 2010;Yao et al., 2004). However, this increase is transitory and is expected to cease when the glaciers have result in a drop in streamflow in summer (Immerzeel et al., 2020) and pose severe threat to access to water resources in the TP and downstream areas (Kehrwald et al., 2008).
Rivers that originate in the TP often suffer from a scarcity of data necessary for hydrological simulation because of the extreme conditions under which measurements must be acquired. Regional hydrological simulation in these ungauged areas 65 therefore faces significant technical challenges due to the complex regional topography. The resolution of most global climate models (GCMs) is insufficient for this region, because regional climate over the TP is strongly affected by the steep topographic gradients that are poorly resolved in global models. This hinders our ability to characterize the water cycle, study regional climate change impacts and find solutions to local water problems. Therefore, this study has dynamically downscaled global climate simulations to better resolve important regional topographic features. 70 The study of past climate conditions is an important tool in any effort to understand the mechanisms influencing present and future climate and its variability. In addition, understanding the range of variability will help to improve climate models that are used to simulate the regional response to a variety of abrupt and long-term changes in temperature, atmospheric circulation, and moisture sources, which in turn regulate the river discharge flowing out of the TP. The mid-Holocene (MH, approximately 6000 years ago) is one of the most widely studied periods of the Quaternary and is a period of focus for the 75 Paleoclimate Modelling Intercomparison Project (PMIP). In the MH, variation in the Earth's orbital parameters led to differences in the latitudinal and seasonal distribution of incoming solar radiation (Berger 1978), resulting in a climate that was remarkably different from that of the present day. The insolation was anomalously strong (weak) in summer (winter) over the Northern Hemisphere (NH), leading to an enhancement of the seasonal cycle. A temperature optimum for the TP during the MH has been indicated by pollen data (Herzschuh et al., 2006;Lu et al., 2011;Ma et al., 2014). Understanding 80 MH climatic changes over the TP will help understand the impact of the recent/modern global warming and predict future climate change (Wanner et al., 2008). During the MH, South and East Asian monsoons are expected to have been enhanced due to the increased boreal summer insolation (Huo et al., 2021;Maher, 2008;Herzschuh, 2006). However, compared with the reconstructions, regional biases and large spreads are found in numerical simulations of the MH climate by the PMIP models, in particular, over and around the TP (Zheng et al., 2013). 85 During the MH, changes in insolation and albedo feedbacks also strengthened the African monsoon, making northern Africa much wetter than today (Skinner and Poulsen, 2016). Lakes, rivers, and wetlands formed and developed in the now arid regions of northern Africa and led to a northward extension of the vegetation cover over the Sahara, i.e., a "Green Sahara" (GS; Chandan and Peltier, 2020;Pausata et al., 2020). Studies have shown that a vegetated Sahara and the strengthened African monsoon can enhanced the Asian monsoon by altering the Walker circulation through changes in tropical Atlantic 90 SSTs (Pausata et al., 2017). The TP, which is situated within the region influenced by the Asian monsoon (Conroy and Overpeck, 2011) is thus also affected by the greening of the Sahara. One of the key goals of the present study is to outline the remote response of the MH TP climate to the GS teleconnection by analysing a set of simulations in which the land cover is changed from desert to shrubland over a large part of northern Africa (Chandan and Peltier, 2020). Since there are few studies that characterize the runoff over the TP under the MH climate conditions, in this study we also analyse the 95 hydrological response of the TP rivers to the MH climate forcing.
The objective of this article is thus to first evaluate the performance of the coupled UofT-CCSM4-WRF-Hydro model over the TP and then analyse the response of hydroclimate over the three rivers over the TP to MH orbital forcing and remote GS feedback, using high-resolution regional climate simulations. Section 2 describes experimental design and the model configurations. The validation of the results is discussed in Sect. 3, before the general patterns of MH anomalies are outlined. 100 The central point of this study is hydroclimatological anomalies in the upper Yellow, Yangtze and Brahmaputra River basins (UYEB, UYAB and UBB, respectively), which is presented in Sect. 4. Section 5 summarizes the main results of the present study.

Experimental design 105
This study employs dynamical downscaling, a modelling methodology based upon the use of a dynamically consistent, physically based high-resolution regional climate model (RCM) to downscale global climate simulations so as to better account for regional climate feedbacks. The global climate simulations used to drive the RCM have been obtained using the University of Toronto version of the National Center for Atmospheric Research (NCAR) Community Climate System Model version 4 (UofT-CCSM4; Peltier and Vettoretti, 2014;Chandan and Peltier, 2018). The global climate simulations are then 110 dynamically downscaled to 10-km resolution using the Weather Research and Forecasting (WRF) model version 4.1 with four different convection parameterization schemes, constituting a small physics ensemble. We have also coupled WRF-Hydro (Gochis et al., 2020) to the WRF model as a hydrological extension package for hydrometeorological (e.g., rainfall, runoff, groundwater flow, streamflow) simulations, which will also be a primary focus of the analyses to be discussed herein.
The synchronously coupled coupled WRF-Hydro modelling system outperforms the WRF-only model and thus has been 115 employed in atmosphere-hydrology simulations over various domains around the globe (Senatore et al., 2015;Kerandi et al., 2017;Somos-Valenzuela and Palmer, 2018). Our hydroclimatic analysis is focused on the upper Yellow, Yangtze and Brahmaputra River basins (UYEB, UYAB and UBB, respectively), the headwaters of which are all located on the TP.
Four global climate simulations were performed using UofT-CCSM4 at approximately 1• resolution, one for the historical period , one for the preindustrial (PI) reference period and two for the MH. UofT-CCSM4 is based on the standard 120 CCSM4 (Gent et al., 2011), but specific modifications have been made to the diapycnal diffusivity of the ocean component to make it more appropriate for paleoclimate simulations (Peltier and Vettoretti, 2014;Chandan and Peltier, 2017). These changes have been taken into account in both the MH and PI simulations to avoid introducing ambiguity related to different diapycnal mixing schemes when these two simulations are compared. Huo et al. (2021) has employed the same PI and MH global simulations for the purpose of studying MH monsoons in South and Southeast Asia. The two MH simulations differ 125 from each other in their land surface boundary conditions over northern Africa (Chandan and Peltier, 2020) but otherwise share the same orbital parameters and trace gases as the PMIP4 recommendations (Otto-Bliesner et al., 2017): MHREF is the "control MH" which uses the same land surface as in PI, while MHGS extends MHREF by adding the GS vegetation.
Compared to the PI, the MHREF experiments are forced by precessionally enhanced boreal summer insolation and slightly lower greenhouse-gas concentrations (Otto-Bliesner et al., 2017). The prescribed vegetation patterns used in the MHGS 130 experiment are similar to those used by Pausata et al. (2016) and are consistent with the proposed MH vegetation sensitivity experiments in Otto-Bliesner et al. (2017), which involves a substitution of the Sahara desert with evergreen shrubs in the south and steppe and savanna in the north. Moreover, based on the results of Hély et al. (2014), the African Guineo-Congolian rain forest has also been extended slightly northward (Chandan and Peltier, 2020).
Each UofT-CCSM4 simulation was separately downscaled over the TP for 15 years, using the WRF model with two nested 135 domains. The outer domain covers most of Asia at 30-km resolution, and the inner domain covers the TP, as well as parts of the surrounding territory, at 10-km resolution. The outlines of the two domains as well as the topography in the WRF outer domain are shown in Fig. 1. This study is based on a four-member physics ensemble of dynamically downscaled climate simulations, using the same physics parameterizations as in Huo et al. (2021), where the MH climate anomalies over the South and Southeast Asia were studied. The simulations in each regional climate ensemble employ different cumulus 140 parameterizations in WRF as in Huo et al. (2021), namely the Kain-Fritsch scheme (KF; Kain, 2004), the Grell-Freitas ensemble scheme (GF; Grell and Freitas, 2014), the Tiedtke scheme (Tiedtke, 1989) and the Betts-Miller-Janjić scheme (BMJ; Janjić, 1994), making a mini-physics ensemble for the TP, that allows us to study the sensitivity of model performance to different cumulus parameterizations and thereby to estimate the uncertainty associated with these parameterizations on the simulated MH climate. 145 The dynamical downscaling methodology applied here is a somewhat further developed version of the dynamical downscaling "pipeline" originally introduced in Gula and Peltier (2012) and then applied in some more recent studies (d'Orgeville et al., 2014;Erler and Peltier, 2016;Peltier et al., 2018;Huo and Peltier, 2019, where the expected climate change impacts over the Great Lakes Basin in North America, along the Rocky Mountains, in western Canada and over the South and Southeast Asia were investigated. In the present version of this pipeline, WRF-Hydro V5.1.1 (Gochis et 150 al., 2020) is coupled to WRF for the purpose of performing hydrometeorological simulations. The WRF-Hydro model is a multi-scale and multi-physics hydrologic community model widely used for flash flood prediction, seasonal simulation of water cycle components and regional hydroclimate impacts assessment both in the short term and long term. In this study, the Noah Land Surface Model (Tewari et al., 2004) is used in the coupled WRF-Hydro modelling system and high-resolution hydrological routing modules have been adopted to represent subsurface lateral flow (Gochis and Chen, 2003). Spin-up 155 simulations of 5 years have also been conducted to achieve numerical stability of the model outputs.

Observational datasets
To evaluate the historical simulations of temperature, the Climatic Research Unit (CRU) dataset V4.05 (Harris et al., 2020) at 0.5• resolution is employed. This widely used climate dataset was derived by the interpolation of monthly climate anomalies from several extensive networks of weather stations. To verify the precipitation results in the historical period, we 160 use both the CRU dataset and APHRODITE's (Asian Precipitation -Highly-Resolved Observational Data Integration Towards Evaluation) gridded precipitation dataset version 1101 (Yatagai et al, 2012), which covers Monsoon Asia (APHRO_MA_V1101) at 0.5•0.5• horizontal resolution. It is a set of long-term (1951 onward) continental-scale daily products that is based on a dense network of rain-gauge data for Asia including the Himalayas, South and Southeast Asia. To compute differences between model results and the observations, the latter have been reprojected and resampled onto the 165 native grid of each model. However, using these two different observational datasets didn't lead to significant differences to the major results of this study in terms of the validation of precipitation simulation ( Fig. 3 and Fig. A1). We have also compared the simulated temperature anomalies during the MH with paleoclimatic reconstructions based on lake sediment (Zhang et al., 2022;Kaufman et al., 2020), and precipitation anomalies with the pollen-based proxy reconstructions by Bartlein et al. (2011) andHerzschuh et al. (2019). Note here for the quantitative reconstructions of Holocene precipitation by 170 Herzschuh et al. (2019), we calculated the averages of all samples from each site within a rather broad time window (5.5-7.45 ka) to represent the MH, in order to be consistent with the method used for calculating the climate anomalies in Bartlein et al. (2011) and to include information from as many points as possible.
In addition, to evaluate the simulated hydrological metrics, mean monthly streamflow data from three stations have been obtained from the Global Runoff Data Centre (GRDC, 2015). These gauging stations (Tanglai Qu, Zhimenda and Yangcun) 175 are located at the basin outlets of the UYEB, UYAB and UBB ( Fig. 1) respectively and observations are available from 1978 to 1997, from 1978 to 1997 and from 1956 to 1982, respectively. There are some missing monthly data in these years of station records; however, all available streamflow data have been used to compute the average discharges discussed in Sect.

Historical and MH climates 180
In this section, we first validate our methodology against the observational dataset before presenting results from the MH simulations. Note that all temperature and precipitation biases and anomalies reported here are average values over the WRF inner domain, which covers the TP, as well as some surrounding regions. Figure 2 shows the June-July-August-September (JJAS) and December-January-February (DJF) temperature biases, 185 compared to the CRU dataset, in the of UofT-CCSM4 historical simulation and in the WRF inner-domain ensemble average for that same time period. The temperature biases are evidently dependent on the season. In JJAS the bias pattern is fairly consistent between the global and the regional models and includes a cold bias over the high elevations of the TP and a warm bias in the northeast of the TP and the lower lands around the TP, especially in the Tarim Basin. Most of the UYEB and UYAB are affected by this warm bias, which will likely reduce summer streamflow due to increased evapotranspiration.

7
Overall, summer temperatures are simulated fairly well by both models and integrated over the entire inner WRF domain, the average temperature bias in both models is less than 0.5•C. However, during winter, the GCM features a strong cold bias over the TP (> 3•C averaged over the inner domain), whereas the regional model has a better representation of DJF temperature with a significantly smaller average bias of approximately 1•C. Among the three river basins, the UBB is most strongly affected by the cold bias in winter, which may decrease streamflow in the cold season due to late snowmelt. Since 195 most of the precipitation over the TP occurs in summer, Fig. 3 shows the distribution of the JJAS average total precipitation in the observational dataset and the biases in the global model and the WRF ensemble. The spatial patterns in both models and observations are strongly dominated by orographic forcing with very high precipitation intensities on the south side of the Himalayas in contrast with the small rainfall intensity in the rain shadows behind the high mountains, especially over the northern TP which covers both the UYEB and UYAB. There is moderate rainfall over the Sichuan Basin east of the TP, and 200 the southern plateau, where the UBB is located. It is evident that the representation of orographic precipitation relies on model resolution: the coarse-resolution UofT-CCSM4 shows a severe underestimation of the peak intensities in front of the Himalayas and an overestimation of the precipitation rates in the northern TP. The representation of precipitation is greatly improved in WRF at 10-km resolution. The wet bias over the northern and eastern TP has also been eliminated in the RCM.
Averaged over the WRF inner domain, the JJAS precipitation bias is 0.9 mm d -1 in the WRF ensemble average and 1.7 mm 205 d -1 in the driving UofT-CCSM4 simulation. The overestimation in precipitation over the western and south-eastern TP in WRF is also accompanied by lower surface temperature, which is likely related to the greater cloud cover reflecting more shortwave flux at high levels. Both models show a similar wet bias in winter (Fig. 3d) and such excessive snow in winter possibly contributes to the lower DJF temperature (Fig. 2) through snow-albedo feedback.

Simulated MH climate 210
We next turn to an attempt to analyse the climate anomalies in MHREF and MHGS and to investigate the impact of the Saharan vegetation on the TP.
The simulated JJAS temperature over the TP in our MH ensembles is shown in Fig. 4. Owing to the variation in the Earth's orbital configuration, during the MH summer insolation in the NH was stronger than today, resulting in a significant warming over the TP in both the GCM and the WRF ensemble (Figs. 4a and 4b). The temperature response also shows 215 cooling south of the TP and over the ocean (Huo et al., 2021). Thus, the land-sea thermal contrast is increased in the MH favouring the enhancement of the monsoon circulations. The winter temperature over the TP shows a general cooling over the TP (Fig. 4e) that is associated with the decrease in winter solar radiation. Comparing the JJAS temperature changes produced by the GCM and RCM, the WRF ensemble shows an overall stronger warming over the TP, particularly in the northern part. Moreover, the regional model also produces a larger warm anomaly than the global model from mid-summer 220 to early-winter (Fig. 4e). Taking the influence of a vegetated Sahara into account results in higher temperature over the TP all year round, except in winter when the difference is negligible (Figs. 4d and 4e). The influence of MHGS in both models is characterized by an overall increase in the annual mean temperature over the TP, which agrees with estimates from pollen records (Shi et al., 1993), while an annual cooling is obtained with MHREF. Averaged over the WRF inner domain, the JJAS temperature anomalies of the WRF ensemble mean are 1.1 •C with MHREF and 1.8 •C with MHGS. Particularly, over the 225 norther plateau, the warming is significantly enhanced from approximately 2 •C to 4 •C, leading to a better agreement with proxy-based estimation of 4-5 •C warming (Shi et al., 1993) over the TP. Moreover, compared to paleoclimatic reconstructions based on lake sediment in Fig. 4 (Zhang et al., 2022;Kaufman et al., 2020), both MH experiments generally capture the summer warming trend over the TP. In south-eastern TP, there is a good agreement between the model simulations and the reconstructions. The simulated temperature anomalies in MHGS fit the reconstructed paleoclimatic 230 records in the north-eastern and south-central TP better but overestimate the warming signal over the central-eastern part of the TP (Fig. 4d). Meanwhile, the two proxy data points located in the central-western TP west of the 85° W indicate strong cooling (< -4 •C) during the MH, which disagrees with all model simulations. Note here these two records are from frozen lakes, where the reconstructed temperature records reveal air temperature changes during the ice-free season (May -September), not just JJAS. Also note here all proxy records are subject to uncertainties that arise from the dating processes 235 (Wang et al., 2021)  Under GS forcing (MHGS), precipitation is simulated to increase over almost the entire domain except a small patch in the southeast (Fig. 5d). The inclusion of GS surface boundary conditions greatly expands and intensifies the wet anomalies in WRF, further reducing the dry bias with respect to the pollen-based reconstructions (Bartlein et al., 2011) by half. The bias produced by the global model, however, only sees a small decrease from -0.73 mm d -1 to -0.71 mm d -1 , mainly due to the 250 dry bias over the eastern part of the TP (Fig. 5c). Compared to the reconstructed precipitation anomalies from Herzschuh et al. (2019), WRF also has a much smaller bias (0.02 mm d -1 ) than the global model (-0.12 mm d -1 ) in MHGS. Reconstructions based on lacustrine sediments (Shi and Song, 2003;Sun et al., 2006) implied an annual precipitation increase of 20-50% over Diaojiao Lake (112• E, 41• N) and Daihai Lake (113• E, 41• N) during the MH. WRF produced precipitation increases of 38% and 30% at these two locations when the influence of a vegetated Sahara is considered, which is consistent with the 255 reconstructed rainfall changes. However, the increase rates at these two locations are significantly underestimated by both UofT-CCSM4 in MHGS ( Fig. 5c; less than 20%) and the WRF simulation without GS forcing ( Fig. 5b; less than 5%).
Comparison of downscaled results including the GS with paleoclimate reconstructions shows significant improvements, especially over the north and west sides of TP. These results highlight the climate sensitivity to Saharan vegetation changes via ocean-atmosphere teleconnections and emphasized the importance of incorporating MH vegetation feedbacks. 260 The spatial patterns of JJAS precipitation anomalies resemble those of the annual mean but with a larger magnitude (Fig. 6), which is expected since annual precipitation largely results from the summer rainfall over the TP, a region heavily influenced by the Asian monsoon circulation. The seasonal cycle of precipitation obtained with MHREF changes in response to the insolation changes at 6 ka ( Fig. 6c): rainfall slightly decreases before May and then strengthens with the onset of the monsoon. Besides, the monsoon withdrawal is also delayed by one month as the increased rainfall persists till October. Such 265 changes are consistent with the results obtained in PMIP models (Zheng et al., 2013). Inclusion of the GS boundary conditions not only leads to a weaker decrease in rainfall than MHREF in spring, but also eliminates the negative rainfall The enhanced precipitation in MH experiment with Saharan vegetation is probably owing to ocean-atmosphere teleconnections as suggested in previous studies (Huo et al. 2021;Pausata et al., 2017). An albedo-induced warming 280 develops over the vegetated Sahara, leading to a strong intensification and northward expansion of the West African Monsoon and a significant tropical North Atlantic SST warming ( Fig. 7b; Pausata et al., 2016), which in turn changes atmospheric circulation and induces a notable intensification and westward extension of the Walker Circulation over the Pacific Ocean in the MHGS (Fig. 7d) compared to the MHREF experiment (Fig. 7c). The changes in the Walker Circulation weakens the low-level easterly winds over the eastern Pacific, but enhances easterly anomalies over the northern Indo-285 Pacific Ocean (Fig. 7b), which suppresses ENSO activity and enhances the Asian monsoon (Pausata et al., 2017).

Three river basins
In this section, the simulated climate for the present day as well as the MH over the three river basins on the TP are discussed. Figures 8-10 present the annual cycles of hydroclimatic variables over the UYEB, UYAB and UBB during the historical period and Figs. 11 and 12 show the anomalies during the MH. In Fig. 8-10, the WRF ensemble averages are 290 plotted in thick black lines while results from individual ensemble members are plotted in colors (see figure legends). The error bars on the observational data in Fig. 8-10 represent the standard error of the mean and the error bands in Figs. 11 and 12 shows the distribution range for all ensemble members. All hydroclimatic variables have been averaged over the corresponding river basin. Observational results in Fig. 8-10 are represented by open circles. Only temperature, total precipitation, and river discharge are available from observations. The observed temperature (precipitation) values are from 295 CRU (APHRODITE) dataset. The river discharges, recorded by gauging stations, have been normalized by the associated basin area. Areas of contribution for UYEB, UYAB and UBB are approximately 286000 km 2 , 138000 km 2 and 153000 km 2 , respectively. Solid precipitation in panel b of Figs. 8-11 is snow and graupel, and net precipitation (panel c in Figs. 8-10 and Fig. 12) is calculated by subtracting evapotranspiration from total precipitation. Note that evapotranspiration, snowfall, and snowmelt are calculated by the Noah land surface model, while streamflow is computed by the WRF-Hydro model. 300 For comparison between the preindustrial and MH simulations, we applied the linear scaling method to the original streamflow results. Due to the large (negative) JJAS runoff bias in WRF during the historical period, especially in the UYEB and UYAB (Figs. 8d and 9d), it was deemed necessary to apply a simple form of bias correction to rescale the monthly river discharge so that the monthly watershed average for the historical period matches the observations (using one factor per month). Linear scaling corrects the ensemble mean values of the historical simulations based on monthly errors (Trambauer 305 et al., 2015). The scaling factor was first obtained through calculating the ratio between the observed and modelled historical values. Then the monthly scaling factor was applied to each uncorrected daily runoff simulations of that month so that the monthly mean values of the historical simulations match those of the observation. The streamflow values of the MH simulations have been bias corrected using the same factors as for the historical results so that they can be usefully compared.

Validation of the average season cycle of hydroclimatic variables 310
The seasonal cycles of hydroclimatic variables over the three river basins are provided in Figs. 8-10. The Yellow River originates in the northern TP and flows eastward through the Loess Plateau and the North China Plain, eventually draining into the Pacific Ocean at the Bohai Gulf on the east coast of China. The UYEB above the Tanglai Qu hydrometric station, with an altitude varying from 4.9 km in the southwest to 1.1 km in the northeast, is situated on the northeast side of the TP but also covers a small section of the western Loess Plateau (Fig. 1). It has over 400 mm annual rainfall according to the 315 APHRODITE dataset and an annual mean temperature of approximately 2 •C based on the CRU dataset. The UYAB lies in the eastern part of the TP and is influenced by the Asian monsoon. The basin, whose altitude varies from 5.7 km in the west to 4.2 km in the east, drains past the Zhimenda hydrological station (Fig. 1). The integrated precipitation in JJAS accounts for more than 85% of the total annual precipitation (Fig. 9b) and the JJAS discharge accounts for more than 70% of the annual total at Zhimenda (Fig. 9d). The Brahmaputra River originates from the Himalayas and passes through the southern hydrometric station ranges between 5.8 km and 3.8 km and drops off from the west to the east. Based on the observation, the basin-averaged annual precipitation is 314 mm and the annual mean temperature is −1 •C.
The seasonal cycle of air temperature is generally accurately reproduced over the three river basins, except for a cold bias during winter, even though the downscaled simulations have already halved the DJF cold bias in the global simulation. The 325 large cold bias in UofT-CCSM4 is likely due to excess winter precipitation , which results in a snowalbedo feedback on temperature. It is interesting to note that the WRF ensemble actually suffers from a stronger cold bias than the UofT-CCSM4 in spring over the two cooler basins (UYAB and UBB), which is likely associated with the poor representation of snowmelt processes by the much simpler snow model employed in WRF. The JJAS temperatures simulated by both the global and the regional models fit the observation well. 330 The Asian monsoon exerts strong influence on the three river basins over the TP: over 70% of annual precipitation occurs in JJAS and only about 10% falls in the winter (DJF) according to the observation. The simulated seasonal cycle of precipitation in WRF is too strong compared to the observation over the UYEB and UBB, largely due to the rainfall overestimation in summer (wet season). This bias is inherited from the GCM and although downscaling greatly reduces the bias, precipitation is still overestimated for JJAS. Over the UYAB, qualitatively the annual precipitation pattern is well 335 reproduced by the WRF ensemble, while UofT-CCSM4 is characterized by a much greater wet bias, with significantly more precipitation in spring and summer (Fig. 9b). The improvement in simulation of both regional temperature and precipitation by WRF indicates the fidelity brought by the use of higher resolution. There is, as expected, little solid precipitation in summer and all precipitation falls in solid form during winter in UYEB, resulting in a streamflow peak driven by snowmelt in late spring and early summer (Fig. 8c). Summer precipitation also primarily falls as liquid rain in UYAB and UBB, but 340 there is still snow even in July and August because of the colder climate over these two basins compared to the UYEB and snowfall in UBB is in fact relatively evenly distributed throughout the year. Note here the spread among the ensemble members is larger for liquid precipitation than for solid precipitation or temperature; different convection schemes greatly affect the JJAS precipitation and the second ensemble member using the GF cumulus scheme exhibits the largest wet bias over the two basins in the northern TP whereas the BMJ scheme produces the largest excess in precipitation over the UBB. 345 The seasonal cycle of net precipitation pattern in UYEB and UYAB has two peaks (June and September), due to the large evapotranspiration in July and August (Figs. 8c and 9c). Snowmelt over these tow basins is also characterized by double peaks, which correspond to the two peaks in snowfall (Figs. 8b and 9b). Therefore, the melting of fresh snowfall is likely the major contributor to the snowmelt peaks simulated by the model. Over the UBB, although both observed and simulated total precipitation starts to increase in June, net precipitation is small in early summer due to evapotranspiration (Fig. 11c). Net 350 precipitation from the first WRF ensemble member with the KF cumulus scheme is in fact negative in May and June. Runoff over the UBB starts to increase in May (Fig. 10d), which possibly comes from the snowmelt contribution that is the largest in May and June (Fig. 10c). In UBB, the runoff maximum over the UBB in August coincides with the largest net precipitation rates in July and August, which indicates that UBB rainfall plays a dominant role in peak runoff during the after-snowmelt seasons. The seasonal streamflow during the historical period for the two basins in the north exhibits a broad 12 shape with flow increasing slowly in April and slowly declining after October. Runoff in spring and early summer clearly follows the timing of snowmelt. Streamflow regime is also characterized by double peaks, in July and September, respectively (Figs. 8d and 9d). Over the UYEB, the good correspondence between the timing of the peaks of runoff and net precipitation in September indicates that the monsoon rainfall is the dominant contributor to streamflow for UYEB. However, over the UYAB, both snowmelting and precipitation contribute to the second runoff peak in September. WRF-Hydro is able 360 to reproduce runoff that broadly agrees with observed river discharge. However, the model tends to underestimate the magnitude of the runoff peaks. This is likely a result of the stronger evapotranspiration, leading to a lower runoff. The mismatch in streamflow peak magnitude against the observation is similar in all WRF simulations based on different physics configurations except the second ensemble member using the GF cumulus scheme; the GF cumulus scheme overestimates the runoff all year round and also produces the largest net precipitation and snowmelt among all four WRF ensemble 365 members.

Simulated MH changes
During the MH, the overall changes in mean temperature (Fig. 11a-c) follow changes of insolation the magnitude of which increases during boreal summer and decreases in winter. The WRF model produces significantly larger positive temperature anomalies in summer compared to UofT-CCSM4. Compared with MHREF, temperature changes in experiments with the GS 370 forcing have very similar characteristics, but the simulated temperature is generally warmer than in MHREF over all river basins in both global and regional models. Warming also starts earlier in May and ends one month later in November and the peak warming also occurs one month later in early autumn in our WRF MHGS experiments relative to the MHREF. In both MH ensembles, different WRF configurations show very similar temperature changes.
The changes in monthly precipitation are shown in Fig. 11d-f. The MHREF simulations for both UYEB and UYAB show a 375 significant decrease in precipitation except in autumn. The absolute regional dry anomaly is the largest in August, while September (October) experiences the largest increase in rainfall in the UYEB (UYAB). When vegetation is imposed over the Sahara, the minimum in August and the maximum in September (October) still exist, but there is another smaller precipitation increase peak in June (July), showing a strong rainfall enhancement during the early monsoon season in response to the GS forcing. However, averaged over the whole monsoon season, the WRF ensemble mean in MHGS only 380 shows a small increase over these two basins (~ 5%) due to the offset between these positive and negative anomalies, while the WRF simulations that only account for the insolation forcing (MHREF) produce a decrease in monsoonal precipitation (9%). Another major change to the hydroclimate over the UYEB and UYAB is the shift of part of the solid precipitation to liquid form in summer in MHGS when compared to MHREF (Fig. 11d and e), which will in turn result in a decreased snowmelt ( Fig. 12a and b). Over the UBB, annual precipitation in the both MH experiments increases (~0.15 mm d -1 or 11% 385 increase for MHREF and 0.35 mm d -1 or 26% increase for MHGS), which is larger than those over the other two basins.
However, although the increasing trend is robust among all WRF ensemble members, the magnitude of the increase is subject to higher model uncertainties than the other two basins. Compared to the increase in total precipitation, there is little change in snow amount for both MH experiments (< 0.04 mm d -1 ) over the UBB, which also leads to very small changes in the annual snowmelt (Fig. 12c). Over the other two basins, however, the reduction in peak snowmelt due to the combined 390 effects of MH insolation forcing and GS forcing is over 30% (Fig. 12a and b). At the same time higher temperatures during the warm season resulting from the MH insolation forcing also lead to larger evaporation, which accounts for the drop in net precipitation and also the substantial decrease in discharge especially in summer ( Fig. 12d and e). The magnitude of streamflow reductions in slightly larger in MHGS than MHREF, a consequence of the higher temperature and less snowfall in MHGS, which leads to greater evaporation and less snowmelt. Compared with UYEB, the reduction in snowmelt plays a 395 more important role in the runoff decrease over the UYAB. Such reduction in the upper Yellow River runoff has also been found over the past decades, which is related to rising temperatures and decreasing wet season precipitation (Cuo et al., 2013). The simulated negative anomaly in JJAS discharge of the upper Yangtze River also agrees with CMIP5-based projections (Krysanova et al., 2017) and hydrological model results forced by PMIP4 data during the last interglacial (LIG; Scussolini et al., 2020). Among the three basins, one major difference in the hydroclimate over the UBB is the increase of 400 net precipitation in JJAS, which also results in an overall increase in runoff. In MHREF, streamflow first decreases very slightly in April and May and then starts to increase in July. In WRF simulations that consider the influence of the Sahara greening, runoff exhibits similar increasing trends from July to September, which is consistent with changes in precipitation, but there is no decrease in May, probably due to the enhanced snowmelt compared to MHREF. Similar river runoff increase has also been found in future projections in a warmer climate using a hydrological model and multi-source spatial data (Cai 405 et al., 2017) and LIG hydrological simulations based on PMIP4 data (Scussolini et al., 2020).

Summary
Using dynamically downscaled high-resolution climate simulations over the TP, this study analyses the hydroclimate in three major river basins, the UYEB, UYAB and UBB. Climate simulations constructed using the UofT-CCSM4 global climate model for the modern and MH time periods were dynamically downscaled to 10-km resolution, using coupled WRF-WRF-410 Hydro with four different cumulus parameterization schemes. Compared to the GCM, WRF is far superior in representing orographic precipitation and the annual cycles of hydroclimatic variables. It captures the sharp contrast between the very high precipitation rates south of the Himalayan Mountain ranges and the dry climate over the interior plateau. Both global and regional models suffer from a winter cold bias owing to excessive snowfall, but the regional model has again significantly reduced this cold bias. 415 The coupled UofT-CCSM4-WRF-WRF-Hydro model has also been employed to investigate the temperature and precipitation changes over the TP during the MH. Both the GCM and the WRF ensemble average show a warming in summer and a cooling in winter over the TP, which is probably induced by insolation change. As far as the JJAS temperature is concerned, the WRF simulations show a larger warming than that in the driving global model, especially in the northern part of TP, and is thereby closer to estimates from proxy records. The strengthened thermal contrast between continent and ocean during summer enhances the monsoon circulation, leading to stronger northward winds, as well as stronger moisture transport from the sea. The spatial patterns for the TP precipitation change are generally similar in sign between the GCM and the RCM, with significant increase in the southern parts of the TP. The annual precipitation cycle over the TP is changed by the insolation anomalies during the MH, such that precipitation weakens slightly before the monsoon season and strengthens in summer. The annual wet anomalies are stronger and more extensive in the downscaled simulations and thus fit 425 better with the proxy data. Since the TP plays a profound role in the development of the Asian monsoon and heavy precipitation forms when the wind from the sea rises along the south and east edges of the plateau, the model's ability to realistically represent the complex steep topography may have a great impact on the coupling between the general monsoon circulations and the regional changes over the TP. Including a GS further strengthens the TP precipitation and significantly reduces the bias against reconstructions over the TP, which highlights the significance of the teleconnections between the 430 Saharan vegetation and the TP in paleoclimate simulations. Saharan vegetation plays a crucial role in intensifying the West African Monsoon and modulating the atmospheric circulation, which alters the Walker circulation and increases Asian monsoon precipitation, through changes in equatorial Atlantic SSTs (Fig. 7). Precipitation changes in MH experiments which assume preindustrial vegetation cover are approximately half the magnitude of those in MHGS in all WRF ensemble members, although there are notable differences in the simulated MH rainfall anomalies using different physics schemes. 435 Both UYEB and UYAB hydrological regimes exhibit changes in the MH as manifested by a decline in streamflow, especially in summer. Such flow decreases are a combined effect of changes in snowmelt and evapotranspiration. A significant amount of solid precipitation shifts to liquid precipitation and the JJAS net precipitation is simulated to shrink in both MH experiments. The GS forcing caused a rise in temperature during the MH, as well as larger rainfall but smaller snowmelt and lager evaporative water losses compared to the MHREF. In the UBB, the simulated annual total precipitation 440 increase in the MH is the largest among three basins, and, unlike the other two basins, the most significant MH hydroclimatic anomaly over the UBB may be an increase in runoff in both MH experiments, particularly in mid-summer, due to greater net precipitation. The greening of the Sahara led to higher temperature and enhanced snowmelt in spring and eliminated the drop in runoff in April in MHREF.
One of the most fundamental uncertainties that limit the validity of the MH simulations in this study is natural variability, 445 including the El Niño-Southern Oscillation. This calls for an initial condition ensemble where the RCM is forced with different periods of the global simulations. Besides, this study has used only one coupled regional modeling system to investigate the influence on the TP climate only due to the MH land cover changes over northern Africa. However, reduced dust emissions from the Sahara during the MH have also been suggested to play a role in modulating the Asian monsoon Pausata et al., 2020). Moreover, reconstructions suggest extensive land cover changes during the MH, in 450 particular higher vegetation coverage on the Eurasian continent (Tarasov et al., 1998;Zhang et al., 2014;Li et al., 2019;Chen et al., 2020b), which has been shown to have global influence by shifting the latitudinal position of the of the intertropical convergence zone (Swann et al., 2014). These all underscore the necessity of further studies over the TP based on different climate models, as well as more realistic dust concentrations and land cover in the paleoclimate simulations.

Code/Data availability
Code and data for reproducing the figures in the article can be obtained from Yiling Huo.