Contrasting runoff trends between dry and wet parts of eastern Tibetan Plateau

As the “Asian Water Tower”, the Tibetan Plateau (TP) provides water resources for more than 1.4 billion people, but suffers from climatic and environmental changes, followed by the changes in water balance components. We used state-of-the-art satellite-based products to estimate spatial and temporal variations and trends in annual precipitation, evapotranspiration and total water storage change across eastern TP, which were then used to reconstruct an annual runoff variability series for 2003–2014. The basin-scale reconstructed streamflow variability matched well with gauge observations for five large rivers. Annual runoff increased strongly in dry part because of increases in precipitation, but decreased in wet part because of decreases in precipitation, aggravated by noticeable increases in evapotranspiration in the north of wet part. Although precipitation primarily governed temporal-spatial pattern of runoff, total water storage change contributed greatly to runoff variation in regions with wide-spread permanent snow/ice or permafrost. Our study indicates that the contrasting runoff trends between the dry and wet parts of eastern TP requires a change in water security strategy, and attention should be paid to the negative water resources impacts detected for southwestern part which has undergone vast glacier retreat and decreasing precipitation.

. Land cover characteristics for each basin. Percentage (%) of each land cover type for each basin was calculated from MCD12 product in 2006.

Evapotranspiration products introduction
Evapotranspiration (ET) models are usually very complicated, requiring remote sensed data, meteorological data, and other ancillary information (such as land cover type, land surface roughness, DEM etc.) 2 . It is expected that uncertainties could be large from a single ET product because of errors in input data, imperfect modeling of physical process, error propagation and failure of achieving global optimum solution in parameterization. To reduce the uncertainties and bias, we obtained four state-of-the-art diagnostic ET products (PML 3 , GLEAM 4 , MOD16 5 , P-LSH 6  any systematic bias for a particular product which could be large due to an absence of the water balance constraint. Summary information of the four ET products is listed in Table S3. Among them, GLEAM is based on Priestly-Taylor (PT) approach, and the other three are based on Penman-Monteith (PM) approach. Their input data and the modeling process are largely different, so it can be assumed that these products are independent and an ensemble of these products can increase the accuracy of ET anomaly estimates. To show spatial pattern of the four ET products in our study region, we mapped mean annual ET for each product over 2004-2009 (i.e., the reference period for anomaly calculation). All the products generally show wet-to-dry gradient from southeast to northwest, but differ greatly in the absolute values and spatial details (Fig.   S1). For example, regional mean value of P-LSH product is almost twice that of PML product; GLEAM and MODIS have similar regional mean values but show very different spatial pattern. Nonetheless, these differences are expected, justifying the necessity of applying the median anomaly analysis for the four ET products.

∆ response to precipitation and snow cover change
In most regions, ∆ and P show strong correlations, especially in northern part, indicating that P drives ∆TWS. For other regions featured with wetter climate, the correlations are weak (Fig. S2). This is consistent with previous findings that ∆ of catchments in dry climate zones are more sensitive to P than catchments in humid climates 7 .  Year on year Changes in minimum snow coverage (km 2 ) TWS (mm year -1 ) information, we used annual minimum snow coverage as a proxy. The MODIS snow cover product data (MOD10) at 0.05 degree spatial resolution were processed to delineate annual minimum snow cover assuming snow that doesn't melt in summer is the permanent snow/ice 8 . We further explored the relationship between ∆ and the year-on-year changes in permanent snow coverage for Jiayuqiao basin of Salween River located in south-western region (Fig. S3). ∆ is strongly non-linear correlated to permanent snow/ice cover change, and becomes more sensitive when the area of permanent snow/ice coverage is reduced by more than 200 km 2 compared to the previous summer value. This reflects the non-linear relationship between snow/ice area change and mass change 9 .

Application of gain factors to ∆
Due to the destriping and smoothing effect caused by filters, gain factor usually need to be applied to the GRACE solutions over land. However, as far as mascon (mass concentration blocks) data is concerned, it is not a mandate to apply gain factors 10 since geophysical constrains are implemented in calculating mascon fields and the gain factors are very close to 1. Nonetheless, we still compared two sets of Q results obtained with and without gain factors. No substantial difference was observed for the five basins, confirming that Q estimation performance cannot be improved by applying gain factors (Fig. S4a). We further investigated the spatial pattern of gain factors with spatial resolution of 0.5 degrees (Fig. S4b), and found that for most of our study site they are close to 1 (i.e., the expected value), yet in the northwest, both extremely small (even negative) and large values exist. This is unreasonable and is likely caused by the inadequacy of the Community Land Model (CLM) in modelling hydrological process over northwest which has widespread permafrost and where two big lakes Lake Zhaling and Lake Eling are located with the areas 528 km 2 and 628 km 2 , respectively 11 . Although these abnormal gain factors don't affect Q results aggregated at basin scale, they could have negative impacts at grid level, so we chose to use raw ∆ data despite of its coarse resolution (3 degrees).

TRMM precipitation validation
We first evaluated TRMM P data against local observations at 87 meteorological stations located within the modelling domain. There was a high consistency in the observed and TRMM P trends (Fig. S5a). For most stations, P observations were strongly correlated TRMM data (r > 0.75), indicating a reasonable match in inter-annual variations (Fig. S5b). The TRMM data also captured dispersion features a b of probability distribution in P observations, which is measured by Coefficients of Variation (CV) (Fig. S5c). Regional averaged TRMM P was higher than the observed average (Fig. S5d). This is expected because TRMM can observe heavy precipitation over high altitudes while most stations are located at lower altitude areas with less precipitation 12,13,14 . Second, we checked the performance of TRMM data in constructing the inter-annual variation in Q. Other precipitation dataset, including west01 15 , MSWEP 16 , nearest interpolation results of the monthly gauge observed P, and thin-plate regression interpolation 17 results using TRMM as a covariate, were also considered for comparison. Results show that the inter-annual variation in Q at five gauges can be best explained with the raw TRMM P data (Table S4). Specifically, TRMM P provided the best validation statistics for Luning, Zhimenda and Jiayuqiao basins, and achieved acceptable accuracy for Tangnaihai and Changdu basins. In comparison, other P dataset produced problematic results for Zhimenda basin and/or Jiayuqiao basin, though the results for Tangnaihai and Changdu basins were slightly better. It is noted that in our study we only focused on inter-annual P data variations, which is insensitive to systematic errors. If the aim is to estimate the absolute values of Q, calibrations of P should be considered.