The Heihe River is an inland river that originates at the northeastern margin of the TP and disappears in the deserts of the downstream region (Figures 1A, B). It is formed by the confluence of the west tributary (YNG) and the east tributary (BBH) in the upper stream region (Gao et al., 2015). The part of the main river located on the TP and monitored by the Yingluoxia (YLX) hydrological station (Figure 1C) is defined as the upstream part (i.e., the UHRB) and has an annual mean precipitation of 200–700 mm and potential evapotranspiration of approximately 700 mm (Ma and Frank, 2006). The UHRB (37.72° N—39.09° N, 98.57° E—101.16° E) has a drainage area of 10,005 km², with elevations ranging from 1,681 m to 5,015 m.

The middle stream of the Heihe River, which encompasses the piedmont lower alluvial fan and fluvial plain, is the primary irrigated region (Song et al., 2020). This region accounts for 97.3% of the total population and generates 96.5% of the GDP of the entire Heihe River Basin (Cheng et al., 2015). However, precipitation in this part is relatively low (100–200 mm) and the potential evapotranspiration is higher (2000–2,500 mm) than that in the upstream part (Cheng et al., 2015). Moreover, the UHRB provides almost 70% of the total river runoff to support the middle and downstream regions (Gao et al., 2015), making it the most crucial water source region for the entire basin. Wu et al. (2015) simulated runoff components in the UHRB and found that glacier melting contributed 8.9% to the streamflow, highlighting the importance of precipitation runoff (including rainfall and snowfall runoff) for the entire basin.

Daily meteorological data covering 62 years, including precipitation amount, maximum air temperature, minimum air temperature, and average air temperature, were collected from six meteorological stations located in and around the UHRB (Figure 1C). These data were obtained from the Resource and Environment Science and Data Center (Table 1). Additionally, to validate the methods to discriminate precipitation type, precipitation-type data covering 24 discontinuous years were sourced from unpublished hydrologic yearbooks of China (Table 1). The observation and recording of the precipitation amount and types by hydrological stations followed the standard procedure published by the Ministry of Water Resources of the People’s Republic of China (MOWR, 1975). The recorded precipitation types included rain, snow, sleet, hail mixed with snow, and hail mixed with rain.

To differentiate between types of precipitation, this study applied three distinct methods, namely, the one-threshold air temperature method (Chen et al., 2014), the two-threshold air temperature method (Zhang et al., 2015), and Ding’s comprehensive model (Ding et al., 2014). The one-threshold air temperature method uses one air temperature to distinguish between rainfall and snowfall. The two-threshold air temperature method uses two different air temperatures in different situations, typically in different seasons. Ding’s comprehensive model considers air temperature, relative wind speed, relative humidity, air pressure, weather station elevation, and wet-bulb temperature. Ding’s method is detailed elsewhere (Ding et al., 2014).

The probability of detection (POD) (Chen et al., 2014; Pan et al., 2015) was calculated to verify the accuracy of the aforementioned precipitation-type estimation methods. The POD was calculated using the following formula: P O D = C A + C , (1) where A represents the number of snowfalls correctly estimated and C denotes the number of snowfalls missed. Thus, A + C is the recorded number of snowfalls.

To investigate the spatial variability in precipitation and air temperature within the UHRB, the contribution of each weather station in the target area was determined using the Thiessen polygon method in ArcGIS 10.5 and the locations of the six weather stations (Supplementary Figure S1). The Thiessen polygon method was proposed for spatial interpolation (Thiessen, 1911), which is based on proximal mapping; i.e., the nearest-distance neighbor. Compared to other interpolation methods, such as inverse distance weighting (IDW), Kriging, and spline interpolation, the Thiessen polygon method is more convenient for interpolating masses of daily data. This method has been successfully applied in many studies applying precipitation interpolations (Li et al., 2018; Malik and Kumar, 2020; Yanos, 2022).

The weighted average precipitation for the whole basin was computed by multiplying the precipitation at each station by its assigned percentage of area (Yanos, 2022). The formula is as follows: P a v e = P a C a + P B C B + · · · + P N C N , (5) where P _ave is the weighted average rainfall for the whole basin; P _A , P _B , … , P _N are the observed precipitation values for stations A, B, … , N, respectively; and C _A , C _B , … , C _N are weighted factors calculated using the Thiessen polygon method for stations A, B, … , N, respectively.

Note that the simple spatial interpolation based on the station data was applied in this study, even though it would potentially not fully capture precipitation gradients in mountainous areas (Zhang et al., 2015).

The Kendall rank correlation test (Kendall, 1934; Mann, 1945) was employed to detect trends in the time series to identify general change trends for all indices. For indices that exhibited significant change trends, the Pettitt non-parametric approach was applied to identify the abrupt change points (Pettitt, 1979).

The classification of precipitation into rainfall and snowfall is important as these types of precipitation have significantly different effects on runoff processes, soil erosion, and disaster response. However, it is often difficult to obtain records on precipitation types from public historical records.

Ding’s method was developed to distinguish precipitation types including rain, snow, and sleet based on wet-bulb temperature (Ding et al., 2014). Deng et al. (2017) applied this method to distinguish snowfall in the Tibetan Plateau and analyze snowfall changes. Han et al. (2019) applied the CREST-snow hydrologic model coupled with Ding’s method to simulate the flow regime in the Lancang River Basin. Su et al. (2022) employed Ding’s method to calculate rainfall and snowfall to evaluate the glacier changes.

The snowfall threshold-temperature method is typically used in most hydrological models, and the parameters of this method can be adjusted within a reasonable range. Han et al. (2010) analyzed precipitation-type records from 643 weather stations in China from 1961 to 1979 and found that the one-threshold temperature method showed the most appropriate performance. Additionally, a spatial distribution pattern has been identified in China, wherein the threshold temperatures are higher in plateau areas (Han et al., 2010; Chen et al., 2014). Chen et al. (2014) indicated that the one-threshold temperature for daily rainfall and snowfall in the UHRB area was in the range of 3.5°C–5.5°C. Han et al. (2010) also suggested that the two-threshold temperature method may be suitable in some arid areas. Zhang et al. (2015) found that the snowfall temperature was negatively correlated with relative humidity in the northern region of the Himalayas, indicating that the threshold temperature might be higher in dry seasons than in wet seasons.

In this study, we first determined the best thresholds for the one-threshold and the two-threshold temperature methods by traversing the threshold temperature between 3.0°C and 5.5°C (Supplementary Figure S2). Except for instances in which there were no differences in extremely cold months (i.e., January, February, November, and December), the best threshold temperature for March and April with the most snowfalls was 4.5°C (i.e., the one-threshold temperature method T_4.5); meanwhile, for May to October, the best snowfall threshold temperature was 3.6°C. Thus, the best solutions for the two-threshold-temperature method were 3.6°C for May to October and 4.5°C for November to April of the next year (T_{3.6_4.5}). We also classified precipitation type based on Ding’s method using data covering 24 years. As most hydrological models require only rainfall and snowfall as inputs, this study focused on classifying these two precipitation types. Moreover, since the number of sleet events in the study area was relatively small, the sleet events were considered rainfall as in other studies (Chen et al., 2014; Deng et al., 2017; Su et al., 2022). The total number of snow days per month and over the 24-year period for the three methods (T_{3.6_4.5}, T_4.5, and Ding’s method) and their probability of detection (POD) values are shown in Figures 2A, B.

Overall, the T_{3.6_4.5} method performed better than the T_4.5 method and Ding’s method, with a POD value as low as 1.9%. Considerable over- and underestimations were noted when employing the T_4.5 (+10.3%) and Ding’s method (−6.9%) for the 24 years of data. Snowfall events constituted the absolute majority during January, February, March, April, November, and December, whereas rainfall events were dominant in June, July, and August. Thus, all three methods performed evenly in these 9 months. However, the accuracies of the three methods differed during the transitional months of May, September, and October. Although the POD values for the T_{3.6_4.5} method for these months were 5.0%, −20.0%, and 4.2%, respectively, those for the T_4.5 method (50.0%, 120.0%, and 12.5%) and Ding’s method (−40.0%, −80.0%, and −10.4%) were not acceptable. In a precipitation-type estimation study based on data from 643 stations across China, Chen et al. (2014) also found that 51% of the stations had POD values <2% when applying the air temperature threshold method. Based on these findings, the T_{3.6_4.5} method was considered the most suitable for precipitation-type differentiation in the UHRB.

Using the T_{3.6_4.5} method for discrimination of precipitation type, the daily rainfall and snowfall events during 1960–2021 were distinguished at the six weather stations located in and around the UHRB. Figure 3 reveals distinct spatial patterns in the average annual precipitation, rainfall, and snowfall in the UHRB. The UHRB has an average annual precipitation of 381.6 mm, with 81.9% rainfall (312.6 mm) and 18.1% snowfall (69.0 mm). Both precipitation and rainfall amounts exhibited a decrease from east to west across the UHRB, which is consistent with the spatial distribution patterns of the WRF model outputs based on station data (Pan et al., 2015) and some precipitation products (e.g., DCMAP, ITP-F, and DCMAP-MicroMet) (Pan et al., 2014). However, the snowfall amount was relatively higher in the west tributary. Owing to the influences of the East Asian monsoon in summer and the Westerlies in winter (Yao et al., 2022), which carry different amounts of moisture, differences were observed in the distributions of total precipitation, rainfall, and snowfall distribution from east to west in the UHRB (Cheng et al., 2015). Consequently, while the climates in the west and east tributary regions of the UHRB are similar, the precipitation amounts and types in these regions exhibit differences due to variations in the locations and underlying surfaces.

As the UHRB is influenced by the East Asian monsoon during summer and the Westerlies during winter, notable seasonal variations were observed in total precipitation, rainfall, snowfall, and temperature (Figure 4). Snowfall accounted for an average of 18.1% of the precipitation amount, with significant variations observed across different months. Snowfall constituted 100% of the precipitation events from November to March. As the Westerlies and the East Asian monsoon shift, both temperature and precipitation rose from May, reaching their peaks in July (94.2 mm of precipitation and a mean temperature of 12.4°C), and then decreasing until the next cycle. Moreover, as the average and minimum air temperatures exceed 0°C during May–September, rainfall events dominate during these months, accounting for 60.0%, 93.8%, 99.4%, 99.0%, and 84.9% of the precipitation amounts in May, June, July, August, and September, respectively. The annual snowfall occurred during April and May (13.6 mm and 16.6 mm). Owing to the accumulation of snowfall in the basin from October to March (25.4 mm in total), spring floods typically occur in April and May (Zhu et al., 2020). Summer floods typically occur in July and August (Luo et al., 2017) because the highest rainfall is received during the East Asian summer monsoon.

The annual total precipitation, rainfall, snowfall, mean temperature (Tmean), maximum temperature (Tmax), and minimum temperature (Tmin) were calculated from daily measurements of these parameters. The change trends of these annual values were assessed using Kendall’s rank correlation tests with the year, and the change rates were estimated using linear regression. The spatial distributions of these parameters are shown in Figure 5. Over the past 62 years, except for station 52,765, which is located in the southeasternmost station and exhibited a non-significant increasing trend, all other stations exhibited significant increases in precipitation and rainfall, with rates ranging from 7.35 to 14.68 mm/10y and 6.42–15.46 mm/10y, respectively. The spatially averaged precipitation and rainfall amounts also showed significant increasing trends, with rates of 13.45 mm/10y (p < 0.001) and 13.02 mm/10y (p < 0.001), respectively. However, we observed no significant change in annual snowfall across all stations, with a slightly increasing trend of 0.43 mm/10y (p = 0.697) for the spatial average. It should be noted the area influenced by both the East Asian monsoon and the Westerlies may experience more severe changes in total precipitation, rainfall, and snowfall. Additionally, in the context of global climate warming, the average, maximum, and minimum temperatures showed notable increases, with spatially averaged increase rates of 0.32°C/10y, 0.26°C/10y, and 0.39°C/10y, respectively, which were three times (Qiu, 2008) the global average warming rate (0.12°C/10y during 1951–2012) (IPCC, 2014) and slightly lower than the average for the TP (0.35°C/10y during 1970–2014) (Yao et al., 2019).

The changes in total precipitation and rainfall in the UHRB during 1960–2021 displayed an abrupt and significant shift in 2002 (p < 0.05) according to the Pettitt (1979) test (Figure 6, the year with the highest Sk value). This abrupt year was in good agreement with the average situation across the entire TP (Wang et al., 2016). No abrupt year was found for snowfall as it did not exhibit any significant change trend (p > 0.05). In addition, the average, maximum, and minimum air temperatures all exhibited stepwise increases around 1996 and 1997, which was consistent with the average for the entire TP (Zhang et al., 2022). Furthermore, the periods of 1960–2001 and 2002–2021 were selected as consolidated periods for comparison before and after the abrupt year (Table 3). The difference in total precipitation in the UHRB was greater, and the decrease in the snow/rainfall ratio (S/R) was negligible compared to those in the entire TP (Wang et al., 2016). Furthermore, snowfall showed a slight increase in the UHRB, whereas it exhibited a significant decrease in the TP (Wang et al., 2016). However, the changes in temperature were comparable with those in the entire TP (Zhang et al., 2022).

The total precipitation, rainfall, snowfall, and temperatures exhibited varying monthly change trends (Figure 7). The Tmax for all months, except for April and May, as well as the Tmean and Tmin, demonstrated significant increases (p < 0.01). While all months exhibited an increase in precipitation, the increases were only significant in January, June, and December, owing to the significant increases in rainfall in June and December, and snowfall in January and December. In contrast, significant decreases in snowfall were observed in June and September, possibly due to the warmer climate and the shift of snowfall to rainfall.

Overall, the changes in precipitation and temperature suggested that despite similar climate warming backgrounds, the changes in precipitation can be more complex due to the influence of general atmospheric circulation, topography, and other factors.

Consecutive dry days (CDD), annual maximum consecutive 5-day precipitation (RX5DAY), annual count of days with precipitation >10 mm (P10 mm), annual total precipitation above the 99th percentile (P99P), annual total precipitation above the 95th percentile (P95P), and simple daily intensity index (SDII) of the precipitation indices computed by ClimDex software were selected for analysis in the present study. However, ClimDex software does not distinguish the precipitation types. Considering that the P99P values in some years were zero, we manually calculated the annual total rainfall above the 95th percentile (R95P) and the annual total snowfall above the 95th percentile (S95P) for analysis of extreme precipitation in the UHRB. The threshold values of P99P, P95P, R95P, and S95P at the six stations are listed in Table 4. Some studies analyzed extreme or erosive precipitation without distinguishing between rainfall and snowfall, which may lead to uncertainties. Therefore, based on the rainfall distinguished, rainfall erosivity (REro) values were also calculated as one of the indices affected by extreme rainfall. The ranges and changes of the indices for each station and spatially averaged values are shown in Figure 8 and Table 5.

Six of the temperature-related indices were also generated by the ClimDex model. The annual change trends of the six stations and the UHRB average are shown in Figure 10. Unlike the precipitation-related indices, the temperature-related indices changed significantly in almost all six stations. The cold day-related indices, namely, the annual count of days with Tmax <0°C (ID0), the annual count of days with Tmin <0°C” (FD0), and the percentage of days when Tmin was below the 10th percentile (TN10P), showed compelling decreasing trends at all sites, indicating significant changes in minimum air temperatures. In contrast, the high-temperature-related indices, namely, the percentage of days with Tmax above the 90th percentile (TX90P) and the annual count of days with at least six consecutive days with TX above the 90th percentile (WSDI), all indicated intensive increasing warmth in the UHRB. The growing season length (GSL) in the UHRB was much shorter than that in other areas outside the TP, ranging from 134.1 to 176.8 days/y, but increased noticeably at 3.6 days/10y, which was almost three times the average rate of change in the Yangtze River basin (1.3 days/10y) (Cheng et al., 2022).

The elevation, latitude, and longitude correlations with precipitation-related indices changes used for the spatial pattern and geographical factor analysis are listed in Supplementary Table S1. Due to the limited number of stations, most of the correlations were not significant. However, P, Rain, Snow, RX5DAY, P10mm, P99P, P95P, R95P, SDII, and REro changes were more pronounced in stations in the UHRB at higher altitudes, indicating that higher elevations were more strongly affected by climate change than lower elevations. Furthermore, all the precipitation-related indices were negatively correlated with latitude and were significantly correlated with total precipitation and rainfall values, which indicated more severe changes in the southern part of the UHRB.

El Niño-Southern Oscillation (ENSO), as an indicator of global climate change, has been shown to modulate the precipitation in the TP (Yao et al., 2022) and to remotely drive extreme precipitation across the TP (Bothe et al., 2010). Xiong et al. (2019) found that the ENSO had a significant positive correlation with P20 mm (annual count of days with precipitation >20 mm) and SDII during 1995–2001 in the TP. A significant positive correlation between the ENSO and P10 mm was also confirmed in the Qilian Mountains (2015). Qin et al. (2010) found that the intensification of the hydrological cycle in the 20th century could be attributable to regional to largescale temperature increases.

The correlation analysis results between Niño 4 SST index and climate indices of this study are listed in Supplementary Table S2. The air temperature-related indices, including Tmean, Tmax, Tmin, ID0, TX90P, WSDI, and GSL, were significantly correlated with the ENSO, indicating that temperature increases and extreme temperatures were largely influenced by the ENSO. Regarding the precipitation-related indices, P, Rain, RX5DAY, P10 mm, R99P, P95P, R95P, SDII, and REro were positively correlated with the ENSO, but the difference was only significant for P10 mm (p = 0.015), which is consistent with the results reported by Cheng et al. (2015).

In addition to global warming, large-scale atmospheric circulation, including the Indian monsoon, westerlies, and the East Asian monsoon, is a main driver of precipitation changes in the TP (Yao et al., 2022). Cheng et al. (2015) proved that rainfall extremes (P10 mm and total precipitation) in the UHRB had a strong negative correlation with the East Asian summer monsoon index and a significant positive correlation with the western Pacific subtropical high-intensity index.

This result indicated that the extreme rainfall events in the UHRB were influenced by both global climate warming and large-scale atmospheric circulation.