CMIP6 Simulation-Based Daily Surface Air Temperature and Precipitation Projections over the Qinghai-Tibetan Plateau in the 21st Century

: The Qinghai-Tibetan Plateau (QTP), the source of many major Asian rivers, is sensitive to climate change, affecting billions of people’s livelihoods across Asia. Here, we developed high-resolution projections of precipitation and daily maximum/minimum temperatures at 0.1 ◦ spatial resolution over the QTP. The projections are based on the output from seven global climate models (GCMs) from the Coupled Model Intercomparison Project Phase 6 (CMIP6) for historical (1979–2013) and projected (2015–2100) climates across four scenarios (SSP1-2.6, SSP2-4.5, SSP3-7.0, and SSP5-8.5). An updated nonstationary cumulative distribution function matching method (called CNCDFm) was used to remove model systemic bias. We verify the necessity of taking into account altitude in downscaling processes and the validity of nonstationary bias correction. Compared to the historical period, the climate in the QTP in the 21st century is warmer (1.2–5.1 ◦ C, for maximum surface temperature) and wetter (3.9–26.8%) according to the corrected GCM projection. For precipitation, the Indus River (IDR), Tarim River (TMR), Inner of Qiangtang Basin (IQTB), Yarlung Zangbo (YLZBR), and Qaidam Basin (QDB) showed growth well above the global average across high radiative forcing scenarios, which could have a profound impact on the regional hydrological cycle. However, there is great uncertainty in precipitation prediction, which is demonstrated by a very low signal-to-noise ratio (SNR) and a large difference between Bayesian model averaging (BMA) and multi-model averages (MMAs). This bias-corrected dataset is available for climate change impact research in the QTP at the subregion scale.


Introduction
The Qinghai-Tibetan Plateau (QTP), called the Asian Water Tower, is the source of many major Asian rivers, such as the Yellow, Yangtze, Mekong, Indus, Salween, and Brahmaputra [1].Climate change in the QTP affects billions of people in Asia by changing regional water resources, agricultural production, energy security, and disaster risk [2][3][4][5].That is why the QTP has attracted the attention of academics and policymakers.In previous studies, the QTP was found sensitive to climate change.Based on historical data and reanalysis data, previous research showed that the QTP area had a more significant increase in temperature than the surrounding areas [6].However, it is difficult to estimate the spatial distribution of changes in meteorological elements.Because of complex topographic conditions, meteorological factors of the QTP have high spatial heterogeneity, especially in near-surface areas [7,8].These distribution characteristics may exhibit different evolutionary trends under different climatic backgrounds [9].Furthermore, low-resolution meteorological data lack representativeness and practical value on the plateau [10].
Numerous downscaling and multi-source data methods have been applied in largescale regions to provide high-resolution data.In recent years, several high-resolution reanalysis datasets have been evaluated and applied in the Qinghai-Tibetan Plateau (QTP), with precipitation and temperature being the most crucial variables [11][12][13][14][15].These studies not only demonstrate the possibility of generating high-resolution meteorological data in ungauged regions but also provide forcing data to drive a variety of terrestrial models, which is essential for climatic impact studies [16][17][18].Nevertheless, due to differences in input data, numerical assimilation models, parameterization schemes, and the spatiotemporal resolutions of final products, these reanalysis datasets exhibit varying performances in different regions [19,20].Therefore, the question of how to apply the reanalysis data remains worthy of discussion.
Global climate models (GCMs) are effective tools for simulating Earth's climate on continental scales and larger [21].They are the most important tools and infrastructure for projecting future climatic conditions and assessing the climate change impact [22].GCMs have been widely used in climate change research [23][24][25], showing growing potential when used as forcing data for land surface models and hydrological models.The Coupled Model Intercomparison Project Phase 6 (CMIP6) indicates an overall improvement in GCMs over the last decade [22,26].Typically, the cold bias of CMIP5 on air temperature is significantly reduced in CMIP6 for the QTP [27].
However, major problems still exist and prevent researchers from using GCM results directly.At first, the models in CMIP6 still have considerable systematic biases in meteorological simulation of regional scale.Lun et al. [28] found that the multi-model ensembles (MME) of 23 GCMs in CMIP6 overestimated the QTP precipitation by about 84.1% in historical simulation, and the wet biases of the same models were even higher than historical observation.Secondly, the outputs of GCMs still feature a low resolution that is difficult to align with the input requirements of land and hydrological models, especially when compared to the more detailed scales needed for accurately simulating land and hydrological processes [29].In general, the resolution of GCM outputs is between 0.5 • and 2.5 • , which is much larger than most land and hydrological processes.For instance, Su et al. [2] developed a hydrological model to simulate the impact of glacier melt in the QTP, and the model was set up on 0.083 • × 0.083 • grids because glacier areas are from 100 km 2 to 10,000 km 2 .Hence, downscaling GCMs has become a pre-processing step in many studies, but existing studies tend to use the linear interpolation method in downscaling [2,5,[30][31][32], ignoring the influence of the underlying surface on the spatial distribution of climate variables, which may cause further systematic errors.Moreover, GCMs still exhibit large uncertainty that shows up in different models and produces widely different results when they simulate the same variable in the same forcing scenario.In the QTP, Zhang and Yang [27] found that GCM-simulated daily precipitation (including both liquid and solid phases, units: mm d −1 ) values have large variation, with most values ranging from 1.4 mm d −1 to 3.5 mm d −1 .Such large variations make it difficult for users to decide which model is the best model or how to ensemble prediction based on several models.These problems become worse when looking at particular uses of GCM simulation output [22].But on the other hand, although it is still challenging to obtain optimal projection over the QTP, smaller deviation forecasts (compared to the raw output from CMIP6) based on effective statistical methods are essential for regional-scale studies.
Therefore, this study attempts to produce a 0.1 • × 0.1 • high-resolution projection over the entire QTP (HRP-QTP) based on reanalysis data and CMIP6 simulations, which includes the historical period (from 1979 to 2013) and future period (from 2015 to 2100), and can be used in watershed scale studies (variables information presented in Table 1).In the present study, combined reanalysis data were applied as historical data for ungauged regions.In order to deal with the complex terrain conditions of the QTP, we adopted a downscaling method different from the linear interpolation commonly used in previous studies to fit the effects of altitude on climate change.Furthermore, an updated bias correction method was used to remove systematic bias from GCM projections.The present study also produced ensemble data from several models.The projection is available for researchers who want to learn about climate change and its impact on the QTP.It also provides a robust and universal method to deal with CMIP6 data, which can be used in other regions confronting climate change risk.

Definition and Division of Study Area
As for the scope of the Qinghai-Tibet Plateau, there are different definitions in different studies.We adopted the high-frequency use scope of the QTP (HF-QTP) as defined by Zhang et al. [33], which included several major river headwaters and can match most research requirements.The QTP ranges from 25 • 59 ′ 37 ′′ N to 39 • 49 ′ 33 ′′ N and from 73 • 29 ′ 56 ′′ E to 104 • 40 ′ 20 ′′ E, covering an area of 2542.30× 10 3 km 2 [34].Its altitudes range from 85 m to 8233 m, and approximately 90.9% belong to high altitude (>3000 m) areas.This range includes the world's highest mountains and massive deep valleys in terms of topography.The complex terrain of the QTP results in a complicated distribution of precipitation and surface temperature.A brief overview can be seen in Figure 1.These regions can be divided into ten subregions based on watershed [35,36], as follows: I. Tarim River (TMR); II.Qaidam Basin (QDB); III.Hexi Basin (HXB); IV.Yellow River (YR); V. Inner of Qiangtang Basin (IQTB); VI.Yangtze River (YZR); VII.Indus River (IDR); VIII.Salween River (SWR); IX.Mekong River (MKR); and X. Yarlung Zangbo (YLZBR).Among these, IQTB and QDB belong to the inland river basin, and some of the river channels are seasonal.We obtained the river distribution of QDB from National Basic Geographic Information [37].Currently, there are only 131 long-observed meteorological stations in the QTP (Figure 1) maintained by the China Meteorological Administration (CMA).They are primarily situated in the eastern part of the QTP, with observations lacking in the western part of the plateau.The digital elevation model (measured from sea level, DEM, units: m) used in this section and Section 2.5 is provided by the National Tibetan Plateau Data Center "https://data.tpdc.ac.cn/ (24 March 2021)" [38].
In Figure 1, the lower Yarlung Zangbo valley (LYZV) is highlighted with red dots.Its boundary is defined by the watershed of lower Yarlung Zangbo and is the joint part of the Yarlung Zangbo River and Brahmaputra River.Unlike other regions in the QTP, LYZV belongs to the south of the Himalayas, is controlled by the Indian monsoon [39][40][41][42][43], and has greatly varied elevations (from 85 m to 6941 m).Because of water vapor transport by monsoons and the uplift effect of the valley, precipitation in LYZV is greater than in other regions [44].Although LYZV is an ungauged region, we can estimate its annual precipitation ranging from 2500 mm to 3500 mm by using observations from similar neighborhoods in the south of the Himalayas and early investigation findings [43,[45][46][47].Thus, it is essential to provide reasonable estimations in LYZV.Special methods are required to deal with meteorological data in LYZV.We designed a method that combines two highresolution grid datasets to reduce the underestimation of precipitation, which can be seen in Section 2.3.

Historical Meteorological Data
(1) Observation data from CMA.
As Figure 1 demonstrates, the Qinghai-Tibet Plateau is indeed a data-sparse region, with limited and unevenly distributed meteorological stations across the area.In the CMA published ground observation data (from "http://data.cma.cn(24 March 2021)"), there are about 147 observation meteorological stations with continuous long series on the Qinghai-Tibet Plateau, and the density of stations is 0.058 per 1000 km 2 .Ma et al. [12] have demonstrated the impact of overly sparse sites on downscaling representativeness, so in previous studies, CMA site observations were generally used to validate reanalysis methods and results [13][14][15]19,20].In the present study, CMA data were used to verify the availability of the downscaling method (Section 2.6).Among them, 8 stations are used as verification stations to represent climate characteristics at different altitudes (white dots in Figure 1) and the remaining 139 stations are interpolation stations (black dots in Figure 1).The precipitation and temperature data of interpolation stations from 1979 to 2013 were interpolated into the verification station grids (the grids are consistent with HRP-QTP) by different methods to compare the fitting effects of different interpolation methods on the verification stations.Further, a suitable interpolation method is selected for downscaling.Basic information on the verification stations is shown in Table 2.

Historical Meteorological Data
(1) Observation data from CMA.
As Figure 1 demonstrates, the Qinghai-Tibet Plateau is indeed a data-sparse region, with limited and unevenly distributed meteorological stations across the area.In the CMA published ground observation data (from "http://data.cma.cn(24 March 2021)"), there are about 147 observation meteorological stations with continuous long series on the Qinghai-Tibet Plateau, and the density of stations is 0.058 per 1000 km 2 .Ma et al. [12] have demonstrated the impact of overly sparse sites on downscaling representativeness, so in previous studies, CMA site observations were generally used to validate reanalysis methods and results [13][14][15]19,20].In the present study, CMA data were used to verify the availability of the downscaling method (Section 2.6).Among them, 8 stations are used as verification stations to represent climate characteristics at different altitudes (white dots in Figure 1) and the remaining 139 stations are interpolation stations (black dots in Figure 1).The precipitation and temperature data of interpolation stations from 1979 to 2013 were interpolated into the verification station grids (the grids are consistent with HRP-QTP) by different methods to compare the fitting effects of different interpolation methods on the verification stations.Further, a suitable interpolation method is selected for downscaling.Basic information on the verification stations is shown in Table 2. (2) Reanalysis data.
It is crucial to select appropriate meteorological data for the historical period because multi-model ensemble prediction and bias correction are highly dependent on historical data.In other words, historical data determine whether the output reanalysis data are reasonable and available.We used two high-resolution reanalysis (0.1 × 0.1 • ) datasets from He et al. [11] (the China Meteorological Forcing Dataset, CMFD, from 1979 to 2013) and Ma et al. [12] (Gridded Precipitation for Quantile Mapping, GPQM, from 1980 to 2009) as historical data in our study.Data were downloaded from the National Tibetan Plateau Data Center "http://data.tpdc.ac.cn (24 March 2021)"; basic information about CMFD and GPQM is presented in Table 3.

Combination of Reanalysis Data
CMFD provides a high-resolution meteorological forcing dataset across China and has been evaluated for its temperature performance in western China.It shows better performance in surface temperature results in the QTP's ungauged regions compared to other reanalysis data [11]; hence, the present study uses CMFD as the historical temperature source directly.However, CMFD precipitation is unverified in the western part of the QTP because of the lack of surface observational data, and it shows obvious bias in specified areas.What stands out is that it underestimates the precipitation in the LYZV compared with observations from ground stations in the Motuo and its immediate surrounding areas [48].Hence, we adopt GPQM as a historical precipitation source in LYZV, which is validated with no obvious bias.GPQM considers observation stations both inside and surrounding the QTP [12].It includes decades of gauge data from the south of the Himalayas, which can comprehensively reflect the precipitation characteristics of the QTP.To align the GPQM with the historical period of the HRP-QTP, we extend it using CMFD, and the extended method is shown in Figure 2.
Figure 2 employs three background colors to denote grid data from distinct sources: orange for GPQM data, purple for CMFD data, and green for the final 0.1 • × 0.1 • grid output of QTP's precipitation data.Initially, we extracted QTP's precipitation data, excluding the LYZV areas, from CMFD to serve as the historical precipitation source for other regions.Then, precipitation data for the LYZV were extracted from GPQM to the same grids as CMFD.Given that GPQM and CMFD share the same resolution, we applied GPQM data grid by grid to correct CMFD within the LYZV region.This correction allowed us to extend the GPQM data, ensuring that precipitation data from both sources matched the time dimension.We utilized an updated nonstationary CDF-matching method (CNCDFm), detailed in Section 2.5, for this correction.Finally, the extended GPQM precipitation data for LYZV were combined with the CMFD-extracted precipitation data for other regions, creating a comprehensive historical precipitation dataset for the entire QTP.This combined dataset will serve as the observed precipitation in subsequent steps.
Atmosphere 2024, 15, 434 6 of 30 the time dimension.We utilized an updated nonstationary CDF-matching method (CNCDFm), detailed in Section 2.5, for this correction.Finally, the extended GPQM precipitation data for LYZV were combined with the CMFD-extracted precipitation data for other regions, creating a comprehensive historical precipitation dataset for the entire QTP.This combined dataset will serve as the observed precipitation in subsequent steps.

Simulation and Prediction Data from GCMs
Seven CMIP6 GCMs were used to produce HRP-QTP (See Table 4).The GCM outputs were obtained from "https://esgf-node.llnl.gov/search/cmip6/(24 March 2021)".These GCMs were selected because precipitation, maximum temperature, and minimum temperature data were available for the 7 GCMs across most future scenarios.Moreover, they represent different resolutions, from 1.12° to 2.5°.In an assessment by Lun et al. [28], the 7 GCMs overestimate the precipitation of the QTP by 0.8-1.2mm d −1 in the historical period, compared with CMFD.We interpolated them to a resolution of 0.1° and used historical data to correct their systemic bias.The CMIP6's global simulation results of these climate models come from institutions worldwide, using different modeling methods.Compared with the representative concentration pathways (RCPs) of CMIP5, CMIP6 developed a new framework called shared socioeconomic pathways (SSPs) to quantitatively describe the relationship between climate change and socioeconomic development [26].In this study, we selected SSP1-2.6,SSP2-4.5, SSP3-7.0, and SSP5-8.5 to represent different development and concentration pathways.These scenarios correspond to a sustainable development scenario, a medium challenge scenario for mitigation and adaptation, a high challenge scenario for mitigation and adaptation, and a high radiative forcing scenario with fossil fuels in the future period [49].The four SSPs present varied global scenarios of radiative forcing, ranging from low to high.SSP1-2.6,SSP2-4.5, and SSP5-8.5 each align with the radiative forcing assumptions of RCP2.6, RCP4.5, and RCP8.5, respectively.In contrast, SSP3-7.0 serves as an alternative to the high radiative forcing scenario.Although it projects lower greenhouse gas concentrations than SSP5-8.5, it suggests that human society will face a greater risk of climate hazards [50].We chose the first 3 SSPs to align our results with those of CMIP5 and SSP3-7.0 to represent a more likely high radiative forcing scenario.

Simulation and Prediction Data from GCMs
Seven CMIP6 GCMs were used to produce HRP-QTP (See Table 4).The GCM outputs were obtained from "https://esgf-node.llnl.gov/search/cmip6/(24 March 2021)".These GCMs were selected because precipitation, maximum temperature, and minimum temperature data were available for the 7 GCMs across most future scenarios.Moreover, they represent different resolutions, from 1.12 • to 2.5 • .In an assessment by Lun et al. [28], the 7 GCMs overestimate the precipitation of the QTP by 0.8-1.2mm d −1 in the historical period, compared with CMFD.We interpolated them to a resolution of 0.1 • and used historical data to correct their systemic bias.The CMIP6's global simulation results of these climate models come from institutions worldwide, using different modeling methods.Compared with the representative concentration pathways (RCPs) of CMIP5, CMIP6 developed a new framework called shared socioeconomic pathways (SSPs) to quantitatively describe the relationship between climate change and socioeconomic development [26].In this study, we selected SSP1-2.6,SSP2-4.5, SSP3-7.0, and SSP5-8.5 to represent different development and concentration pathways.These scenarios correspond to a sustainable development scenario, a medium challenge scenario for mitigation and adaptation, a high challenge scenario for mitigation and adaptation, and a high radiative forcing scenario with fossil fuels in the future period [49].The four SSPs present varied global scenarios of radiative forcing, ranging from low to high.SSP1-2.6,SSP2-4.5, and SSP5-8.5 each align with the radiative forcing assumptions of RCP2.6, RCP4.5, and RCP8.5, respectively.In contrast, SSP3-7.0 serves as an alternative to the high radiative forcing scenario.Although it projects lower greenhouse gas concentrations than SSP5-8.5, it suggests that human society will face a greater risk of climate hazards [50].We chose the first 3 SSPs to align our results with those of CMIP5 and SSP3-7.0 to represent a more likely high radiative forcing scenario.

Methods of Processing CMIP6 Data
Methods used to process CMIP6 data are divided into three steps: downscaling, multimodel ensemble forecasting, and bias correction.At first, we used a statistical method to downscale 7 GCM outputs from different resolutions.Secondly, a prior-based method was used to generate an ensemble forecast from GCMs.Finally, all results from the downscaling and multi-model ensemble were corrected by a nonstationary method.Details of the three steps are as follows: (1) Downscaling.
Most previous studies used linear interpolation methods to downscale GCM data and ignored the relationship between altitude and meteorological elements [2,4,25].However, orographic studies demonstrated an approximately linear correlation between precipitation and temperature with elevation in high-altitude regions [51].Several studies showed that it is necessary to consider altitude when generating high-resolution meteorological forcing [52,53].Hence, we utilized ANUSPLIN (V4.36) to downscale the GCM data.ANUS-PLIN is a software package that is based on the theory of thin-plate smoothing splines, and is an extension of linear regression, which considers distance, altitude, and other variables.Ma et al. [12] confirmed that ANUSPLIN is more accurate than inverse distance weighted and ordinary kriging methods in monthly precipitation interpolation over the QTP.In this study, ANUSPLIN was used to downscale daily precipitation and surface temperature, and altitude was considered a covariate in ANUSPLIN.More details about ANUSPLIN can be seen in the paper from M.F.Hutchinson and P.E.Gessler [54].
Due to the uncertainty in GCMs, the results of different climate models are widely distributed, and it is necessary to provide certain projection results for specific studies [55].In this research, Bayesian model averaging (BMA) was used to produce a multi-model ensemble forecast result based on 7 GCM outputs in the historical period and 4 SSPs.BMA is a statistical method that is based on the posterior probability density functions (pdfs) of model predictions [56].BMA, which is different from the multi-model average (MMA), could be described as a weighted average, as follows: where M BMA is the output result of BMA, M k is the output of model k, and w k is the weight of model k, denoting the posterior probabilities of model k conditioned.If we assume that the actual climate change process, y, is a combination of k model outputs, w k could be estimated based on fitting the data to the following relationship: where the left side of Equation ( 2) denotes the pdf of the truth value in the historical period, and p(y|M k ) denotes the conditional probability distribution of M k ; it represents the probability that y follows the pdf of the output from model k.When p(y|M k ) is assumed to follow a normal distribution, it could be estimated as follows: with y(t) and M k (t) representing the y and M k values at time, t, respectively, and σ 2 is the variance of y − M k .In this study, a Markov Chain Monte Carlo (MCMC) method proposed by Jasper et al. [57] was used to sample the weights of 7 GCM outputs.The weights were calculated based on historical GCM outputs in every 0.1 • grid to maintain the spatial distribution characteristics of meteorological variables.The grid's daily temperature and precipitation data obtained in Section 2.3 were used as truth values.BMA, which differs from taking the mean value of several models, considers the model's performance during the training period, which could reduce the uncertainty of the ensemble result effectively [56][57][58].MMA was also calculated in this study to compare with BMA.
Because climate processes are complex and cannot be described accurately by existing methods, all GCMs still have considerable biases in their output results.Previous studies have adopted many statistical methods to eliminate the systematic errors of the model, such as linear and nonlinear corrections or cumulative distribution functions (CDFs) [31,59,60].Moreover, it is important to retain nonstationary information in future scenarios and avoid the occurrence of abnormality during bias correction.In this research, we adopted a nonstationary bias correction technique proposed by Miao et al. [61], called updated nonstationary CDF-matching (CNCDFm), to correct downscaled GCM data and BMA results.CNCDFm model methods can be presented as follows: where x denotes the meteorological variable from historical observations (o, denoted here as the reanalysis data obtained from Section 2.3) or models (m) during the historical period (h), current climate (c), or future projection (p) periods, F(•) is the CDF for a variable, and F −1 (•) is its inverse.CNCDFm is an updated version of the CDF technique.By combining the method by Li et al. [62] with the method from Wang and Chen [63], it effectively retained nonstationary information in future scenarios and avoided the occurrence of abnormality.For instance, this approach could avoid negative values for precipitation or unreasonably high values (when the denominator on the right-hand side of Equation ( 6) is very small).Miao et al. [61] provided validation within the CMIP5 to approve the method's efficiency.
Estimating the probability distribution of CDF is the key point of CNCDFm.We used the empirical cumulative distribution function (ECDF) to estimate the CDF in all periods for robust application.We used the 1979-2013 period as a training period to estimate F o (•) and F m-p (•) by the ECDFs of observed data and downscaled model output (obtained from Section 2.5 "Downscaling").The wet day threshold was set at 0.1 mm/day to correct for the drizzle effect in precipitation data [31,64].Model simulations in the historical period and the output of BMA were corrected in this step.

Effectiveness Evaluation of Methods
(1) Applicability test of downscaling method.
This study designed an experiment to compare the effectiveness of different interpolation methods.Using CMA historical observations, 8 stations (white dots in Figure 1) were used for validation due to their long and complete series of observations.Observation data from the remaining 139 stations (black dots in Figure 1) were interpolated to the validation stations using three downscaling methods, namely, ANUSPLIN, inverse distance weight (IDW) [65], and original Kriging (OK) [65]; the latter two methods do not take altitude into consideration, which is used to compare with ANUSPLIN.
Here, IDW is used to represent linear two-dimensional interpolation, defined as follows: where Z(x 0 ) represents the interpolation result at verification station x 0 , Z i represents the observation from interpolation station i; n is equal to 139, and d i is the planar Euclidean distance in kilometers between x 0 and station i.
OK is used to represent non-linear two-dimensional interpolation, which is equivalent to thin plate spline interpolation without considering covariates.In this paper, OK is also based on ANUSPLIN; for technical details, please refer to Hutchinson and Gessler [54].
The interpolation results were compared with ground observation data by the Wasserstein distance (WD) [66] and quantile-quantile (Q-Q) plot.WD, also named the earth mover's distance (EMD), is able to measure the distribution similarity between two samples very efficiently [67]; it could be defined as follows: where γ ∼ ∏(P 1 , P 2 ) denotes all possible joint probability distributions between P 1 and P 2 , in f indicates the greatest lower bound, and E (x,y)∼γ [∥x − y∥] denotes the expectation of distance between samples x and y (which follow P 1 and P 2 distributions, respectively, representing the interpolated result and the observation at a verification station) when their joint distribution follows γ.In other words, WD indicates the minimum distance needed to move x into the same distribution as y.The smaller the WD, the closer the probability distribution between the interpolation results and the observations.The interpolation method with the minimum WD means that it most effectively captures the statistical distribution characteristics of the verification station.Due to the inherently high randomness and uncertainty of weather systems, traditional Euclidean metrics, such as the mean absolute error (MAE) and root mean square error (RMSE), are limited to evaluating the effectiveness of interpolation methods [12].This is because almost all interpolation methods struggle to capture the daily variations at unknown points [68].In contrast, evaluating the fitting distribution characteristics through WD is more meaningful [69].Long-term climate change studies may not necessarily require precise daily precipitation or temperature data for specific days, but accurate probability distribution characteristics within specific time periods are essential.This is crucial for research on extreme statistics and is important for both historical simulations and future predictions.Therefore, we adopt the WD as the evaluation metric.
To evaluate the absolute mean bias (AMB) of each model before and after correction, we verified the validity of the correction method.The evaluation was performed as follows: First, we randomly selected 20 years from the historical period of 34 years; GCM outputs from this 20-year period were divided into two equal-length parts, x m_h and x m_p , as seen in Equations ( 5) and (6).Secondly, x m_p was corrected by x m_h and reanalysis data in the corresponding period; the pre-correction and post-correction AMBs against these observations are calculated as follows: where x o denotes the values extracted from the reanalysis data obtained in Section 2.3, and n denotes the sample size in the validation period.The above steps were repeated 10 times in each 0.1 • resolution grid to ensure the performance of the multi-model ensemble forecast.Bias correction was evaluated by the mean value of the results from these 10 iterations to guarantee the robustness of the evaluation [50], and the AMB was taken using absolute values to prevent the positive and negative from canceling out.
Quantile bias (QB) is used as a supplementary test to evaluate the ability of bias correction to extreme values.The QB of the specified percentage p can be defined as follows: where F −1 m−h and F −1 o are the same as in Equations ( 5) and ( 6), and F −1 m−h is estimated by ECDF before and after correction.
Uncertainty is a key factor affecting the usability of model prediction results; the signal-to-noise ratio (SNR) was used to quantify the proportion of valid information amid uncertainty.SNR is defined as follows: with where, x means the prediction result of MMA, x h denotes the historical simulation result, n denotes the number of climate models, and x i denotes the prediction result of model i.It indicates that reliability is higher than uncertainty when SNR is greater than 1.Furthermore, because SNR is dimensionless, it can be used to compare different variables.

Performance of Downscaling Methods in the Historical Period
Table 5 and Figure 3 display the effects of different interpolation methods.ANUSPLIN shows the best-fitting results at all eight validation stations for precipitation as well as the minimum temperatures.For maximum temperature, it also presents the best-fitting results at six validation stations.This suggests that considering altitude as a covariate can significantly reduce statistical distribution errors in downscaling.According to Table 5, ANUSPLIN achieves the smallest average WDs (0.49 mm d −1 for precipitation, 1.18 K for maximum temperature, and 1.39 K for minimum temperature, respectively) among the three downscaling methods for all variables.Only in the cases of maximum temperature at Dingri and Bomi and minimum temperature at Bomi and Maqu do IDW and OK methods perform slightly better than ANUSPLIN, respectively.likely because they do not account for the linear relationship between temperature and elevation.Overall, ANUSPLIN (indicated by purple dots in Figure 3) outperforms the other two methods in addressing the above issues, especially when fitting precipitation extremes, which is crucial for regional hydrological research.Figure 3 demonstrates that ANUSPLIN provides a better fit to the 'x = y' line at most stations.For precipitation, OK consistently overestimates light rain and underestimates heavy rainfall, as indicated by the red crosses in Figure 3. IDW exhibits the same issue at Gaize and Qingshuihe stations, represented by blue prisms in Figure 3a,j.At other stations, IDW tends to overall overestimate precipitation.For temperature, both IDW and OK consistently overestimate temperatures at high-altitude stations (Qingshuihe, Gaize, and Dingri, whose altitudes are higher than 4000 m) and underestimate at relatively lowaltitude stations (Basu, Derong, and Bomi, whose altitudes are lower than 3000 m).This is likely because they do not account for the linear relationship between temperature and elevation.Overall, ANUSPLIN (indicated by purple dots in Figure 3) outperforms the other two methods in addressing the above issues, especially when fitting precipitation extremes, which is crucial for regional hydrological research.
Overall, the results in this section suggest that ANUSPLIN is available for the meteorological element interpolation of the QTP; that is, it is more reasonable than IDW and OK in the performance of precipitation extremes and temperatures.This advantage is likely attributed to its superior ability to model the nearly linear relationship between climate variables and altitude, taking elevation into account as a significant covariate.It also highlights that neglecting altitude's impact in downscaling efforts within the QTP can lead to significant statistical distribution biases, a finding that aligns with prior studies [8,12,[70][71][72].

Multi-Model Ensemble and Bias Correction Performance
The AMBs before and after correction are shown in Figures 4-6, showing a substantial reduction in the AMBs after bias correction across all three variables from the seven GCMs and BMAs for the historical period .The outputs of these three variables from the seven GCMs all showed large biases before correction, as depicted in Figures 4, 5 and 6b-h.For precipitation, the AMBs of seven GCMs ranged from 0.87 to 1.73 mm d −1 , some of them were even higher than the mean precipitation over the QTPs (1.24 mm d −1 approximately); the high bias value (≥3.00 mm d −1 ) occurred in the southeast of the QTP, covering the southern parts of IDR, YLZBR, YZR, MKR, and SWR, before bias correction (Figure 4b-h).For temperature, seven GCM AMBs ranged from 1.35 to 6.84 K (maximum surface temperature, Figure 5b-h) and from 1.23 to 10.63 K (minimum surface temperature, Figure 6b-h).For maximum surface temperature, the high bias value (≥8.50 K) occurred at the edge of QDB and most of TMR, IDR, and YLZBR; for minimum surface temperature, the high bias value (≥8.50 K) occurred in the western part of the QTP (TMR, IDR, and western parts of IQTB and YLZBR) and small parts of QDB, before bias correction.Therefore, such uncorrected data cannot be used for climate research and decision-making at the regional scale [31].After bias correction, the GCMs' absolute mean bias was reduced to an acceptable range.For precipitation, it was from 0.13 to 0.19 mm d −1 (about an 85% to 89% reduction Figure 4j-p).For temperature, they were from 0.40 to 0.64 K (maximum surface temperature, Figure 5j-p) and from 0.37 to 0.55 K (minimum surface temperature, Figure 6j-p).High values only occurred in smaller areas compared with uncorrected data.For precipitation, only NorESM2-MM exhibited high AMB (≥3.50 mm d −1 ) in a small part of LYZV (Figure 4k).The reduction in bias for mean annual precipitation, as well as the maximum and minimum temperatures, demonstrates the effectiveness of our bias correction based on CNCDFm.To verify the correctness of the corrected data, we separately analyzed the seasonal cycles before and after bias correction (Figure 7) and the mean annual values after bias correction (Figures S1-S3) for three variables over the historical period.The data show a wet bias for precipitation and a relatively cold bias for temperature in the GCMs' seasonal cycles compared to observations from before bias correction (Figure 7a-c), typically due to systematic errors overlooking the monsoon effect and surface meteorological processes in GCMs [27,28,73].Compared to the uncorrected data, the seasonal cycles of the corrected data more closely align with observations (Figure 7d-f S3 show the spatial distributions of mean annual values for the corrected data, indicating a strong agreement with observations.Precipitation data, in particular, show a consistent decline from the southeast to the northwest of the QTP (Figure S1c-i), aligning with prior studies [7,8,74].Regarding temperature, all seven GCMs indicate lower values in the IQTB due to its high altitudes and latitudes, with higher temperatures primarily found in the southeastern part (Figures S2c-i and S3c-i), which is consistent with orographic findings [51].Overall, our results indicate that the CNCDFm method has produced a reasonably corrected dataset based on CMIP6 GCMs, which is suitable for climate impact studies in the QTP. ) between individual models in CMIP6 and observational data over the QTP, AMB before correction on the left of the blue dotted line (a-h), AMB after correction on the right side (i-p). Te mean value of AMB across the entire QTP is denoted in each subgraph as m; m is defined as the sum of all grid point values divided by the total number of grid points over the Tibetan Plateau, as shown in the following figure.The numbering (I-X) of each subregion (consistent with Figure 1) is also indicated in each subfigure.
To verify the correctness of the corrected data, we separately analyzed the seasonal cycles before and after bias correction (Figure 7) and the mean annual values after bias correction (Figures S1-S3) for three variables over the historical period.The data show a wet bias for precipitation and a relatively cold bias for temperature in the GCMs' seasonal cycles compared to observations from before bias correction (Figure 7a-c), typically due to systematic errors overlooking the monsoon effect and surface meteorological processes in GCMs [27,28,73].Compared to the uncorrected data, the seasonal cycles of the corrected data more closely align with observations (Figure 7d-f), as CNCDFm has corrected biases for both high and low values [61].The uncertainty in GCM outputs has been significantly reduced post-bias correction, as evidenced by the narrowed ranges (shaded areas in Figure 7d-f).Figures S1-S3 show the spatial distributions of mean annual values for the corrected data, indicating a strong agreement with observations.Precipitation data, in particular, show a consistent decline from the southeast to the northwest of the QTP (Figure S1c-i), aligning with prior studies [7,8,74].Regarding temperature, all seven GCMs indicate lower values in the IQTB due to its high altitudes and latitudes, with higher temperatures primarily found in the southeastern part (Figures S2c-i and S3c-i), which is consistent with orographic findings [51].Overall, our results indicate that the CNCDFm method has produced a reasonably corrected dataset based on CMIP6 GCMs, which is suitable for climate impact studies in the QTP.At the same time, we calculate the mean biases, seasonal cycles, and mean annual values for the BMA (as shown in Figures 4-6a,i, the blue line in Figures 7 and S1-S3b).Before the bias correction, BMA exhibited the lowest absolute mean bias (AMB) compared to the outputs from the seven GCMs, attributed to the consideration of the models' posterior probability distributions (Figures 4-6a).After bias correction, BMA showed the lowest  At the same time, we calculate the mean biases, seasonal cycles, and mean annual values for the BMA (as shown in Figures 4-6a,i, the blue line in Figures 7 and S1-S3b).Before the bias correction, BMA exhibited the lowest absolute mean bias (AMB) compared to the outputs from the seven GCMs, attributed to the consideration of the models' posterior probability distributions (Figures 4-6a).After bias correction, BMA showed the lowest the performance of bias correction in extremes, the QBs for uncorrected and corrected BMA 95th percentiles were compared against observations for the historical period (1979-2013) for three variables, as shown in Table 1 (Figure S4).The analysis revealed that biases in the 95th percentiles were successfully mitigated, with the corrected BMA data bias approaching 0 (Figure S4b,c,f) compared with the uncorrected.Thus, the corrected data also demonstrate improved performance for extreme values.

Future Climate Change Projection in the QTP
Projected changes in precipitation and maximum/minimum surface temperatures in the QTP for the 21st century, based on outputs from the corrected BMA and seven GCMs, span from 2025 to 2100.We estimated these changes for the near (2025-2049), medium (2050-2074), and far (2075-2100) future periods, using a historical baseline of 1985-2013, across all four scenarios (SSP1-2.6,SSP2-4.5,SSP3-7.0, and SSP5-8.5),as shown in Figures 8-10.Overall, the projections indicate consistent warming and wetting trends for the QTP At the same time, we calculate the mean biases, seasonal cycles, and mean annual values for the BMA (as shown in Figures 4, 5 and 6a,i, the blue line in Figures 7 and S1-S3b).Before the bias correction, BMA exhibited the lowest absolute mean bias (AMB) compared to the outputs from the seven GCMs, attributed to the consideration of the models' posterior probability distributions (Figures 4, 5 and 6a).After bias correction, BMA showed the lowest AMB for precipitation (Figure 4i) and avoided the high-temperature AMB that occurred in every individual model (Figures 5i and 6i).In terms of seasonal cycles, the corrected BMA data align well overall with observations after bias correction (Figure 7d-f).Furthermore, the spatial distribution of the corrected BMA output is consistent with the observations (Figures S1-S3b).Overall, BMA provides a robust multi-model ensemble forecast result for precipitation, as well as maximum and minimum temperatures across the Tibetan Plateau.
Studying extreme events is a crucial aspect of climate change research.To evaluate the performance of bias correction in extremes, the QBs for uncorrected and corrected BMA 95th percentiles were compared against observations for the historical period (1979-2013) for three variables, as shown in Table 1 (Figure S4).The analysis revealed that biases in the 95th percentiles were successfully mitigated, with the corrected BMA data bias approaching 0 (Figure S4b,c,f) compared with the uncorrected.Thus, the corrected data also demonstrate improved performance for extreme values.

Future Climate Change Projection in the QTP
Projected changes in precipitation and maximum/minimum surface temperatures in the QTP for the 21st century, based on outputs from the corrected BMA and seven GCMs, span from 2025 to 2100.We estimated these changes for the near (2025-2049), medium (2050-2074), and far (2075-2100) future periods, using a historical baseline of 1985-2013, across all four scenarios (SSP1-2.6,SSP2-4.5,SSP3-7.0, and SSP5-8.5),as shown in Figures 8-10.Overall, the projections indicate consistent warming and wetting trends for the QTP in the 21st century, aligning with historical trends [6,75], except for SSP1-2.6,where the rate of increase in precipitation and temperature is projected to slow down significantly in the far future (Figures 8, 9 and 10c).Compared to the low radiative forcing pathway (SSP1-2.6),high radiative forcing pathways (SSP3-7.0 and SSP5-8.5)have considerably large increases in precipitation and temperature.For instance, under SSP5-8.5,mean precipitation is projected to rise by 26.8% in the far future (Figure 9i), with mean maximum and minimum surface temperature increases of 5.1 • C and 5.7 • C, respectively (Figures 9 and 10i).Excessive warming may accelerate the retreat of plateau glaciers, especially in the TMR and IQTB regions.Furthermore, it shows a consistent spatial pattern in the different SSPs.For instance, a higher precipitation increase was found in the northern parts of the IQTB and IDR, TMR, and QDB under SSP3-7.0 and SSP5-8.5 in the far future.
Atmosphere 2024, 15, 434 16 of in the 21st century, aligning with historical trends [6,75], except for SSP1-2.6,where t rate of increase in precipitation and temperature is projected to slow down significant in the far future (Figures 8-10c).Compared to the low radiative forcing pathway (SSP 2.6), high radiative forcing pathways (SSP3-7.0 and SSP5-8.5)have considerably large i creases in precipitation and temperature.For instance, under SSP5-8.5,mean precipitatio is projected to rise by 26.8% in the far future (Figure 9i), with mean maximum and min mum surface temperature increases of 5.1 °C and 5.7 °C, respectively (Figures 9 and 10 Excessive warming may accelerate the retreat of plateau glaciers, especially in the TM and IQTB regions.Furthermore, it shows a consistent spatial pattern in the different SSP For instance, a higher precipitation increase was found in the northern parts of the IQT and IDR, TMR, and QDB under SSP3-7.0 and SSP5-8.5 in the far future.Future precipitation and temperature were also estimated in 10 subregions based o corrected BMA (Tables S1-S3), as seen in the supplementary materials, Figures 8-10.T projections indicate increases in the mean annual precipitation and maximum and min mum temperatures across these subregions under middle and high radiative forcing pat ways (SSP2-4.5,SSP3-7.0, and SSP5-8.5),with significantly higher increases across hig radiative forcing scenarios compared to low radiative forcing ones, aligning wi  (d-f), (g-i) and (j-l), the same as (a-c) but for SSP2-4.5,SSP3-7.0, and SSP5-8.5, respectively.The mean change across the entire QTP (m) is shown in each panel as m.The numbering of each subregion (consistent with Figure 1) is also indicated in each subfigure.The numbering (I-X) of each subregion (consistent with Figure 1) is also indicated in each subfigure.
projections for South Asia [31].Notably, under SSP3-7.0 and SSP5-8.5, all three variabl are projected to continue to rise, which is consistent with the overall trend across the QT Under SSP1-2.6, the precipitation trend in HXB, IDR, IQTB, QDB, and TMR in the far f ture is lower than in the middle future, a trend likely attributed to the mitigating effec of reduced greenhouse gas concentrations on climate change.Overall, the bias-correct dataset produced a reasonable projection for subregions in the QTP, and it is available f hydroclimatic research at the regional scale.Figure 11 illustrates the time series for BMA (represented by solid lines) and MM (represented by dotted lines), demonstrating a more continuous growth trend in tempe ature (both maximum and minimum) in SSP2-4.5,SSP3-7.0, and SSP5-8.5, as observed a ter applying a 30-year moving average.Conversely, a distinct turning point occurs mi century in SSP1-2.6,where warming shifts to cooling across all subregions, highlightin the impact of carbon emission limitations on regional temperature trends.This warmi trajectory parallels the global warming process projected by the IPCC's Sixth Assessme Report (AR6), which predicts a global surface temperature increase of 4 °C by the ce tury's end under SSP5-8.5 compared to 1995-2014 [76].Notably, the rate of warming the QTP exceeds the global average, underscoring the region's sensitivity to clima change, as historically observed [6].
Furthermore, noteworthy changes in the spatial distribution of precipitation an temperature across the Tibetan Plateau merit attention.Specifically, future projections i dicate a significant increase in precipitation rates in IDR, QDB, TMR, and the northwe of IQTB (Figures 8 and 11(b1,c1,e1,g1), Table S1), which is not only far ahead of glob Future precipitation and temperature were also estimated in 10 subregions based on corrected BMA (Tables S1-S3), as seen in the Supplementary Materials, Figures 8-10.
The projections indicate increases in the mean annual precipitation and maximum and minimum temperatures across these subregions under middle and high radiative forcing pathways (SSP2-4.5,SSP3-7.0, and SSP5-8.5),with significantly higher increases across high radiative forcing scenarios compared to low radiative forcing ones, aligning with projections for South Asia [31].Notably, under SSP3-7.0 and SSP5-8.5, all three variables are projected to continue to rise, which is consistent with the overall trend across the QTP.Under SSP1-2.6, the precipitation trend in HXB, IDR, IQTB, QDB, and TMR in the far future is lower than in the middle future, a trend likely attributed to the mitigating effects of reduced greenhouse gas concentrations on climate change.Overall, the bias-corrected dataset produced a reasonable projection for subregions in the QTP, and it is available for hydroclimatic research at the regional scale.
Figure 11 illustrates the time series for BMA (represented by solid lines) and MMA (represented by dotted lines), demonstrating a more continuous growth trend in temperature (both maximum and minimum) in SSP2-4.5,SSP3-7.0, and SSP5-8.5, as observed after applying a 30-year moving average.Conversely, a distinct turning point occurs midcentury in SSP1-2.6,where warming shifts to cooling across all subregions, highlighting the impact of carbon emission limitations on regional temperature trends.This warming trajectory parallels the global warming process projected by the IPCC's Sixth Assessment Report (AR6), which predicts a global surface temperature increase of 4 • C by the century's end under SSP5-8.5 compared to 1995-2014 [76].Notably, the rate of warming in the QTP exceeds the global average, underscoring the region's sensitivity to climate change, as historically observed [6].
drastic change in precipitation is bound to have a profound impact on the regional hydro logical cycle.For temperature, there is no obvious spatial differentiation in the near future change (Figures 9 and 10a,d,g,j), while in the far future, spatial differentiation is promi nent.For instance, IQTB and TMR have shown higher warming than other regions in SSP5-8.5, which may be due to the underlying positive feedback of the radiation absorp tion and will accelerate the ablation of regional glacial permafrost [81][82][83][84].These spatia differences highlight the necessity of high-resolution projection for research in the QTP.Table S4 provides statistics on the mean rate and change rate of the annual maximum five-day precipitation for future years based on BMA.Notably, IDR exhibits the highes change rate among the 10 sub-regions in the vast majority of scenarios and periods, with growth rates reaching 22.9 to 41.2% in the distant future.The significant increase in the annual maximum five-day precipitation, as an indicator of extreme precipitation, implies an elevated risk of hydrogeological disasters, such as floods, landslides, and debris flow [85][86][87].Most areas in IDR are above 4000 m, given its higher elevation (Figure 1).Geolog ical disasters in this region will likely cause more significant damage downstream [88] highlighting the need for enhanced observation and monitoring in this area.Furthermore, noteworthy changes in the spatial distribution of precipitation and temperature across the Tibetan Plateau merit attention.Specifically, future projections indicate a significant increase in precipitation rates in IDR, QDB, TMR, and the northwest of IQTB (Figures 8 and 11(b1,c1,e1,g1), Table S1), which is not only far ahead of global future growth but also far ahead of the rest of the QTP.The future moistening trend in the northern region is consistent with historical trends.Zhang et al. estimated a historical precipitation growth rate of 1.5-3.0mm yr −1 in this area [77].The significant growth rate may be partly due to the historically low precipitation in this area (Table S1), and partly due to the enhanced atmospheric circulation caused by warming [78], which accelerates the local moisture recycling and increases the input of external moisture [77,79,80].However, the drastic change in precipitation is bound to have a profound impact on the regional hydrological cycle.For temperature, there is no obvious spatial differentiation in the near future change (Figures 9 and 10a,d,g,j), while in the far future, spatial differentiation is prominent.For instance, IQTB and TMR have shown higher warming than other regions in SSP5-8.5, which may be due to the underlying positive feedback of the radiation absorption and will accelerate the ablation of regional glacial permafrost [81][82][83][84].These spatial differences highlight the necessity of high-resolution projection for research in the QTP.
Table S4 provides statistics on the mean rate and change rate of the annual maximum five-day precipitation for future years based on BMA.Notably, IDR exhibits the highest change rate among the 10 sub-regions in the vast majority of scenarios and periods, with growth rates reaching 22.9 to 41.2% in the distant future.The significant increase in the annual maximum five-day precipitation, as an indicator of extreme precipitation, implies an elevated risk of hydrogeological disasters, such as floods, landslides, and debris flow [85][86][87].Most areas in IDR are above 4000 m, given its higher elevation (Figure 1).Geological disasters in this region will likely cause more significant damage downstream [88], highlighting the need for enhanced observation and monitoring in this area.
Atmosphere 2024, 15, 434 20 of 31 Geological disasters in this region will likely cause more significant damage downstream [88], highlighting the need for enhanced observation and monitoring in this area.

Discussion
In the present study, the effectiveness of downscaling and bias correction methods were evaluated.The projection results are presented at a basin scale and compared with previous research.Such results could verify the rationality of HRP-QTP, but some issues remain that deserve further discussion.

Availability of Reanalysis Precipitation
In numerous areas worldwide, the scarcity of observation data significantly hampers the execution of various environmental research initiatives, especially for precipitation [89].He et al. [11] evaluated the accuracy of CMFD in simulating conditions in regions across western China that are poorly monitored, using additional data from measuring

Discussion
In the present study, the effectiveness of downscaling and bias correction methods were evaluated.The projection results are presented at a basin scale and compared with previous research.Such results could verify the rationality of HRP-QTP, but some issues remain that deserve further discussion.

Availability of Reanalysis Precipitation
In numerous areas worldwide, the scarcity of observation data significantly hampers the execution of various environmental research initiatives, especially for precipitation [89].He et al. [11] evaluated the accuracy of CMFD in simulating conditions in regions across western China that are poorly monitored, using additional data from measuring stations in IDR, IQTB, YZR, and SWR.The findings indicate that CMFD provides a more accurate representation of temperature in areas of the QTP that lack direct observational data, leading us to adopt CMFD temperature data as a proxy for historical observations.However, the dataset's precipitation measurements have not been similarly validated in the QTP.
Through the extension and correction of GPQM and CMFD data, we derived a combined dataset aimed at correcting the CMFD's underestimation of precipitation in the LYZV area, using additional station information from GPQM beyond the QTP.The historical annual mean precipitation across 12 subregions of the QTP (as shown in Figure 12) was estimated based on the data generated in Section 2.3.The indicate that the QTP's historical annual mean precipitation is 456 mm/year, with a slight increase over historical periods, aligning closely with ground station measurements [75].Spatial analysis reveals that the eastern and southern subregions, such as MKR, SWR, YZR, and YR, exhibit higher precipitation levels (602 mm/year, 653 mm/year, 598 mm/year, and 500 mm/year, respectively) compared to the central and northern subregions (HXB, IQTB, QTB, and TMR, with annual mean precipitation of 241 mm/year, 261 mm/year, 168 mm/year, and 145 mm/year, respectively), reflecting the patterns of water vapor transport in the QTP [8,47].Notably, due to its high precipitation, at approximately 3120 mm/year, the YLZBR region records the highest precipitation among all subregions at 1108 mm/year.
Atmosphere 2024, 15, 434 20 of 30 However, the dataset's precipitation measurements have not been similarly validated in the QTP.Through the extension and correction of GPQM and CMFD data, we derived a combined dataset aimed at correcting the CMFD's underestimation of precipitation in the LYZV area, using additional station information from GPQM beyond the QTP.The historical annual mean precipitation across 12 subregions of the QTP (as shown in Figure 12) was estimated based on the data generated in Section 2.3.The results indicate that the QTP's historical annual mean precipitation is 456 mm/year, with a slight increase over historical periods, aligning closely with ground station measurements [75].Spatial analysis reveals that the eastern and southern subregions, such as MKR, SWR, YZR, and YR, exhibit higher precipitation levels (602 mm/year, 653 mm/year, 598 mm/year, and 500 mm/year, respectively) compared to the central and northern subregions (HXB, IQTB, QTB, and TMR, with annual mean precipitation of 241 mm/year, 261 mm/year, 168 mm/year, and 145 mm/year, respectively), reflecting the patterns of water vapor transport in the QTP [8,47].Notably, due to its high precipitation, at approximately 3120 mm/year, the YLZBR region records the highest precipitation among all subregions at 1108 mm/year.Generally, gauge data are considered the most accurate precipitation data [90].However, it is difficult to estimate and validate LYZV precipitation due to a lack of observations (Figure 1), and reanalysis based on multi-source data is the only method currently available to estimate precipitation in LYZV.Due to differences in methods and data sources, different reanalysis datasets have demonstrated a wide range of estimated results, ranging from 1000 mm/year (for instance, He et al. [11] estimated the LYZV annual mean precipitation at 1001 mm/year by combining upstream ground station data with satellite reanalysis data) to 4000 mm/year (the estimation by Ma et al. [12] is 3011 mm/year using interpolation of downstream gauge data; the estimation by Li et al. [8] is from 1000 to 4000 mm/year at different altitudes using grid data from the Global Precipitation Measurement, GPM).We tend to use the higher precipitation estimates for LYZV because they are more consistent with precipitation characteristics in the southern Himalayas, which are close to the same altitude regions in Bhutan [45] and corroborate early findings by Yang et al. [91].Nonetheless, given the QTP's extensive area and complex topography, different reanalysis datasets offer distinct advantages across various regions [20].Accurately estimating precipitation across the QTP remains an unresolved challenge that requires further observational evidence.Generally, gauge data are considered the most accurate precipitation data [90].However, it is difficult to estimate and validate LYZV precipitation due to a lack of observations (Figure 1), and reanalysis based on multi-source data is the only method currently available to estimate precipitation in LYZV.Due to differences in methods and data sources, different reanalysis datasets have demonstrated a wide range of estimated results, ranging from 1000 mm/year (for instance, He et al. [11] estimated the LYZV annual mean precipitation at 1001 mm/year by combining upstream ground station data with satellite reanalysis data) to 4000 mm/year (the estimation by Ma et al. [12] is 3011 mm/year using interpolation of downstream gauge data; the estimation by Li et al. [8] is from 1000 to 4000 mm/year at different altitudes using grid data from the Global Precipitation Measurement, GPM).We tend to use the higher precipitation estimates for LYZV because they are more consistent with precipitation characteristics in the southern Himalayas, which are close to the same altitude regions in Bhutan [45] and corroborate early findings by Yang et al. [91].Nonetheless, given the QTP's extensive area and complex topography, different reanalysis datasets offer distinct advantages across various regions [20].Accurately estimating precipitation across the QTP remains an unresolved challenge that requires further observational evidence.

Availability of Downscaling Methods
In Section 3.1, we present a simple comparative experiment to demonstrate the effectiveness of ANUSPLIN, which takes elevation into account, in downscaling climatic elements in the QTP.Herein, IDW represents linear two-dimensional interpolation, and OK represents non-linear two-dimensional interpolation, both of which are compared with the non-linear three-dimensional interpolation of ANUSPLIN.
The interpolation result of IDW equals the weighted average of all known stations, with the weights being inversely proportional to the distance between the known points and the interpolation point.Therefore, once the locations of the sites are determined, the interpolation relationship is also established.OK and ANUSPLIN perform interpolation by constructing non-linear relationships to fit the climate surface of known stations [54].The fitting relationship is time-varying (in this study, both use annual fitting).The difference lies in that ANUSPLIN considers covariates.
Compared to ANUSPLIN, the limitation of two-dimensional interpolations lies in their focus on spatial relationships between stations, leading to a reliance on information from nearby stations and overlooking the decisive role of elevation on plateau climates.Taking temperature as an example, a common rule is that temperature decreases with the increase in elevation [7].The elevation of stations near the Basu station is higher than that of Basu (Figure 1), with the nearest five having an average elevation of 3082 m, significantly higher than Basu's 2589.2 m.This results in the two-dimensional interpolations significantly underestimating the temperature at Basu (Figure 3e,f).This bias is corrected in ANUSPLIN.Furthermore, the fitting performed by ANUSPLIN is time-varying, meaning that if the correlation between climate elements and elevation changes in future scenarios, the new correlation can also be captured [92].This reduces the systematic bias caused by limitations in fitting capability during downscaling.
However, downscaling climate models continues to be an exceedingly complex issue.Identifying the optimal downscaling method and understanding how systematic errors are transmitted during the downscaling process urgently require further investigation.This study merely demonstrates the viability of applying ANUSPLIN in plateau areas.

Uncertainty and Reliability in Future Projections
Previous studies have indicated significant divergence in climate projections for the QTP from GCMs [27,28,30], leading to low SNRs.As shown in Figures 13-15, bias correction significantly increases the SNR for three variables: for precipitation, the SNR range increased from 0.23-0.40 to 0.43-0.90;for temperature, it rose from 0.37-1.12 to 3.67-5.52.This suggests that bias correction effectively eliminates systemic model biases, making future model projections more reliable.However, reliability varies greatly among different variables.For temperature, both maximum and minimum, the SNR exceeds 1 in all subregions (Figures 14 and 15) after correction, indicating a clear trend in the future temperature rise and reduced uncertainty.In contrast, precipitation shows an overall SNR below 1, with several regions, including LYZV, MKR, SWR, YZR, and YR (Figure 14b,d,f,h), exhibiting SNRs far below 1.This suggests that future changes in precipitation for these regions are not pronounced and are subject to considerable uncertainty, a conclusion supported by Figure 11(d1,f1,h1,i1,j1).
We compare our findings with those from Cui et al. [30], who projected precipitation for several major river basins using a broader set of GCMs, albeit at a lower spatial resolution.Their results, showing SNRs in Table 6, reveal that SNRs exceeding 1 only occur in the Indus, Mekong, Salween, and Yellow River basins under SSP2-4.5 when future temperatures increase by 3 • C. Contrasting these findings with ours (Figures 13-15) indicates that the SNRs for multi-model projections remain consistently low, suggesting that merely increasing the number of predictive models does not necessarily enhance prediction reliability.Mechanistically, precipitation predictions depend on factors such as wind, humidity, and aerosols [93][94][95].The approximation and parameterization of these factors in GCMs might introduce more uncertainty into precipitation forecasts than temperature forecasts.Furthermore, while current studies often prioritize MMA, we advocate for BMA based on our belief that historical simulations can inform and constrain future climate predictions.Unlike BMA, MMA does not incorporate prior information from historical periods, which is particularly relevant for temperature projections where the differences between BMA and MMA are minimal.However, for precipitation, especially in regions like IQTB, MKR, and YZR post-2050, BMA projections are significantly higher than those from MMA, as illustrated in Figure 11.Precipitation, subject to considerable model uncertainty, poses challenges.Ferguglia et al. [96] demonstrated that the emergent constraint approach-relying on a physically explainable link between historical and future data-does not robustly address systematic uncertainties, indicating the limitations of depending solely on multimodel precipitation outcomes.However, due to the lack of historical data and the unclear precipitation mechanism, statistical correction is still a widely used reanalysis method in existing research.We retained the correction results of each model in the final forecast dataset provided for further study.Furthermore, while current studies often prioritize MMA, we advocate for BMA based on our belief that historical simulations can inform and constrain future climate predictions.Unlike BMA, MMA does not incorporate prior information from historical periods, which is particularly relevant for temperature projections where the differences between BMA and MMA are minimal.However, for precipitation, especially in regions like IQTB, MKR, and YZR post-2050, BMA projections are significantly higher than those from MMA, as illustrated in Figure 11.Precipitation, subject to considerable model uncertainty, poses challenges.Ferguglia et al. [96] demonstrated that the emergent constraint approach-relying on a physically explainable link between historical and future data-   Table 6.SNRs for precipitation projection results from Cui et al. [30].The bold terms indicate SNRs greater than 1; SNRs were calculated using Equation ( 11) ignoring models with no data, and the division of subregions is slightly different from our study.

Name Scenarios
Temperature Rise Range 1.5

Conclusions
The QTP serves as Asia's most crucial water source, with climate change in the QTP poised to profoundly affect water security and geopolitical relations across the continent.

Conclusions
The QTP serves as Asia's most crucial water source, with climate change in the QTP poised to profoundly affect water security and geopolitical relations across the continent.The region's complex terrain necessitates high-resolution climate projections to underpin studies on regional climate change and guide adaptation strategies.To address this need, our research developed a high-resolution projection by utilizing reanalyzed meteorological datasets alongside outputs from seven GCMs.We employed a downscaling method that accounts for altitude variations and implemented an updated nonstationary bias correction technique, both of which were rigorously tested for validity.Key findings from our study include the following: 1.
There is an approximately linear relationship between climatic factors and altitude in high-altitude areas.Therefore, the ANUSPLIN software, which incorporates DEM as a covariate, can more accurately fit the probability distributions of meteorological elements.Ignoring altitude's impact on the distribution of meteorological elements during downscaling leads to the overestimation or underestimation of temperatures and overlooks extreme precipitation events.Thus, incorporating elevation data as prior information in downscaling processes is crucial.Despite significant systematic biases in GCM simulations and predictions for the QTP climate, these can be substantially mitigated by employing an advanced nonstationary bias correction method, CNCDFm.For example, this method has been shown to minimize the deviation of extreme values to nearly zero.Consequently, our research enhances the applicability and accuracy of GCM results for regional studies in the QTP.

2.
Based on high-resolution projections, the QTP is expected to continue experiencing warming and wetting trends.The projected temperature change range remains consistent across the QTP subregions, mirroring the historical trend where the QTP temperature increase surpasses the global surface temperature rise.Under the Shared Socioeconomic Pathway SSP1-2.6, the region is anticipated to reach a warming inflection point by the mid-21st century.Furthermore, future precipitation projections indicate significant spatial variations, with growth rates in regions such as IDR, IQTB, QDB, and TMR far exceeding the global average.Such marked increases in temperature and precipitation are likely to fundamentally transform the regional hydrological cycle.

3.
Two significant challenges hinder the accurate forecasting of future precipitation trends in the QTP.Firstly, the historical period lacks global high-resolution precipitation data that are convincing enough for detailed analysis, especially in ungauged regions such as the IQTB, YLZBR, TMR, and IDR.This scarcity complicates the analysis of historical precipitation changes and the calibration of GCMs.Secondly, there remains substantial uncertainty in GCM projections.While bias correction techniques can mitigate some systematic biases inherent in GCMs, the reliance on historical data to constrain future precipitation projections proves to be less effective than for temperature, rendering precipitation forecasts less reliable.This uncertainty also leads to significant discrepancies between the precipitation forecasts of BMA and MMA, both of which indicate an increasing trend across all basins.Further research is imperative to enhance the reliability of multi-model precipitation predictions.
Despite the limitations outlined above, this research provides a high-resolution future prediction for regional-scale studies of the QTP that is openly accessible for further testing, research, and decision-making purposes.

Figure 1 .
Figure 1.Location and geographical domains for HF-QTP.Black dots and white dots denote the meteorological stations from CMA, white dots are used as verification stations, and black dots are used as interpolation stations in Section 2.6.The QTP's major river basin boundaries are in black and the main streams are in blue.Due to the absence of surveying and mapping data, the rivers within the Qiangtang region are not depicted on the map.Names of the major river basin boundaries are written within each basin boundary.DEM is shown as the color scale in the background.The boundary of the low Yarlung Zangbo Valley (LYZV) is shown with a red dotted line.

Figure 1 .
Figure 1.Location and geographical domains for HF-QTP.Black dots and white dots denote the meteorological stations from CMA, white dots are used as verification stations, and black dots are used as interpolation stations in Section 2.6.The QTP's major river basin boundaries are in black and the main streams are in blue.Due to the absence of surveying and mapping data, the rivers within the Qiangtang region are not depicted on the map.Names of the major river basin boundaries are written within each basin boundary.DEM is shown as the color scale in the background.The boundary of the low Yarlung Zangbo Valley (LYZV) is shown with a red dotted line.

Figure 2 .
Figure 2. Flowchart for producing historical precipitation data over the QTP based on GPQM and CMFD.

Figure 2 .
Figure 2. Flowchart for producing historical precipitation data over the QTP based on GPQM and CMFD.

Figure 3 .
Figure 3.The Q-Q plot between downscaling results and observations from CMA in validation stations.The subfigures are numbered in sequence from (a-x) and labeled in the top left corner of each subfigure according to their order.Different colors and dot types denote different downscaling methods.Horizontal axes denote observation values and vertical axes denote downscaling results; they have the same units.The gray dotted line, "x = y", indicates the 1:1 line.Values closer to this line indicate a more accurate statistical distribution.

Atmosphere 2024 ,
15, 434 13 of 30 maximum and minimum temperatures, demonstrates the effectiveness of our bias correction based on CNCDFm.

Figure 4 .
Figure 4. Absolute mean bias (AMB) for precipitation (mm d −1) between individual models in CMIP6 and observational data over the QTP, AMB before correction on the left of the blue dotted line (a-h), AMB after correction on the right side (i-p).The mean value of AMB across the entire QTP is denoted in each subgraph as m; m is defined as the sum of all grid point values divided by the total number of grid points over the Tibetan Plateau, as shown in the following figure.The numbering (I-X) of each subregion (consistent with Figure1) is also indicated in each subfigure.
), as CNCDFm has corrected biases for both high and low values [61].The uncertainty in GCM outputs has been significantly reduced post-bias correction, as evidenced by the narrowed ranges (shaded areas in Figure 7d-f).Figures S1-

Figure 4 .
Figure 4. Absolute mean bias (AMB) for precipitation (mm d −1) between individual models in CMIP6 and observational data over the QTP, AMB before correction on the left of the blue dotted line (a-h), AMB after correction on the right side (i-p).The mean value of AMB across the entire QTP is denoted in each subgraph as m; m is defined as the sum of all grid point values divided by the total number of grid points over the Tibetan Plateau, as shown in the following figure.The numbering (I-X) of each subregion (consistent with Figure1) is also indicated in each subfigure.

Atmosphere 2024, 15 , 434 14 of 30 Figure 5 .
Figure 5. AMB is the same as in Figure 4 but for the maximum surface temperature (K).

Figure 6 .
Figure 6.AMB is the same as in Figure4but for minimum surface temperature (K).

Figure 5 . 30 Figure 5 .
Figure 5. AMB is the same as in Figure 4 but for the maximum surface temperature (K).

Figure 6 .
Figure 6.AMB is the same as in Figure4but for minimum surface temperature (K).

Figure 6 .
Figure 6.AMB is the same as in Figure4but for minimum surface temperature (K).

Figure 7 .
Figure 7.The seasonal cycle of precipitation, as well as maximum and minimum temperature, before ((a-c), on the left side) and after ((d-f), on the right side) bias correction.Comparison of the Bayesian model averaging (BMA, blue) mean seasonal cycle of bias-corrected (a,d) precipitation, (b,e) maximum temperature, and (c,f) minimum temperature against the observations for the 1979-2013 period (red).The shaded area represents the uncertainty (from minimum to maximum) of all 7 CMIP6 GCMs.

Figure 7 .
Figure 7.The seasonal cycle of precipitation, as well as maximum and minimum temperature, before ((a-c), on the left side) and after ((d-f), on the right side) bias correction.Comparison of the Bayesian model averaging (BMA, blue) mean seasonal cycle of bias-corrected (a,d) precipitation, (b,e) maximum temperature, and (c,f) minimum temperature against the observations for the 1979-2013 period (red).The shaded area represents the uncertainty (from minimum to maximum) of all 7 CMIP6 GCMs.

Figure 8 .
Figure 8. Future change of precipitation (%) against historical baseline (1985-2010) based on co rected BMA.(a-c) The projected change in mean annual precipitation for the SSP1-2.6 in the ne (2025-2049), medium (2050-2074), and far (2075-2100) future; (d-f), (g-i) and (j-l), the same as ( c) but for SSP2-4.5,SSP3-7.0, and SSP5-8.5, respectively.The mean change across the entire QTP ( is shown in each panel as m.The numbering of each subregion (consistent with Figure 1) is al indicated in each subfigure.The numbering (I-X) of each subregion (consistent with Figure 1) is al indicated in each subfigure.

Figure 8 .
Figure 8. Future change of precipitation (%) against historical baseline (1985-2010) based on corrected BMA.(a-c) The projected change in mean annual precipitation for the SSP1-2.6 in the near (2025-2049), medium (2050-2074), and far (2075-2100) future; (d-f),(g-i) and (j-l), the same as (a-c) but for SSP2-4.5,SSP3-7.0, and SSP5-8.5, respectively.The mean change across the entire QTP (m) is shown in each panel as m.The numbering of each subregion (consistent with Figure1) is also indicated in each subfigure.The numbering (I-X) of each subregion (consistent with Figure1) is also indicated in each subfigure.

Figure 9 .
Figure 9. Future changes based on BMA; the same as in Figure 8 but for the maximum surface temperature.

Figure 9 .
Figure 9. Future changes based on BMA; the same as in Figure 8 but for the maximum surface temperature.

Figure 10 .
Figure 10.Future changes based on BMA; the same as in Figure 8 but for the minimum surface temperature.

Figure 10 .
Figure 10.Future changes based on BMA; the same as in Figure 8 but for the minimum surface temperature.

Figure 11 .
Figure 11.Future change time series of BMA and multi-model average (MMA) in precipitation, as well as maximum and minimum temperatures in different subregions (Numbered as a-j, marked in the top left corner of each subfigure) of the QTP.Changes in the BMA and MMA for mean annual precipitation (%), maximum temperature (°C), and minimum temperature (°C) were estimated using a 30-year moving window against the historical reference period of 1979-2013.The shaded region shows uncertainty (estimated using one standard deviation) based on 7 GCMs.Different colors denote the 4 SSPs, and different types of lines denote BMA and MMA, respectively.

Figure 11 .
Figure 11.Future change time series of BMA and multi-model average (MMA) in precipitation, as well as maximum and minimum temperatures in different subregions (Numbered as (a-j), marked in the top left corner of each subfigure) of the QTP.Changes in the BMA and MMA for mean annual precipitation (%), maximum temperature ( • C), and minimum temperature ( • C) were estimated using a 30-year moving window against the historical reference period of 1979-2013.The shaded region shows uncertainty (estimated using one standard deviation) based on 7 GCMs.Different colors denote the 4 SSPs, and different types of lines denote BMA and MMA, respectively.

Figure 12 .
Figure 12.Annual mean from 1979 to 2013 (a) and seasonal cycle mean (b) of precipitation from combined data in different subregions of the QTP; the black lines represent the mean values of the QTP, and the colored lines denote multiple subregions.

Figure 12 .
Figure 12.Annual mean from 1979 to 2013 (a) and seasonal cycle mean (b) of precipitation from combined data in different subregions of the QTP; the black lines represent the mean values of the QTP, and the colored lines denote multiple subregions.

Atmosphere 2024 ,
15, 434 22 of 30GCMs might introduce more uncertainty into precipitation forecasts than temperature forecasts.

Figure 13 .
Figure 13.Signal-to-noise ratio (SNR) of MMA precipitation in different SSPs.The mean values of SNR across the entire QTP are denoted in each subgraph as m.The numbering (I-X) of each subregion (consistent with Figure 1) is also indicated in each subfigure.

Figure 13 .
Figure 13.Signal-to-noise ratio (SNR) of MMA precipitation in different SSPs.The mean values of SNR across the entire QTP are denoted in each subgraph as m.The numbering (I-X) of each subregion (consistent with Figure 1) is also indicated in each subfigure.

Figure 14 .
Figure 14.SNR, as in Figure 13 but for maximum temperature.Figure 14.SNR, as in Figure 13 but for maximum temperature.

Figure 14 .
Figure 14.SNR, as in Figure 13 but for maximum temperature.Figure 14.SNR, as in Figure 13 but for maximum temperature.

Figure 15 .
Figure 15.SNR, as in Figure 13 but for minimum temperature.

Figure 15 .
Figure 15.SNR, as in Figure 13 but for minimum temperature.

Table 1 .
Definition of variables in the HRP-QTP.

Table 2 .
Basic information on verification stations.

Table 3 .
Basic information about CMFD and GPQM.

Table 4 .
7GCMs used in this study.

Table 5 .
WD for different methods at 8 independent validation stations.The bold terms indicate the best method for statistical distribution validation.

Table 6 .
[30] for precipitation projection results from Cui et al.[30].The bold terms indicate SNRs greater than 1; SNRs were calculated using Equation (11) ignoring models with no data, and the division of subregions is slightly different from our study.