Oceanic climate changes threaten the sustainability of Asia’s water tower – Nature

Environment input data

All data used in this study are accessible from public databases. For the simulation of water vapour transport (see next section) by the FLEXPART v10.4 model^39,40, ERA-Interim atmospheric files covering 61 vertical levels from 0.1 to 1,000 hPa were downloaded and pre-processed by the flex_extract software (https://www.flexpart.eu/downloads). The spatial and temporal resolutions were 0.75° × 0.75° and 6 h, respectively. Monthly precipitation reanalysis data for 2003–2017 were sourced separately from ECMWF Reanalysis 5th Generation (ERA5)⁴¹ and from the Global Precipitation Climatology Centre (GPCC)⁴². These have a spatial resolution of 0.25° × 0.25°. Reanalysis of monthly evaporation and air temperature for 2003–2017 were sourced from ERA5. Monthly multimodel Coupled Model Intercomparison Project Phase 6 (CMIP6)-based precipitation, evaporation and air temperature data for 2003–2099 were acquired from the Center for Environmental Analysis. The considered CMIP6 models included ACCESS-ESM1.5, BCC-CSM2-MR, CanESM5, GFDL-ESM4, IPSL-CM6A-LR, MIROC6, MRI-ESM2.0 and NorESM2-LM under the r1i1p1f1 format. They were coded as CMIP6 models 1 to 8. The TWS anomaly for 2003–2017 was obtained from the GRACE RL05 Mascon dataset⁴³ (0.5° × 0.5°), whereas the TWS anomaly during 2020–2021 was taken from the GRACE-FO RL06 Mascon dataset⁴⁴ (0.25° × 0.25°). The GRACE RL05 data were interpolated to 0.25° × 0.25° using bilinear interpolation⁴⁵. Monthly mean ocean current data for 1993–2019 were acquired from the Ocean Surface Current Analysis Real-time (OSCAR) dataset⁴⁶ and have a spatial resolution of 0.33° × 0.33°. The snow-cover data (0.05° × 0.05°) were acquired from the National Snow and Ice Data Center⁴⁷. The continental map data used in the study have been generated on the basis of the continental region groups in the country-scale world map data³⁴. There is no conflict on the country boundaries because the map shown in the study is continental. The map data for the TP³⁵ were acquired from the National Tibetan Plateau Data Center, China.

Backward simulation of water vapour source regions for the TP

The Lagrangian particle dispersion model FLEXPART v10.4 simulates forward or backward trajectories of air masses (including their water vapour content) by representing a given air mass as a series of small-scale particles³⁹. At the particle scale, the balance between evaporation (e) and precipitation (p) determines the variation in specific humidity (equation (1)). When e − p > 0, a particle recharges its water content by evaporation from the surface; when e − p < 0, a particle loses its water through precipitation. For a given area, variation in water vapour (E − P; uppercase letters denote aggregated quantities) can then be determined as the sum of e − p for all particles in the area (equation (2)),

in which \(\frac{{\rm{d}}q}{{\rm{d}}t}\) denotes the temporal variation in the specific humidity of a particle, m is the mass of the particle and n is the number of particles over the considered area.

We used FLEXPART v10.4 and ERA-Interim climate reanalysis data to perform backward simulations of water vapour trajectories to the TP. This allows us to determine water vapour source regions. The model was then used to perform monthly simulations of trajectories for up to 500,000 particles. The number of the particles depends on the scale of the target region and computational resources^39,48. Individual simulations were performed at 6-h time steps for the period 2003–2017, enabling the determination of water vapour trajectory from source regions to the TP. Given that the residence time of the water vapour in the atmosphere is approximately 10 days (refs. ^49,50), we performed backward simulations from day 1 to day 10 of each month. A k-means method⁵¹ was then used to categorize all particle trajectories into 19 out of 20 main water vapour trajectories based on the spatial locations of particles in individual simulations. For the determination of the number of cluster trajectories, see Supplementary Note 1.

Because FLEXPART v10.4 is driven by ERA-Interim and does not adapt to the more recent ERA5 (this choice is due to the fact that ERA5 data supporting the run of FLEXPART v10.4 is not yet accessible to the public user; for details, see https://www.flexpart.eu/flex_extract/Ecmwf/access.html), a comparison was performed between the simulated total water vapour release and ERA5 and GPCC precipitation over the TP to investigate the sensitivity of simulation results⁵². The total water vapour release for particles residing over the TP was calculated by equation (2) when E − P < 0 (ref. ⁵³) and compared with the standardized regional average precipitation obtained from ERA5 and the GPCC. In both cases, comparisons were performed using correlation analysis. It should be noted that the total water vapour release is not equal to precipitation theoretically but it has the positive effect on the formation of the precipitation^39,48,52. The total released water vapour over the TP is highly correlated to both the GPCC and ERA5 precipitation, with correlation coefficients between 0.85 and 0.86 (P-value < 0.01; see Supplementary Fig. 11). We consider this to be an acceptable accuracy for the FLEXPART results (Supplementary Fig. 11).

On the basis of the initial simulation, 13 primary water vapour source regions were identified: Asia, Indian Ocean, North Africa, North Atlantic, Mediterranean Sea, Europe, Pacific Ocean, Red Sea, Tibet, North America, Caspian Sea, Black Sea and Arctic Ocean. However, because particles can lose or gain water vapour along their track, the water vapour released in the TP can reflect the contribution from both the source regions and the on-route transitions. Thus we calculated the relative contributions in two ways: (1) accounting for or (2) disregarding the on-route changes in water vapour. When accounting for water vapour variations (that is, recharge or loss) along the transitional routes, the considered air particles show positive variations in water vapour inside the source region (equation (4)), because the negative variations inside the source region indicate that the air particles release water vapour instead of acquiring and transporting water vapour from source regions to the TP. Meanwhile, the considered air particles should also be positive variations in water vapour before reaching the TP when accounting for on-route changes (equation (5)).

in which (E − P)_route and (E − P)_source refer to the water vapour variations along the transition routes and over the source regions, respectively. Similarly, r and s refer to the number of particles along the transition routes (r) and over the source regions (s). Because the air particles might gain or lose water vapour, the absolute variations in water vapour along the transition route should be calculated according to equation (5). For particles originating within the TP, (E − P)_source should be less than zero, thus recharging water vapour in the TP.

After the selection of particles reaching the TP (equations (3)–(5)), relative water vapour contributions from individual source regions were calculated as follows (equations (6) and (7)).

in which (E − P)_i refers to the total water vapour released over the TP from particles originating from the ith source region, (E − P)_Tibet refers to the total water vapour released in the TP and RC_i is the relative water vapour contribution obtained by the TP from the ith source region. When accounting for the on-route variation in water vapour, k_i is the number of selected particles (see equations (3)–(5)) from the ith source region. When disregarding water vapour on-route variations, k_i refers to the number of all particles from the ith source region.

However, there is another method⁵² that determines the relative contribution by dividing the total water vapour release over the target region by the moisture gains over the source regions. To avoid impacts induced by the differences between methods on the results, we also calculated the relative contribution using this method. It suggested that both methods come to the same conclusion that Asia, the Indian Ocean, North Africa and the North Atlantic are the four main water vapour source regions for the TP (Supplementary Fig. 12 and Fig. 1c). However, given that the moisture gains over the source regions are far larger than the total water vapour release in the target region, we applied the first method in this study.

Standardized trend item of the index

To avoid the impact of periodicity on correlations between variables, a time series of the variable during January 2003 to June 2017 was first smoothed using a moving average as implemented in the R function decompose. The derived trend series of the variable covers the period July 2003 to December 2016. Then, to ensure that different indices can be compared, the trend series of every variable was standardized over the period 2003–2016 by equation (8).

in which x denotes the trend series of the variable, X denotes the standardized version of index x, mean(x) and sd(x) are the arithmetic average and standard deviation, respectively, of x over the period 2003–2016 and i refers to the ith month during 2003–2016. Here the variables include PME, TWS, T or Snowcover, in which Snowcover and T denote the snow-covered area of the southern TP mountains and the air temperature, respectively. The relative variation for a standardized index was calculated by dividing the index variation by the maximum amplitude of all standardized indices (see equation (9)). Furthermore, we also provide mass/area change accordingly for every standardized index.

in which X_m and X_n denote the values at moments m and n in a time series X, RV_m,n is the relative variation between X_m and X_n, and maxV denotes the maximum amplitude of all standardized indices in this study. Because the maximum scale for all standardized indices is in the range −3 to 3, we applied maxV as 6 here.

Furthermore, to determine the spatial patterns of the trends, we used the modified Mann–Kendall trend analysis as implemented in the R package modifiedmk. This analysis was conducted at the scale of individual grid cells and covers the period 2003–2016.

Correlation analysis and maximum covariance analysis

Pearson’s correlation analysis was used to quantify the relation between standardized indices during 2003–2016. Cross-correlation analysis⁵⁴ was applied to evaluate the relation between (1) PME over the North Atlantic and (2) TWS in the 14 small-scale subregions across mid-latitude Eurasia with consideration of the optimum lags. The optimum lag is defined as the lag when the maximum correlations between two indices are detected. Both the Pearson and the cross-correlation analyses were performed using the R package stats.

Maximum covariance analysis was performed to further verify the covariation between spatial matrixes for any pair of indices. The spatial and temporal coefficients were calculated as follows:

in which L and R refer to the spatial matrixes of (1) PME over both the North Atlantic and the Indian Ocean and (2) TWS across Eurasia, respectively, and indices m and q refer to a given cell in the considered field, whereas the index k refers to the time step. C_l,r is the covariance matrix between L and R. U and V are the spatial modes, respectively, for L and R. PC_l,m and PC_r,q are the temporal coefficients for cells m in L and q in R. The superscript T in equations (10)–(14) denotes the transpose and Σ in equation (12) denotes the diagonal matrix. When the maximum correlation analysis is applied to evaluate the relation between T over ocean and land, L and R also refer to the spatial matrixes of (1) T over the related oceans and (2) T across Eurasia, respectively.

Attribution analysis of TWS variation in the central TP

Previous studies reported that both snowmelt in the southern TP mountains and increased precipitation over the central TP contributed to increased TWS in the central TP^5,31. However, the balance between precipitation and evaporation was not included in those studies. Here we applied two linear regression models to quantify the contributions of both snowmelt in the southern TP mountains and PME changes in the central TP to TWS changes in the central TP. The models are of the form:

in which TWS_r refers to the TWS in a given region r, which is either the southwest TP mountains (TPM1) or the southeast TP mountains (TPM2), TWS_TPS and PME_TPS denote the TWS and PME in the central TP (TPS), respectively, and a_r and b_r are the regression coefficients determined by linear regression.

Projections of the TWS non-deficit and deficit areas in the TP

The change in area not affected by TWS deficits (hereafter ‘TWS non-deficit area’) in the TP was determined by identifying grid cells with annual TWS > 0. On the basis of approximately synchronous variations in the TWS non-deficit area and (1) TWS (Fig. 3b) and (2) snow-covered area (Extended Data Fig. 6) in the southwest and southeast TP mountains, a further linear regression model was developed:

in which A is the TWS non-deficit area in the TP, Snowcover_r is the snow-covered area over region r, a_r and b_r are the regression coefficients and ε is the model residual.

As there are only very few CMIP6 models that simulate TWS_r and Snowcover_r explicitly, three random forest models (equation (17), m = 1, 2, 3) were developed for TWS_r and Snowcover_r projections. The models were first trained and validated on the basis of the GRACE-based TWS or the ERA5-based snow-cover area over region r during 2003–2016 and have the form:

in which Y_k denotes the ERA5-based PME (index k = 1) over the southeast North Atlantic, ERA5-based T (index k = 2) over the southwest TP mountains or ERA5-based T (index k = 3) over the southeast TP mountains during 2003–2016. RF denotes the random forest model. The first random forest model was formed by the dependent variable TWS_r (D_r,1) and independent variable Y₁ when m = 1. The second random forest model was formed by the dependent variable TWS_r (D_r,2) and independent variables Y₁ and Y₂ in southwest TP mountains or Y₁ and Y₃ in the southeast TP mountains when m = 2. The third random forest model was formed by the dependent variable Snowcover_r (D_r,3) and independent variable Y₂ in the southwest TP mountains or Y₃ in the southeast TP mountains when m = 3. The variable TWS_r was projected with (m = 1) and without (m = 2) considering T, because TWS_r is partially attributed to snow-cover reduction (Supplementary Fig. 7) resulting from local warming (Supplementary Fig. 8).

Equations (16) and (17) were both trained and validated using historical data, reflecting the actual relations between variables (Supplementary Figs. 13–15). These relations were based on the global atmospheric circulation systems (Figs. 1–4) including the westerlies and the Indian monsoon that are generated by the rotation of the Earth and its energy-balance system⁵⁵. Here we project future variations in the TWS non-deficit area in the TP, TWS and snow-covered area in the southern TP mountains using equations (16) and (17) with the CMIP6-modelled data as input.

The distributions of the CMIP6 climatic variables agree well with distributions of the ERA5 climatic variables when P < 0.05 (Supplementary Fig. 16). However, given the poor relationships between CMIP6 climate projections from individual model–scenario combinations and ERA5 climatic variables (Supplementary Figs. 17 and 18), (1) average (equation (18)) and (2) weighting sum (equations (19) and (20)) of climatic indices of CMIP6 models were first determined. These can strongly reduce the uncertainties⁵⁶ of CMIP6 projections (Supplementary Figs. 19 and 20) and the average and weighting sum of the climatic indices are obtained as follows:

in which s denotes the scenario (either SSP245 or SSP585), \({\bar{Y}}_{k,s}\) denotes the average index for all Y_k,i,s, W_k,s is the weighting sum of Y_k,i,s and Z_k,i,s denotes either the group of T (index k = 2) over the southwest (NATO1) and southeast (NATO3) North Atlantic or the group of T (index k = 3) over the four regions in the Indian Ocean (IO1–IO4) derived from the CMIP6 model i under scenario s because air temperature variations over oceans are related to the air temperature over the southern TP mountains (Extended Data Fig. 5). Equation (19) denotes the linear regression model with Y_k (see equation (17)) as the dependent variable and Y_k,i,s and Z_k,i,s as the independent variables. The weight w_k,i,s for Y_k,i,s, the weight c_k,i,s for Z_k,i,s and model residuals ε_k,s in equation (20) are determined by equation (19). For the selection of the index i in equation (19), see Supplementary Note 2. Furthermore, the (1) average and (2) weighting sum of TWS_r or Snowcover_r during 2020–2098 under scenario s are projected by equation (17) with \({\bar{Y}}_{k,s}\) (equation (18)) and W_k,s (equation (20)) as input, and are defined as (1) \({\bar{D}}_{r,m,s}\) and (2) WD_r,m,s, respectively.

The annual A during 2020–2098 under scenario s is projected by equation (16) with four combinations of \({\bar{D}}_{r,m,s}\) and WD_r,m,s (Supplementary Fig. 21). The average A (see equation (16)) of these projections (A_mean,s) is determined by equation (21), which has the form:

in which A_n,s denotes A_1,s projected by equation (16) with \({\bar{D}}_{r,1,s}\) and \({\bar{D}}_{r,3,s}\) as inputs, A_2,s projected by equation (16) with \({\bar{D}}_{r,2,s}\) and \({\bar{D}}_{r,3,s}\) as inputs, A_3,s projected by equation (16) with WD_r,1,s and WD_r,3,s as inputs or A_4,s projected by equation (16) with WD_r,2,s and WD_r,3,s as inputs. NegA_mean,s is the average area affected by TWS deficit (hereafter ‘TWS deficit area’) during 2020–2098 under scenarios s, computed as the difference between the A_mean,s and the total area of the TP (A_TP).

Evolution of the area under TWS deficit in the TP

Previous studies suggested that the TWS in the central TP is unpredictable owing to the strong interannual variability^5,31. By contrast, the TWS in the southwest TP and southeast TP are predictable based on the relation proposed in equation (17). Grid-scale TWS in the southern TP mountains is projected by equation (17) when m = 1. The centre points of grid cells, at which the annual sum of the monthly average TWS during 2020–2098 is negative (Fig. 4), on the north margin of the southern TP mountains forms the initial northern border of the TWS deficit area in the TP. Because we only study the spatial variation in the TWS deficit area inside the TP, the southern border of the TP is fixed and the area surrounded by the initial northern border and the southern border of the TP is defined as iniA (Extended Data Fig. 9). The variation in the TWS deficit area (ΔA) was determined by calculating the difference between NegA_mean,y,s and iniA:

in which NegA_mean,y,s is the TWS deficit area in the TP in year y under scenario s. When iniA > NegA_mean,y,s, the area with TWS < 0 shrinks southward by area |ΔA_y,s|. By contrast, when iniA < NegA_mean,y,s, the area with TWS < 0 expands northward by area |ΔA_y,s|.

However, ΔA alone is not enough to determine the areal evolution of the northern border of the TWS < 0 area. Given the stable radial expansion pattern of the TWS < 0 area with two negative centres under the minimum annual TWS sums in the southwest and southeast TP mountains since 2009 (Fig. 3 and Extended Data Fig. 6), future negative centres in these regions are first identified as the cells with the minimum annual sum of the monthly average TWS during 2020–2098 (see Fig. 4c,e) under SSP245 and SSP585 scenarios.

With the aforementioned centres as starting points, we determine the evolving paths as lines between the points forming the initial northern border and the centres in the southern TP mountains (Extended Data Fig. 9). Along the evolving paths, the evolving directions of points forming the initial northern border are thus determined as either northward (for iniA < NegA_mean,y,s) or southward (for iniA > NegA_mean,y,s) inside the TP. Furthermore, by setting the expansion steps Δd as 0°, 0.25°, 0.5°, 0.75°, 1°, 2°, 3°, 4°, 5°, 6° and 7°, a fitted line between the pairs of Δd and ΔA is identified in the scatter plot (Extended Data Fig. 9). For a given ΔA_y,s, determined by equation (23), a specific Δd could be derived by looking up the scatter plot (Extended Data Fig. 9) using the interpolation method. Once Δd was determined, by moving all of the points forming the initial northern border along the evolving directions, the northern border of the TWS deficit area inside the TP for a specific year and a specific scenario could be determined (Extended Data Fig. 9). We verify this area expansion model by the observed TWS deficit areas during 2009–2016 and find overlapping ratios of 75–91% (Supplementary Fig. 22). Performing the same test for 2020–2021 results in overlapping ratios of 75–86% (Supplementary Fig. 23). We considered these results to be sufficiently accurate for our analysis.