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

Water resources sustainability in High Mountain Asia (HMA) surrounding the Tibetan Plateau (TP)—known as Asia’s water tower—has triggered widespread concerns because HMA protects millions of people against water stress1,2. However, the mechanisms behind the heterogeneous trends observed in terrestrial water storage (TWS) over the TP remain poorly understood. Here we use a Lagrangian particle dispersion model and satellite observations to attribute about 1 Gt of monthly TWS decline in the southern TP during 2003–2016 to westerlies-carried deficit in precipitation minus evaporation (PME) from the southeast North Atlantic. We further show that HMA blocks the propagation of PME deficit into the central TP, causing a monthly TWS increase by about 0.5 Gt. Furthermore, warming-induced snow and glacial melt as well as drying-induced TWS depletion in HMA weaken the blocking of HMA’s mountains, causing persistent northward expansion of the TP’s TWS deficit since 2009. Future projections under two emissions scenarios verified by satellite observations during 2020–2021 indicate that, by the end of the twenty-first century, up to 84% (for scenario SSP245) and 97% (for scenario SSP585) of the TP could be afflicted by TWS deficits. Our findings indicate a trajectory towards unsustainable water systems in HMA that could exacerbate downstream water stress.

The authors present a very comprehensive study to possible explanations for the spatial variations in trends of terrestrial water storage on the Tibetan Plateau. They use a combination of the FLEXPART particle dispersion model with reanalysis climate data and GRACE remotely sensed TWS data, and argue that TWS decline in the southern TP can be attributed to a deficit in precipitation minus evapotranspiration from the southeast North Atlantic. They further argue that the high mountain ranges in Asia block the propagation of the deficit further onto the central TP which leads to increasing TWS there. Lastly they demonstrate using CMIP6 GCM data how this blocking could weaken in the future, causing also the central TP to witness future TWS deficits. To my knowledge this is an original approach, it is mostly well described and components of the analysis thoroughly validated. The findings are of broad interest and provide new knowledge about the possible mechanisms of variability in climate change trends across the TP. In some cases I am not fully convinced about the robustness of intermediate conclusions, which propagate into the further analysis. I believe these need further clarification. These are included in the specific comments below.
Title: The title seems somewhat strange in that it could be read as that the ocean is drying up. Consider rewording. L84-85: I do not understand the statement 'Since glacier melting dominated the TWS changes in the high mountains instead of elevation'. Figure S1: Suggest to change the y-axis for panels b-d to be able to see more details, like also done in other rows in the figure. L139-140: Looking at Figs S1 and S2, the contributions from North Africa are about equal to the North Atlantic. Therefore, I wonder how the authors exactly come to the conclusion that contributions from the North Atlantic dominate over contributions from North Africa? This conclusion is important as it propagates further in the analysis when focus is going to the North Atlantic region. I would guess that the ocean being a body of water supplies more water vapor than the northern African land mass with the dry Sahara, but I cannot deduct that from Figs S1-S2. In Figs S3-S9 correlations for the North Atlantic and Indian Ocean are shown and North Africa is not further investigated whereas the data in Figs S1 and S2 suggest it should. I am no expert in particle dispersion modeling, so the explanation could be something obvious, but then it requires a bit more explanation for the reader. Also, while Figs S3-S9 indeed show better correlations for North Atlantic compared to Indian Ocean, it seems that Indian Ocean still plays a large role, in particular IO4 (Fig S8, panel c).
L198-201: I find the preceding argumentation and interpretation of Extended Data Figures 3 and 4 for the point that the high mountains block the propagation of the PME deficit hard to follow. Therefore it also does not fully convince me and I am not sure how valid the point is. I suggest to revisit lines 186-201 and refer to particular panels in the referred figures and elaborate explanation for this important part of the paper.
L324-327: Why were these particular CMIP6 models included and others not? I also believe MICRO6 is a typo and should be MIROC6. : Comparing timeseries of CMIP6 output directly year-to-year to reanalysis data in my experience never provides good correlation. I believe it is more valid to compare climate characteristics averaged over longer periods, like distributions, degree of variability.
Supplementary Table 1: I suggest to omit this table, the statement in the Methods section is sufficient.

Summary comments and reply
Dear authors, I have reviewed your manuscript entitled ''Sustainability of Asian water tower threatened by drying and warming ocean''.
Reply: Thank you for your time in reviewing our manuscript. To make it easy to follow, under each comment, we have included appropriate reference to the revised manuscript (e.g., line numbers), wherever appropriate. Quoted text from the manuscript is provided in italics.
The results are of interest to the scientific community but also as information for taking decisions for the benefit of the economy and the population. The authors found interesting and unexpected results, such as the relative importance of the North Atlantic over the Indian Ocean as a moisture source for the TP. However, there are some details I would like to share with you. First, improve the introduction of the manuscript in order to make the objective of the study clearer. Indeed, it is not clear, at least for me, when the authors focus on the relationships between PME and TWS over 3 subregions (SR1, SR2 and SR3) and the contribution from oceanic sources. The title of the study is then, in my opinion, not consequent with these analyses, because the study should be focused on the role of the sources on the precipitation and TWS in the TP. I would appreciate your opinion about it, because to be published in this journal, the study must be concise and focused on the objective of the study without rambling on collateral topics.
Reply: Thank you for your positive comments on our work. We improved the interpretation parts of the results in the revised manuscript and further analysis was also added.
The reason why we analyzed "the relationships between PME and TWS over 3 subregions (SR1,SR2 and SR3) and the contribution from oceanic sources" is explained as follows.
Firstly, it is mainly because that we want to demonstrate the transition process of the changes in oceanic climate to TP. The TP is the part of the SR1-3 but physically separated from the oceans. And the oceans (e.g., Indian Ocean and North Atlantic) are demonstrated in this revision that recharge the water vapor of air particles over lands (e.g., Asia and Africa) that further transit water vapor into the TP (Supplementary Fig.   1). Thus, the changes in oceanic climates might transit to TP by atmospheric circulations leading to consistent TWS variations on route (e.g., SR1 and SR2). If the transition process truly exists, the consistent relationships between PME over oceans and TWS in SR1-3 (including TP) are expected. In this study, we verified the transition process by combining the backward trace analyses and relation analyses. It suggests that the TWS in SR1-3 and southern TP mountains are all significantly related to the PME over northwest (NATO2) and southeast (NATO3) North Atlantic (Supplementary Secondly, it is because that we wanted to verify whether the TWS in lands are influenced by the disturbance on the local PME induced by changes in oceanic climate. We analyzed the relationships between PME over SR1-3 and PME over oceans and the relations between PME and TWS over SR1-3. It suggests that the TWS in lands are not Overall, the analysis for the "relationships between PME and TWS over 3 subregions (SR1,SR2 and SR3) and the contribution from oceanic sources" are necessary for the study focusing on the impacts on TWS in TP induced by changes in oceanic climates.
In this line, I also expected to find the analysis of the optimum lag relationship between the PME and TWS, which was not assessed, or an explanation about it. And I mention this because TWS may be influenced in the medium and long term by the moisture supply from the external sources, taking into account that TWS includes surface but also underground freshwater.
Reply: Yes, we agree with you on this comment. And we had performed the cross- Focusing on regions within 20°N-50°N, considering optimum lags, we show high cross correlations between the PME deficit over the southeast North Atlantic and the TWS changes within 12 out of 14 sub-regions discerned along the water vapor propagation routes towards the TP Supplementary Fig. 6).
We also analyzed the cross-correlation between TWS in southern TP mountains and PME over North Atlantic with consideration of the optimum lags (see Lines 161-164).

Lines 161-164
And the cross-correlation analysis suggests the decreased TWS over the southern TP mountains can be attributed primarily to the PME deficit from the southeast North Atlantic with considering optimum lags .
The authors confirmed the crucial importance of Asia as moisture source for the TP.
But along the paper there is no explanation about the possible role of recycling in this region, particularly the recycling of precipitation attributed to the moisture transport from the Indian Ocean, and the relationship of the moisture export from Asia to the TP, and the moisture contribution from the IO to Asia. This chain could increase the importance given to the IO in the study, and must be discussed.
Reply: We have performed additional analysis and demonstrated the process of the moisture transport from oceans to TP. Given that the SR3 covers most area of the TP, the additional backward trace analyses were performed for particles initially residing over the SR1 and SR2 in this revision. According to Supplementary Fig. 1, it suggests that the oceans (e.g., Indian Ocean and North Atlantic) recharge the water vapor of air particles over lands (e.g., Asia and Africa) that further transit water vapor into the TP by the westerlies and the intermedia current (Figure 1a-b).
We agree with you on that the IO also plays a key role on impacting the climatic conditions in Asia and TP, and have attached discussions to the revised version of the manuscript (Lines 101-109).

Lines 101-109: Revised discussion
Distance to the TP has negative influence on the contributions from source regions to TP. Even though, it should be noted that the ocean is the primary continental water source region in a large-scale circulation. The oceans (e.g., Indian Ocean and North Atlantic) recharge the water vapor of air particles over lands (e.g., Asia and Africa) that further transit water vapor into the TP (Supplementary Fig. 1) via the two main trajectories. The first is from the North Atlantic to the western TP transported by the westerlies (Figure 1a-b), and the other is the intermedia current (see Figure 1a-b) that brings water vapor from North Atlantic and Indian Ocean into the southern TP.
In the title of the manuscript, you mention: Ocean warming. I did not find in the manuscript the rate of warming of the oceanic moisture sources (under the historical and future scenarios)?
Reply: We agree with you on that the title might cause some misunderstanding. We have now revised the title to "Changes in oceanic climate threaten the sustainability of   The approach is relatively novel and has been previously used for similar studies.
Considering the above comments, I suggest major revisions must be done to be acceptable for publication. I invite the authors to revise their manuscript to address the concerns (some listed below) before a final decision by the Editorial is reached. I strongly ask you to clarify some aspects of the methodology, see the comments below.
Reply: Thank you for your positive comments on our work. We have revised the manuscript by following your specific suggestions.
Reply: We replaced 'the' with 'a in the revised manuscript. In Figure 1a, is confusing the almost similar colors (yellow-reddish) for some of the arrows that you use to denote the backward direction of air parcels from the regions (e.g AS, EU, etc.) Reply: Thank you for your suggestions. We also apply the red-yellow colors for trajectories to distinguish them from background information in Figure  For your second puzzle that "decreasing as well the importance of longer regions as moisture sources for the target region", we also agree with you on this point. In this revision, we added the discussion about the impact on the relative contribution induced by the distance (Lines 101-109). And when the on-route changes are considered, the relative contributions from Asia, Indian Ocean, Africa and North Atlantic Ocean ranks from the largest to the smallest while their distances to TP ranks from the smallest to the largest (Figure 1c). This result complies with the aforementioned statement on the impacts of distance. However, when the on-route changes are disregarded, we find that the contribution from the North Atlantic is larger than the one from the Indian Ocean.
This is because that the number of particles from North Atlantic is far larger than the number of particles from Indian Ocean disregarding on-route changes, which can be demonstrated by the number of clustered trajectories ( Figure 1b). The more particles, the more vapor brought into TP. Thus, the North Atlantic contributes more than Indian Ocean when the on-route changes are disregarded.

Lines 95-101
When accounting for on-route changes (see Methods), we find long-term average RCs of ~42% from Asia, ~7% from the Indian Ocean, ~4% from North Africa, and ~3% from the North Atlantic (Fig.   1c). Disregarding on-route changes does not dramatically alter the outcomes (Fig. 1c). Crucially, during 98% of the study period (77% when disregarding on-route changes), the above four source regions account for over 75% of the total water vapor input to the TP (Fig. 1d-

Comment 4
Lines 135-137: 135 The TWS changes along the trajectories followed by water vapor transported towards the TP are influenced more prominently by the PME deficits in the North Atlantic than by both PME variations in the Indian Ocean and local PME variations. According to this sentence, you assessed the TWS over the continents along the trajectories, but, you are following the air parcels initially residing over the TP and not from the rest of Eurasia, so, the relationship between TWS and PME outside the TP deserves further analysis. Just explain it to me and improve the explanation in the manuscript if necessary.
Reply: We agree with you on this comment, and we analyzed the relations between TWS and PME over the 3 subregions which are outside the TP in Lines 121-136.
Besides, we also further dispersed the 3 subregions as 14 small-scale subregions, and analyzed the cross correlations between the TWS in detailed 14 sub-regions and PME over North Atlantic with considering optimum lags (see Lines 144-147).
In addition, by following your suggestion, we also performed the backward trace analyses for particles initially residing over the SR1 and SR2 and added discussion about it in this revision (see Lines 103-106). It demonstrates that the oceans (e.g., Indian Ocean and North Atlantic) recharge the water vapor of air particles over the lands (e.g., Asia and Africa) that further transit water vapor into the TP by the two main trajectories.
( Supplementary Fig. 1). The reason why we did not perform backward trace analysis for SR3 is that the SR3 covers most area of the TP for which backward trace analysis has already been done.

Lines 103-106
The oceans (e.g., Indian Ocean and North Atlantic) recharge the water vapor of air particles over the lands (e.g., Asia and Africa) that further transit water vapor into the TP (Supplementary Fig. 1 Supplementary Fig. 6).

Comment 5
Lines 182-183. Drying over the North Atlantic induced a TWS decrease between midlatitude Eurasia (20°N-50°N) and the southern TP mountains, but not in the central TP.
Regarding this sentence, I can understand the analysis for the TP taking into account that backward analysis was performed from this region, but for the rest of Eurasia, how can you explain this?
Reply: By following your suggestion, we have performed additional backward trace analysis for particles initially residing over the SR1 and SR2. Since the SR3 covers most area of the TP, thus, we did not perform additional backward trace analysis for the SR3. Additional backward trace analysis demonstrates that the oceans (e.g., Indian Ocean and North Atlantic) recharge the water vapor of air particles over lands (e.g., Asia and Africa) that further transit water vapor into the TP (Supplementary Fig. 1) by the two main trajectories. The first is from the North Atlantic to the western TP transported by the westerlies (Figure 1a-b), and the other is the intermedia current (see Figure 1a-b) that brings water vapor from North Atlantic and Indian Ocean into the southern TP.
In addition, based on the cross-correlation analyses (Supplementary Fig. 6 and Extended Data Fig. 3), it suggests that the PME deficit over the southeast North Atlantic mainly transited by westerlies leads to coherent declines in TWS across 14 subregions across the mid-latitude Eurasia and the southern TP mountains but not in TWS over the central TP.
Reply: This was a mistake from our side: we meant to refer to Extended Data Figure 5 (now is Extended Data Figure 4). We corrected the corresponding sentence as follows: Lines 184-186 in revised manuscript    Fig. 6-7).

Comment 8
The analysis of TWS under future scenarios seems very interesting.
Reply: Thanks for your encouraging comments!

Comment 9
Which version of FLEXPART are you using? You should mention it.
Reply: Thanks for your note. We are using FLEXPARTv10.4 in this study. And we have noted this information in the revised manuscript.

Comment 10
I wonder why the authors chose 500 000 particles? Is there any criterion for this? Can you explain what is the relation between the number of particles and the resolution of the input data? I have seen the world atmosphere divided into near 2 millions of parcels, but not a half of million for a small region like the TP; just clarify it to me and in the manuscript.
Reply: According to the guideline of the FLEXPARTv10.4 model 39 , there is no relation between the number of particles and the resolution of the input data 39 . On the other hand, the resolution of the input data depends on the data source 39 . Besides, the guideline points out that the choice of the number of particles depends on the scale of the target region and computation resource 39 . For example, the guideline suggests 900,000 particles for a small-scale simulation 39 . Besides, another regional study also applies thousands of particles for the determination of the moisture sources for the Western North American Monsoon 47 .
We have added a discussion regarding this in the revised manuscript (Lines 457-458).
Lines 457-458 The number of the particles depends on the scale of the target region and computation resources 39,47 .

Comment11
Lines 350-351. You argue that according to 39, … Given that it approximately takes 10 days for an air particle to travel around Earth. Reference 39 is about mining, so, from where did you take this metric? Secondly, I have read and used 10 days according to Numaguti (1999). There are also other studies that assess the residence time of the water vapor in the atmosphere, which is a better criterion than a particle travel (in days) around Earth.
Reply: Thank you for your insightful review. It is a kind of typo mistake. The original reference 39 is an example illustrating the bilinear interpolation of the GRACE data.
Following your suggestion, we improved the interpretation of the method by citing more appropriate reference here in the revised manuscript.

Lines 460-461 Revised interpretation
Given that the residence time of the water vapor in the atmosphere is approximate 10 days 48,49 , we perform backward simulations from day 1 to day 10 of each month.

Comment 13
Why 19 regions in the K-mean analysis? Which criterion did you utilized for determining the optimum number of clusters?
Reply: Based on the determination standard 56 for the optimal number of the clusters using k-means method, we initially determined the optimal number of the clusters for the trajectories is 3, and we further derived the 3 clusters of the trajectories from all trajectories ( Supplementary Fig. 24). It suggests that there are only 3 eastward trajectories into TP without denoting the trajectory from the Indian Ocean ( Supplementary Fig. 24). Although 3 is the optimal number of clusters, the result could not illustrate enough information about the transition path from the source regions to the TP. As for the criterion for determining the optimum number of clustered trajectories, we do not find relative standards in the handbook of the FLEXPARTv10.4 model and we also find various number of clusters used in other studies 47,51 . Thus, the optimum number of clusters is generally determined on a case-by-case basis.
In this study, we attempt to enrich the information of the trajectories and avoid far

Comment 14
In line 346 you describe the backward analysis for determining the moisture sources for the TP region. However, from line 359 to 364 you describe the high relationship between the total water vapor released from the TP and precipitation data from ERA5 and GPCC. So, it is not clear for me because in the backward analysis for identifying the moisture sources is calculated the moisture gain by the parcels and the total average on the total column. If you calculate the moisture release in the backward analysis, did you compute the moisture gain? It is not clear for me. Also, did you calculate the moisture release in the backward experiment computing PME within the TP and correlated it with the precipitation from ERA5 and GPCC averaged over the TP? Another doubt, did you run the (E-P) just in the backward mode?
Reply: In response to your first question, in the Equation 4, we have calculated the net moisture gain over source regions in the backward simulation and applied the net moisture gain >0 as the selection rule for particles. However, we did not directly determine the water vapor relative contribution from source regions by dividing the total vapor release over the TP with the moisture gain over the source regions. This is because that total moisture gain from all source regions is far larger than the total release over the TP ( Supplementary Fig. 12).
In this study, to determine the relative contributions from source regions, we firstly calculated the moisture gain over source regions in Equation 4 and  , and then determine the relative contributions of water vapor by dividing the total water vapor releases over TP by the release over the TP from every source region in Equation 6. Above is the calculation process when the on-route changes are accounted for. When the on-route changes are disregarded, the calculation process is the same as above but for all particles into TP.
We also calculated the relative contribution based on moisture gain from sources by following your suggestion. Noteworthy is that the relative contributions calculated based on (i) moisture gain from sources ( Supplementary Fig. 12) and (ii) moisture release over TP from sources ( Figure 1c) Fig. 12 and Figure 1c).
However, given the moisture gains over the source regions are far larger than the total water vapor release in the target region, thus we apply the first method in this study.
The second question in concern 14: Also, did you calculate the moisture release in the backward experiment computing PME within the TP and correlated it with the precipitation from ERA5 and GPCC averaged over the TP?
Reply: Thanks for the comments. Yes, we did calculate the moisture release in the backward experiment computing PME within the TP. In this study, we calculate the total water release for all particles residing over the TP when the (E-P) < 0. The original description "The total water vapor released from the TP is calculated by Eq. 2 (in this case, n=500,000)" is inappropriate, which might arouse misunderstandings. Thus, we revised the description in the manuscript accordingly in Lines 473-476.
Above calculation is based on previous study 52 , which demonstrated that "Eq. 2 can diagnose E-P, but not E or P individually". It also stated that "When rain falls, it normally clearly exceeds evaporation. Therefore, by assuming that E and P cannot coexist in the same location at the same time, instantaneous rates of evaporation Ei = E − P when E − P > 0, or precipitation Pi = P − E when E − P < 0, can be diagnosed.
(where i denotes time i)". Thus, we calculated the total moisture release for all particles residing over the TP when the (E-P) <0, and compared it with the GPCC and ERA5 precipitation data. We added relative reference to this description.

Lines 473-476 Revised description
The total water vapor release for particles residing over the TP is calculated by Eq. 2 when E-P<0 52 , and is compared to the standardized regional average precipitation obtained from ERA5 and GPCC.
The third question in concern 14: Another doubt, did you run the (E-P) just in the backward mode?
Reply: Yes. According to previous studies 39,47,51 , both forward-and backward-only simulations can successfully facilitate the determination of the water vapor source regions. In this study, we only run the FLEXPATv10.4 model for the determination of the moisture source regions for TP, SR1 and SR2 in the backward mode.

Comment 15
The moisture loses (E-P) <0 is not exactly the precipitation. Is has been named as 'Lagrangian precipitation' in previous studies (e.g. Nieto and Gimeno, 2019;Gimeno et al., 2020). In this study you intend to make a validation between the moisture lose computed with FLEXPART and the precipitation from ERA5 and GPCC, only through Pearson correlations, but it is not correct to name is as 'validation', you assess the relationship. Authors must clarify that moisture release is not the precipitation, and remove the word 'validation'. Indeed, for a real validation a correlation is not enough, as you should know.
Reply: We agree with you on this comment. The moisture loses (E-P) <0 over an area is not exactly the precipitation, but has positive effect on the formation of the precipitation 39,47,51 over the area. Thus, the positive relation between them could be a tool to evaluate the relative accuracy of the simulation results. If the relation is negative, the modelling results might be inaccurate. In our case, the total released water vapor over TP is highly related to both GPCC and ERA5 precipitation, with correlation coefficients between 0.85-0.86, which demonstrates the reliability of the simulated results here.
By firmly following your instruction, we have clarified that the moisture release is not the precipitation and removed the word "validation" in the revised manuscript (Lines 471-482).

Lines 471-482
A comparison is performed between the simulated total water vapor release and ERA5 and GPCC precipitation over TP to investigate the sensitivity of simulation results 51 .  Supplementary Fig. 11). We consider this to be as an acceptable accuracy for the FLEXPART results ( Supplementary Fig. 11).

Comment 16
In section 5 the authors describe the regression analysis performed. The TWS is considered as the total water terrestrial water resources, which consider the surface water bodies, so, I wonder why the authors did not make a lag in time analysis to determine the optimum lag in which TWS is affected by moisture contribution from the sources to precipitation over the TP. Indeed, in Supplementary Figure 9 is observed the low correlations.
Reply: Thanks for the suggestion. We performed the cross-correlation analyses between the TWS in 14 detailed subregions (Fig. 2b-c and Supplementary Fig. 6) across the mid-latitude Eurasia and the PME over the North Atlantic with consideration of the optimum lags (see Method 5). The optimum lag is defined as the lag when the maximum correlations between the TWS and the PME are detected. The interpretation of the cross-correlation analysis can be found in Lines 144-147.
Lines 144-147: Revised manuscript Focusing on regions within 20°N-50°N, considering optimum lags, we show high cross correlations between the PME deficit over the southeast North Atlantic and the TWS changes within 12 out of 14 sub-regions discerned along the water vapor propagation routes towards the TP (Fig. 2a-d, Supplementary Fig. 6).
We also analyzed the cross-correlation between TWS in southern TP mountains and PME over North Atlantic by considering the optimum lags (see Lines 161-164).

Lines 161-164
The cross-correlation analysis suggests the decreased TWS over the southern TP mountains can be attributed primarily to the PME deficit from the southeast North Atlantic with considering optimum lags .
Reply: Here, the ERA5 reanalysis data denotes the ERA5-based snow cover area over the region r during 2003-2016 (Lines 593-595).
Lines 593-595 in revised description

Comment 18
There is no conclusions?
Reply: Thanks for the suggestion. According to the "Guides for the Author" of Nature, there is no requirement on the Conclusion section but strict word limitation (less than 3200) on the Main Text. Since there are limited rooms for Conclusion section, we collectively presented the main findings of the study in the Abstract section. Thanks again for these highly constructive comments!  Hereafter, the TPM1-2 refers to southwest and southeast TP mountains, respectively. The TPS refers to the central TP surface and the NTP refers to northern TP. The continental map data 31 and the map data of TP 32 in panels a-b are from the public data source and plotted using R 33 . And results used to generate the figure are available in Zenodo 34 .

Supplementary Figure 4. Relations between PME over oceans and PME and TWS over SRs.
(a-c) denote spatial locations of the subregions for North Atlantic (a), Indian Ocean (b) and midlatitude Eurasia along water vapor trajectories to TP (c). (d-e) denote relations between PME over oceans and PME over SR1-3. (g-i) denote relations between PME over oceans and TWS over SR1-3. The continental world map data 31 in panels a-c and the map data of Tibet Plateau 32 in panel c are acquired from public data source and plotted by R 33 . And results used to generate the figure are available in Zenodo 34 .

Summary comments and reply
Review of "Sustainability of Asian water tower threatened by drying and warming ocean", by Zhang, Shen et al.
The authors present a very comprehensive study to possible explanations for the spatial variations in trends of terrestrial water storage on the Tibetan Plateau. They use a combination of the FLEXPART particle dispersion model with reanalysis climate data and GRACE remotely sensed TWS data, and argue that TWS decline in the southern TP can be attributed to a deficit in precipitation minus evapotranspiration from the southeast North Atlantic. They further argue that the high mountain ranges in Asia block the propagation of the deficit further onto the central TP which leads to increasing TWS there. Lastly, they demonstrate using CMIP6 GCM data how this blocking could weaken in the future, causing also the central TP to witness future TWS deficits.
To my knowledge this is an original approach, it is mostly well described and Reply: Thank you for your positive comments on our work. We have thoroughly checked the manuscript and revised it based on your comments and concerns. To make it easy to follow, we have included appropriate references to the revised manuscript (e.g., line numbers). Quoted text from the manuscript is provided in italics.

Comment 1
Title: The title seems somewhat strange in that it could be read as that the ocean is drying up. Consider rewording.
Reply: Thank you for your suggestion. We have modified the title as "Changes in oceanic climate threaten the sustainability of Asian water tower". Hopefully, the modified title better reflects the contents and is free of misunderstanding.

Supplementary Note 1
Based on the determination standard 56 for the optimal number of the clusters using kmeans method, we initially determined the optimal number of the clusters for the trajectories is 3, and we further derived the 3 clusters of the trajectories from all trajectories ( Supplementary Fig. 24). It suggests that there are only 3 eastward trajectories into TP without denoting the trajectory from the Indian Ocean ( Supplementary Fig. 24

Comment 4
Figure S1: Suggest to change the y-axis for panels b-d to be able to see more details, like also done in other rows in the figure.
Reply: Thank you for the comment. We changed the y-axis scale fitting all panels to display more details in Figure S1 and Figure S2.

Comment 5
L139-140: Looking at Figs S1 and S2, the contributions from North Africa are about equal to the North Atlantic. Therefore, I wonder how the authors exactly come to the conclusion that contributions from the North Atlantic dominate over contributions from North Africa? This conclusion is important as it propagates further in the analysis when focus is going to the North Atlantic region. I would guess that the ocean being a body of water supplies more water vapor than the northern African land mass with the dry Sahara, but I cannot deduct that from Figs S1-S2. In Figs S3-S9 correlations for the North Atlantic and Indian Ocean are shown and North Africa is not further investigated whereas the data in Figs S1 and S2 suggest it should. I am no expert in particle dispersion modeling, so the explanation could be something obvious, but then it requires a bit more explanation for the reader. Also, while Figs S3-S9 indeed show better correlations for North Atlantic compared to Indian Ocean, it seems that Indian Ocean still plays a large role, in particular IO4 (Fig S8, panel c).
Reply: We summarize the above questions as two questions and answer them accordingly.
Question 1: Why did not the study focus on the contribution of water vapor from the North Africa and other land source regions?
Answer: It is because that the ocean is the water vapor source into the continent by a large-scale atmospheric circulation. We performed additional backward trace analyses for the particles initially residing over the SR1 and SR2. Since SR3 covers most of area of TP, we did not perform further backward trace analysis for the SR3. According to Supplementary Fig. 1, it demonstrates that the oceans (e.g., Indian Ocean and North Atlantic) recharge the water vapor of air particles over the lands (e.g., Asia and Africa) that further transit water vapor into the TP via the two main trajectories. The first is from the North Atlantic to the western TP transported by the westerlies (Figure 1a-b), and the other is the intermedia current (see Figure 1a- Answer: The reason why the PME over Indian Ocean was not included in the further analysis is that the relations between TWS over SR1-3 and PME over IO4 is becoming weaker with the distance closer to TP (from SR1 to SR3). This does not apply to the directions of the two main transition paths (westerlies and intermedia current) into TP ( Figure 1b). Thus, we excluded the PME over Indian Ocean in the further analysis.
However, we found that the air temperature over the southeast TP mountains (TPM2) are coherent with the air temperature over the Indian Ocean (Extended Data Figure. 5).
Thus, in the further analysis, air temperature variations over the Indian Ocean have been included when we developed the models for the projection of the air temperature over the TPM2 (Equations. 18-19). The reason that you feel puzzled here might be the lack of the link between context. To address your concerns and avoid potential misunderstanding by readers, we have included the link at the end of the second section in the revised manuscript (Lines 176-181).

Lines 159-176
Our linear attribution analysis suggests that the increased TWS over the central TP can be attributed to decreased TWS across the southern TP mountains instead of increased PME (Extended Data Fig. 2). And the cross-correlation analysis suggests the decreased TWS over the southern TP mountains can be attributed primarily to the PME deficit from the southeast North Atlantic with considering optimum lags . Besides, significant relations between TWS depletion (e.g., glacier mass loss) and the reductions in snow cover area over southern TP mountains are detected ( Supplementary Fig. 7). The latter is mainly caused by regional warming ( Supplementary Fig. 8

). Even though glacier/snow meltwater from the southern TP mountains could increase TWS in the central TP 35 , such meltwater also replenishes TWS in the basins surrounding the TP. A comparative analysis for the surrounding basins (HSR8-12) of the TP and the central TP
shows that the surrounding regions are directly exposed to the westerlies and intermedia current that largely terminate at the southern margin of the central TP (Figure 2b). Accordingly, the TWS changes in surrounding basins  are consistent with the PME deficit transited by the westerlies from southeast North Atlantic but it is not the case for the TWS in the central TP (Extended Data Figures 3o).

Lines 176-181
Since the southern TP mountains belong to the HMA, it is reasonable to assume that blocking by the HMA's high topography dampens the propagation of the PME deficit that emerged in the southeast North Atlantic towards the TP, which in turn could have resulted in the observed increase in TWS over the central TP. We will further investigate this assumption in the next section based on the FLEXPART simulation.

Comment 7
L324-327: Why were these particular CMIP6 models included and others not? I also believe MICRO6 is a typo and should be MIROC6.
Reply: The selected CMIP6 models here include ACCESS-ESM1.5, BCC-CSM2-MR, CanESM5, GFDL-ESM4, IPSL-CM6A-LR, MIROC6, MRI-ESM2.0 and NorESM2-LM. The main reason why we select these particular CMIP6 models is that they are developed by the major organizations around the world and have been applied in the previous studies (Gillett et al., 2021;McKinnon et al., 2021). Besides, since the multiple average and weighting sum of CMIP6 outputs (e.g., PME and T) are finally determined for the future projections and organizations always provide several CMIP6 models by themselves, to avoid impacts on the weights in calculations of the weighting sum of CMIP6 outputs, we select only one CMIP6 model from each organization. Additionally, we selected the CMIP6 models under the r1i1p1f1 format. However, not all of CMIP6 models have the format, e.g., GISS-E2-1-G, HadGEM3-GC31-LL and CESM2. The r1i1p1f1 format is also one of the criteria for the selection of the CMIP6 models.
And the MICRO6 is a typo, we have corrected this typo through the whole manuscript. Besides, by following your suggestion, we compared the distributions of the CMIP6 outputs and the distribution for the ERA5 reanalysis data using the Kolmogorov-Smirnov test (Supplementary Fig. 16). It demonstrates that distributions of outputs from CMIP6 do not have significant bias from the distribution of the ERA5 reanalysis data when p<0.05. Given that the projected model is based on the CMIP6 output series, we determined both multiple mean and weighting sum of CMIP6 outputs to improve the capability of the CMIP6 models on reflecting the temporal variations in precipitation and temperature.

Comment 9
Supplementary Reply: Thanks for the note. We cited this relevant paper in the main text as the 20 th reference. Different from this relevant paper, we explored and clarified the mechanism behind the heterogeneities in the trends in terrestrial water storage (TWS) over the TP from the perspective of the large-scale atmospheric circulation. In our study, we evidenced that the blocking effect of the westerlies-carried changes in ocean climates by the High Mountain Asia as the cause for the abnormal increase in TWS. In addition, the blocking effect is demonstrated in weakening. Further, we also projected the area under TWS deficit in TP in future by considering the weakening blocking effect by the High Mountain Asia and oceanic climate changes.

Lines 73-75
The unprecedented changes in the HMA's water systems have been reported in numerous studies 3,8,19,20 . However, the atmospheric mechanisms causing the distinct and regionally-varying TWS changes are not well understood 5 .
L162-164: The correlation observed is clear, but I am not certain to which extent there is really a case of attribution. Could there be another mechanism that causes the correlation?
Reply: The comment was related to the sentence (L162-164) "Our linear attribution analysis suggests that the increased TWS over the central TP can be attributed to decreased TWS across the southern TP mountains instead of increased PME (Extended Data Fig. 2).". We agree that the observed correlation could also be due to other mechanisms (see our answer to the Question 2, Comment 4) and that the word "attribution" might have been too strong in this case. In the revised version, we thus choose a more careful wording (Lines 154-157): Lines 154-157 in the Track Changes version of the revised manuscript

Our linear attribution analysis suggests that the increased TWS over the central TP might be related to decreased TWS across the southern TP mountains instead of local increased PME
(Extended Data Fig. 2).

Comment 4
L179-182: How reasonable this assumption of blocking is remains rather vague to me. Is it not more likely that the TWS loss is mainly driven by increasing temperatures and the TWS in the central TP behaves differently because there are mostly endorheic basins and they do not drain the water out, in contrast to the surrounding areas, where loss of water stored in glaciers/snow drains out to the sea? I acknowledge there is not a straightforward answer and the exploration of the mechanism explored in the paper is very useful, I am just not completely convinced by the paper that this is indeed the driving mechanism. I suggest to use very careful wording that this could be one possible mechanism.
Reply: Thank you for these critical comments. We answer them on the basis of the two main questions here below.
Question 1: Is it not more likely that the TWS loss is mainly driven by increasing temperatures?
Reply: We concur that increasing temperatures in the southern TP mountains have an impact on the regional TWS, and indeed, we do consider "temperature" as one of the variables for predicting the future TWS in the region (see Equation 17). However, our analysis suggests that this is not the main driving factor. Although the warming induced a decrease in the seasonal snow cover ( Supplementary Fig. 8a,e), we noticed that TWS has only a poor relation to the temperature in the southern TP mountains (Supplementary Figure 8c,g). This is because by definition, TWS is the sum of the surface water, groundwater, soil moisture, canopy water and water stored in glacier ice. The reason why temperature is not the main driving factor is that glacier melt alone is not representative for the total TWS variations, thus explaining the poor correlations between TWS and temperatures in the southern TP mountains (Supplementary Figure 8c,g). This consideration also explains why we could not draw the conclusion that increasing temperatures mainly drive the decrease in TWS.
In our revision (Line 160-163), we better point out that the observed snow cover decrease is indeed induced by warming but that changes in TWS are less clearly affected by that (Supplementary We also detect a significant relation between TWS depletion (e.g., glacier mass loss) and the reductions in snow cover area over the southern TP mountains (Supplementary Fig. 7). In contrast to TWS declines, the reduction in snow cover area is mainly caused by regional warming ( Supplementary Fig. 8).
Question 2: Is it not more likely that the TWS in the central TP behaves differently because there are mostly endorheic basins and they do not drain the water out, in contrast to the surrounding areas, where loss of water stored in glaciers/snow drains out to the sea?
Reply: We concur with this possibility. However, some of our results are in contradiction to this.
Unlike the increased TWS observed over endorheic basins in the central TP (Figure 2b and Extended Data Figure 3o), the TWS over the Tarim basin in HSR12 (the largest endorheic basin in China) declines despite of the same meltwater supply from the southern TP mountains as for the central TP ( Figure 2b and Extended Data Figure 3l). This contrasting behavior suggests that the endorheic nature of the basins is not sufficient for explaining the observed increase in TWS. That said, we now provide a more balanced discussion about this, and instead of presenting a change in blocking effects as the sole hypothesis, we also discuss the potential for endorheic basins to play a role in the TWS variation. This additional discussion is found at Lines 165-184 of the revised manuscript: Lines 165-184 in the Track Changes version of the revised manuscript The endorheic character of the central TP differs from the exorheic nature of the basins situated in the south of the TP (e.g., Indus basin in HSR11, Figure 2b) Supplementary Figs. 9). Such a mechanism would also help explaining the correlation identified between (i) TWS changes in the northern TP and (ii) the variations in temperature and snow cover area in the southwest TP mountains . L216: The snow cover area cannot decline by more than 100%. This is probably related to the use of values of indices and confusing. Similar issue in L220-221. It is explained in the manuscript and Methods that these large percentual changes are because of use of indices, but still it is confusing. I suggest to think of an alternative way than percentual changes in index values, to indicate the changes.

Reference in this reply
Reply: We agree that values larger than 100% are confusing in this context. The reason for such values was that we originally determined the relative variation for a standardized index by dividing the variation in the index with the value at the beginning (Eq. 1 in the reply). In our revision, we Lines 220-222 in the Track Changes version of the revised manuscript Further, we detect an annual TWS decrease by ~11Gt (~66%, or ~13 kg/m 3 ) and ~12Gt (~55%, or ~12kg/m 3 ) in the southwest and southeast TP mountains during 2005-2016, respectively (Fig. 3b).

Lines 229-230 in the Track Changes version of the revised manuscript
Indeed, the area of the TP affected by a TWS deficit (i.e., the area with TWS ≤0) increases by ~167×10 4 km 2 (~62%) (Fig. 3b)  Reply: For the connection between snow cover decrease and blocking effects, please see our previous answer to Comment 5. For the specific sentence (now Lines 223-227) we adjusted the wording to: Lines 223-227 in the Track Changes version of the revised manuscript Indeed, the annual average snow cover areas decline by ~11× 10 4 km 2 (~70%) and ~7× 10 4 km 2 (~64%) in these two regions, respectively (Extended Data Fig. 6). We suggest that the reductions in snow cover area and TWS initiate the upslope wind in HMA and thus weakens the blocking effect by assisting air flows (i.e., westerlies) to cross the HMA.

Comment 8
Extended Data Figure 6: It is confusing that now the suffixes _Moun1 and _Moun2 are used where TPM1 and TPM2 were used before (and also in the caption).
Reply: Thank you for your professional suggestion. We now replaced "_Moun1" and "_Moun2" with "TPM1" and "TPM2", respectively, in the Extended Data Figure 6 and its legend. Reply: Yes, the corridor is defined as the path along which the area with negative TWS propagates from the southern TP mountains to the central TP. In the revision, we have denoted the corridor in the Figure 3e-h with arrows, and defined it in the legend of Figure 3. We also rephrased the description of the corridor in the revised manuscript (Lines 232-234): Lines 232-234 in the Track Changes version of the revised manuscript We also find corridors that act as paths for the area with negative TWS to transit from the southern TP mountains to the central TP (Figure 3e-h). These corridors roughly align with the westerlies or the Indian monsoon (Figure 3e-h).

Comment 10
L246-247: Here also I think specification is required for the reasoning why accelerated melting weakens the blocking effect.
Reply: For the connection between snow cover decrease and blocking effects, please see our