Aerosol forcing regulating recent decadal change of summer water vapor budget over the Tibetan Plateau

Domain: earth_science
Published: 2024-03-12
Editorial outcome: major_revision

Abstract
The Tibetan Plateau (TP), known as the Asian water tower, has been getting wetter since the 1970s. However, the primary drivers behind this phenomenon are still highly controversial. Here, we isolate the impacts of greenhouse gases (GHG), aerosols, natural forcings and internal climate variability on the decadal change of summer water vapor budget (WVB) over the TP using multi-model ensemble simulations. We show that an anomalous Rossby wave train in the upper troposphere travelling eastward from central Europe and equatorward temperature gradient in eastern China due to the inhomogeneous aerosol forcing in Eurasia jointly contribute to anomalous easterly winds over the eastern TP. Such anomalous easterly winds result in a significant decrease in water vapor export from the eastern boundary of the TP and dominate the enhanced summer WVB over the TP during 1979-2014. Our results highlight that spatial variation of aerosol forcing can be used as an important indicator to project future WVB over the TP.

Introduction
The Tibet Plateau (TP) and its surrounding high mountain areas are the second largest glacier-gathering area in the world except for polar ice sheets, storing a large amount of water resources in the form of glaciers, snow, lakes, and rivers1–3. Meanwhile, the TP is the birthplace of twelve important rivers in Asia, such as the Yangtze River, the Yellow River, the Yarlung Zangbo River, the Indus River, and the Ganges River. It provides necessary water resources for almost two billion people, and is known as the Asian water tower3–5. Therefore, change in water resources over the TP will seriously affect the natural ecosystems and socio-economic development in the TP itself and its surrounding countries.

The unique high-altitude terrain and atmospheric circulation dominated by monsoon and upper westerly bring abundant water resources to the TP6,7. However, climate change has greatly affected the hydroclimatic changes in the TP1,8–10. In the past decades, the TP has shown significant warming and humidification3,11–15. Meanwhile, the average precipitation over the TP has increased, mainly manifested as an increase in the northwest but a decrease in the southeast16,17. This phenomenon and the mechanism behind it have received widespread concerns from the scientific community and policymakers.

As one of the most active components of the water cycle, atmospheric water vapor transport determines the spatial pattern and trends of water resources in the Asian water tower4,18. Summer is the season with the highest water vapor budget (WVB) over the TP19. Previous studies reported a significant increasing trend in summer net WVB over the TP (especially in the north) since the late 1970s4,20. However, the dominant factors of increase in WVB over the TP are still highly controversial. On the one hand, enhanced atmospheric water content capacity driven by global warming may lead to the increase of water vapor content over the TP4,17. On the other hand, changes in water vapor transport driven by large-scale circulation are identified to be the main drivers of the interdecadal changes in WVB over the TP3. The weakening of the mid-latitude westerlies in the upper troposphere in Eurasia and the resulting decrease in water vapor export from the eastern boundary of the TP may dominate the interdecadal increase in net WVB over the TP20,21. Oceanic internal variability such as changes in sea surface temperatures in North Atlantic and Pacific may partially explain the westerly jet anomalies and changes in WVB in the northern TP16,20. However, an early study suggested that the high Asian mountains prevented further propagation of precipitation minus evaporation (PME) deficit originating from the southeastern North Atlantic to the central TP region, which in turn may have contributed to the wetting of north-central TP22.

As the second largest climate forcing factor, aerosol forcing has played a crucial role in regional climate change over the industrial era23–27. A recent study showed that the inhomogeneous changes in anthropogenic aerosol emissions in Eurasia since the 1970s (i.e., decreasing in Europe and increasing in Asia) was likely the main driver of the interdecadal weakening of the Eurasian subtropical westerly jet during summer28. This study raises the question of whether aerosol forcing is linked to the decadal changes of summer WVB over the TP. Here we examine the changes of summer WVB over the TP during the past four decades (1979–2014) using a climate reanalysis dataset. We disentangle the contributions of greenhouse gases (GHG), aerosols, natural forcings, and internal climate variability to the decadal trend of summer WVB over the TP using the multi-model simulations from the Coupled Model Intercomparison Project Phase 6 (CMIP6) (Supplementary Table 1), and identify the dominant role of aerosol forcing. Our results provide valuable insights into the mechanism of change in WVB over the TP.

Results
Observed changes in WVB over the TP
The water vapor sources over the TP mainly come from the water vapor transports of mid-latitude westerly winds in the upper troposphere in northern Hemisphere and Asian monsoon6,29. An amount of moisture is carried into the TP from its western, northern and southern boundaries by the westerly winds and Indian monsoon, respectively, while part of it is exported from its eastern boundary (Supplementary Fig. 1a). The ERA5 reanalysis shows that the summertime WVB (represented by the difference between precipitation and evaporation) over the entire TP is increasing significantly in the past four decades (1979–2014), especially in the western and northern regions, where the maximum increased rate exceeds 5 mm month−1 decade−1 (Fig. 1a). The increase of WVB has greatly contributed to the wetting trend of the TP in recent decades22,30,31 (Supplementary Fig. 2c), although there is a slight increase in evaporation due to global warming32 (Supplementary Fig. 2b).

The increasing summer WVB over the TP can be attributed to the changes of atmospheric moisture storage and moisture flux divergence (Details in Eq. 1 in “Method” section). During 1979–2014, the vertically integrated moisture flux divergence shows a significant increase (p < 0.1), with a seasonal average rate of 4.7 kg s−1 decade−1. In contrast, there are limited changes of atmospheric moisture storage, with a seasonal average rate of −0.05 kg s−1 decade−1. This comparison demonstrates the dominated contribution of moisture flux divergence to increasing WVB over the TP. Furthermore, we decompose the changes of vertically integrated moisture flux divergence to the changes of water vapor fluxes in four boundaries (Details in Eqs. 2–3 in “Method” section). The results show that there is a significant positive trend of 6 × 106 kg s−1 decade−1 (p < 0.1) in summer water vapor flux in the eastern boundary of the TP during 1979–2014 (Fig. 1b), indicating a marked decrease of water vapor export from this boundary. However, the summer water vapor input from the northern boundary significantly decrease at a rate of 2.8 × 106 kg s−1 decade−1 (p < 0.1) during 1979–2014. There are limited changes in water vapor inputs (not passing statistical significance at 90% confidence level) from the western and southern boundaries of the TP. Therefore, the decrease in water vapor export from the eastern boundary dominates the increase in net WVB over the TP in summer during 1979–2014.

Attribution of changes in water vapor flux in eastern boundary of the TP
The changes of summer WVB over the TP in the last four decades are likely influenced by multiple external forcings and internal climate variability33,34. Isolating these influences can only depend on climate models. The CMIP6 multi-model ensemble (MME) results show that the historical all forcing (ALL) simulations can reproduce the observed positive trend of water vapor flux in the eastern boundary of the TP in the past four decades (Supplementary Table 2). Notably, the CMIP6 models underestimate the increase of summer water vapor flux in the eastern boundary of the TP, with a MME averaged trend of 0.8 kg s−1 decade−1. This could be primarily related to the systematic underestimation of climatological vertically integrated water vapor and WVB over the TP and its surrounding areas by the climate models34,35 (Supplementary Fig. 1b vs. 1a).

Furthermore, we quantify the impacts of multiple external forcings and internal climate variability on the enhanced summer water vapor flux in the eastern boundary of the TP during 1979–2014 using ALL, single-forcing (AER, GHG, and NAT) and piControl simulations (Fig. 2). Model results reveal that the aerosol forcing dominates the enhanced water vapor flux in the eastern boundary of the TP in the past four decades (red point in Fig. 2) in the historical all forcing simulations, with trends of 1.1 vs. 0.8 kg s−1 decade−1 in the MME. In addition, internal climate variability also contributes to a wetter TP (orange point in Fig. 2), but the value is less than half of that from aerosol forcing. In a contrary, greenhouse gas and natural forcings show a negative contribution to the enhanced water vapor flux in eastern boundary of the TP (blue and green points in Fig. 2), which partly offset the impacts of aerosol forcing and internal climate variability. Our attribution analysis is generally consistent with two earlier studies33,34, which conclude that the recent decadal changes of WVB over the TP are mainly attributed to combined effects of anthropogenic forcings and internal climate variability. Nonetheless, our study further isolates the impacts of greenhouse gas and aerosol forcings and overweight the dominated contribution from aerosol forcing rather than internal climate variability to enhanced WVB over the TP.

Potential mechanisms
The changes of water vapor flux are determined by both wind speed and specific humidity, representing dynamic and thermodynamic processes, respectively. By fixing specific humidity (wind speed) and varying wind speed (specific humidity), we decompose anomalous summer water vapor flux in the eastern boundary of the TP to the contributions from the changes in dynamic and thermodynamic processes (Fig. 3). It is seen that the changes of zonal wind in the upper troposphere (i.e., anomalous easterly winds) dominate the decrease in summer water vapor export from the eastern boundary during 1979–2014 (red bar vs. gray bar in Fig. 3), which primarily contributes to the increase in summer WVB over the TP. In a contrary, the changes of specific humidity slightly decrease the summer water vapor flux in the eastern boundary of the TP (blue bar in Fig. 3a). The dominated contribution from dynamic process to enhanced summer water vapor flux in the eastern boundary of the TP is also captured by climate models from CMIP6, despite differences in amplitude (Fig. 3b).

By calculating the linear trends of global summer winds at 500 hPa, we further examine the possible mechanisms of the interdecadal change in anomalous easterly winds driven by external forcings, such as greenhouse gases and aerosols (Fig. 4). Observation shows an anomalous wave train in the upper troposphere propagating eastward from central Europe to East Asia during 1979–2014 (Supplementary Fig. 3). This wave train leads to significant positive anomalies in geopotential height at 500 hPa near Lake Baikal, indicating an anomalous anticyclone here. Anomalous easterly winds have increased markedly in the south of the anomalous anticyclone, which contribute to the decrease in water vapor export from the eastern boundary of the TP caused by the climatological westerly winds and increase in water vapor over the TP. Such anomalous anticyclone near Lake Baikal and easterly winds over the TP are clear in ALL and AER simulations (Fig. 4a, c), but not seen in GHG and NAT simulations (Fig. 4b, d). The anomalous wave train in the upper troposphere in Eurasia is also reproduced in AER simulation (Fig. 4c), highlighting the leading role of aerosols in regulating interdecadal change in easterly anomalies.

In the past decades, the anthropogenic aerosol emissions in Europe show large decline to meet the European Union’s target for human health36,37. From 1979 to 2014, the aerosol optical depth (AOD) in Europe (0–40°E, 40–60°N) decreases at a rate of 0.07 per decade in ALL simulation (Fig. 5a), causing a large disturbance on regional climate change26,38. The reduced AOD greatly increase the surface solar radiation (SSR) (Fig. 5b), resulting in a local warming in the lower atmosphere (Fig. 5c). The diabatic heating can induce anomalous ascending motion and then excite an anomalous Rossby wave train in the upper troposphere, which propagates energy eastward from Europe (Fig. 4). As a consequence, an anomalous anticyclone circulation at 500 hPa is formed near Lake Baikal, which continuously weakens the westerlies in the eastern boundary of the TP, leading to weakened water vapor export from the eastern boundary of the TP in the past four decades. Therefore, the decrease in water vapor export from the eastern boundary driven by atmospheric circulation anomaly dominates the increase in summer precipitation over the TP. The more latent heat release from increasing precipitation may strengthen convection development, which may in turn intensify anomalous anticyclone in the upper troposphere of the northern TP39, thus forming a local positive feedback effect. Such an anomalous Rossby wave train driven by aerosol reductions in Europe was also reported in an early study40.

Additionally, the anomalous anticyclone circulation can cause local warming in the upper-middle troposphere northeastern of TP (Fig. 5d). Meanwhile, a local cooling is found to the east of TP, which is mainly attributed to increased aerosol emissions in eastern China during 1979–2014 (Fig. 5a). Thus, an anomalous temperature gradient towards equator forms in the upper troposphere over the eastern China (Supplementary Fig. 4a), which further weakens the westerly winds in the eastern boundary of the TP based on the thermal wind relationship (Supplementary Fig. 4b). Such an anomalous temperature gradient and weakening westerly over TP are reproduced again in ALL and AER simulations (Fig. 6a, b, e, f). In a contrary, the westerly over the TP is increased slightly in NAT simulations due to the enhanced poleward temperature gradient (Fig. 6g, h), posing a negative effect on enhanced water vapor flux in the eastern boundary of the TP in the past few decades. Noticeably, although slightly weakened westerlies in the eastern boundary of the TP due to the anomalous temperature gradient towards equator is shown in GHG simulations (Fig. 6c, d), GHG forcing also shows a negative effect on enhanced water vapor flux in the eastern boundary of the TP (Fig. 2). Such discrepancy is attributed to the impact of thermodynamic process rather than dynamic process in GHG simulations (Supplementary Fig. 5).

Methods
Observation and model outputs
We choose the fifth-generation climate reanalysis of the European Centre for Medium-Range Weather Forecast (ERA5) as benchmark to evaluate the performance of climate models from the CMIP6 and quantify the impacts of various factors on the WVB, due to the lack of long-term continuous multivariate observations in the TP. Because of its high reliability, the ERA5 reanalysis has been widely used in climate change research across regions, including the TP46–48. Monthly variables, including specific humidity, zonal and meridional winds, surface pressure, and geopotential height, during 1979–2014 are used to investigate the decadal change of WVB in the TP.

Considering the availability of all required variables in this study, eleven climate models in CMIP6 are selected to quantify the contributions of anthropogenic and natural forcings to the decadal change of WVB over the TP (Supplementary Table 1). These models provide both historical all forcings simulations (ALL) and single-forcing simulations, including GHG forcing only (GHG), anthropogenic aerosol forcing only (AER), and natural forcing only (NAT, solar, and volcanic combined) from the Detection and Attribution Model Intercomparison Project (DAMIP). For a specific scenario, the outputs of the first three ensemble member (r1i1p1f, r2i1p1f, and r3i1p1f) in each model are used in this study. Moreover, 44 chunks of nonoverlapping 36-year timeseries from preindustrial control simulations (piControl) with ten models, except CNRM-CM6-1 are used to estimate the impacts of internal climate variability on the decadal change of WVB in the TP. Both model outputs and observations are re-gridded to a common resolution of 1.5°×1.5° using the bilinear interpolation method. In this study, we first calculate each model average with three ensemble members, which is then used to generate the multi-model mean.

Calculation of WVB
In this study, we calculate the vertically integrated WVB at monthly time-scale based on Trenberth and Guillemot49:

<P-E> = - <∂q/∂t> - <∇·(qV)> + R (1)

Where P, E, q, V and R represent precipitation, evaporation, specific humidity, horizontal wind, and residual term, respectively. < > represents the vertical integration from the surface to top of atmosphere. The first term of - <∂q/∂t> represents the time change of atmospheric moisture storage, which can be negligible at monthly time-scale. The second term of - <∇·(qV)> represents the vertically integrated moisture flux divergence, which can be calculated through water vapor fluxes at four boundaries of the TP based on the 2D divergence theorem (Supplementary Fig. 9):

∫∫ <∇·(qV)> dA = ∮ <qV> dl (2)

Where A represents the area of the TP; l represents the length of boundaries. Here, we further calculate the net WVB across the boundaries of TP as follows:

∮ <qV> dl = B_W + B_S - B_N - B_E (3)

Where B_E, B_W, B_S and B_N (kg s−1) present the water vapor flux in eastern, western, southern, and northern boundaries, respectively, which can be calculated as follows:

B_E = ∫ Q_u dl
B_W = ∫ Q_u dl
B_S = ∫ Q_v dl
B_N = ∫ Q_v dl (4)

Where Q_u (kg m−1 s−1) and Q_v (kg m−1 s−1) represent the vertically integrated zonal and meridional water vapor flux, respectively, which can be calculated as follows:

Q_u = 1/g ∫_p_t^p_s q u dp
Q_v = 1/g ∫_p_t^p_s q v dp (5)

Where g represents gravitational acceleration; q represents specific humidity; u represents zonal wind vector; v represents meridional wind vector; p represents atmospheric pressure; p_s represents surface pressure; p_t represents pressure at upper bound (setting to 50 hPa).

Data availability
The monthly meteorological variables from the ERA5 reanalysis are open access at https://www.ecmwf.int/en/forecasts/datasets/reanalysis-datasets/era5. The multi-model outputs from historical, hist-aer, hist-GHG, hist-nat, and piControl experiments can be obtained at https://esgf-node.llnl.gov/search/cmip6/.

Code availability
Data analysis in this study was conducted using open-source programming languages and software, including Climate Data Operators (CDO; https://code.mpimet.mpg.de/projects/cdo) and NCAR Command Language (NCL; http://www.ncl.ucar.edu/). All figures were produced using NCL.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "In this study, the authors address two important topics: climate change over the Tibetan Plateau (TP) and\naerosol forcing. ERA5 reanalysis data and several CMIP6 models were used to idenfify the reason for the\nincreased wefting of the TP in recent decades. It was found that a significant decrease in the flux of\nmoisture from the eastern boundary of the TP is the main contributor to the recent increase in total\ncolumn water vapor over TP. This decrease, in turn, is related to a circulafion anomaly induced by\nasymmetric aerosol forcings over Europe and East Asia. The conclusions are interesfing. But I hope they\nwould have been more convincing. My foremost reason for worry is the methodology. The study uses\nreanalysis and model-derived data sets at pressure level over a complex topography where the data\nitself could be erroneous due to inter/extrapolafion, and the effect of verfical flux cannot be ignored.\nSecondly, it is not shown if the moisture budget is closed. Without a closed budget, it is difficult to\nconclude if a parficular component is responsible for an observed change. My detailed comments are\nbelow.\nMajor comments\n1. The moisture budget at the fime scale of a season of longer usually has very small components of the\nstorage term (TPW). This insignificant term is at the center of this study. An increase in TPW with year\ndue to flux, that is a couple of orders of magnitude higher than the d(TPW)/dt must be jusfified through\na budget that is closed.\n2. The study region is covered by complex topography. Over such regions. Moisture flux through\nmountain slopes is important. Can pressure level data capture such fluxes? Can it be shown that those\nfluxes are not important?\n3. As a confinuafion of the previous comment, the role of evaporafion is not described. Due to global\nwarming, the retreat of glaciers and the increase of water bodies at that alfitude could increase\nevaporafion. What effect the surface evaporafion have on the budget of moisture and its trend?\n4. It is not obvious if the large-scale circulafion anomaly (Fig 1a, 3) is due to east-west asymmetric\naerosol forcing or due to a change in rainfall over the TP itself. A change in rainfall and related diabafic\nheafing could change the circulafion paftern and reduce the outgoing flux of moisture through the\neastern boundary on account of increased convergence.\nMinor comments:\n1. It is not clear what the analysis is done for which season. Is it the annual average of the boreal\nsummer average? It’s not menfioned anywhere in the text (am I missing something?).\n2. It is also not clear why the study is restricted from 1979 to 2014. The ERA5 data is available from 1940\nup to recent fime and a longer data set should provide more robust results, along with long CMIP6 runs.\n3. In Fig 4 (as well as Supp Fig S6), it is menfioned that the data is averaged between 90-120E. That\nmeans this dimension is reduced. As yet, the x-axis ficks are longitude.", "reviewer_2": "The authors present an invesfigafion into the driving causes behind decadal scale changes in the\nsummerfime water vapor budget over the Tibetan plateau. They use the standard driver aftribufion\nmethod of analysing the changes in relevant physical proxies in single forcer climate model simulafions,\nand conclude that there is an outsized influence from the changes in anthropogenic aerosol emissions\nsince 1979.\nTibetan Plateau moisture budgets is clearly an important topic, and understanding the underlying drivers\nof recent changes is crucial. Overall, the authors seem to have performed a thorough analysis, and\nreached a number of insighfful conclusions. I do have a few methodological concerns, however. Also,\nimportantly, the presentafion of the results makes it difficult to fully follow the logic of the paper. I\ntherefore recommend major revisions before this paper can be considered for publicafion.\nAlso, I note that a very similar paper appeared in Nature very recently. While the authors were likely\nunaware of this work, it would be important to include, discuss and compare to the results of Jiang et al.\nin the revised manuscript:\nhftps://www.nature.com/arficles/s41586-023-06619-y\nMethodological points:\n- L125: I am very surprised that the authors chose to exclude the two CMIP6 models that show a\ndiffering trend to observafions. Since the objecfive of this study is to disentangle the effects of various\ndrivers, these models would be key for understanding whether one parficular driver has a strong\ninfluence. As an example, the high climate sensifivity of CanESM5 could be driving a strong moisture\ninfluence due to CO2 changes. Excluding models \"to improve model confidence level\" (L128) is not\ndefensible in this case.\n- The authors use three ensemble members from each model in the single forcer cases. I know this is at\nfimes all that is available, but for weak changes like aerosol influces on a regional scale it is clear that\neven three ensemble members may be too liftle. Is there one or more models where there are more\nmembers available, where the authors could test whether three members can be expected to give a\nreasonable representafion of the forced signal?\n- Would it be relevant to include a process schemafic that illustrates what the different drivers do to\nmoisture transport into/out of the TP, and why? This would, I think, make the dynamical explanafions\nprovided much easier to follow.\nPresentafion points:\n- My main concern here is that so much of the key material resides in supplementary figures, that\nfollowing the logic of the paper is virtually impossible without them. Also, the figures are not really\ndiscussed in order, but are brought into the discussion at various fimes and without clear explanafion.\nI strongly encourage the authors to expand the number of figures in the main paper, and to slightly\nrestructure the presentafion so that it links closer to the graphics. Examples are Figures S4 and S5, which\nare key to understanding the forcing that would be driving the changes described.\n- The current main figures also need units on the axes, and improved clarity. Notably figure 3 is more or\nless impossible to read, as the green arrows are hard to separate.\nMinor points:\n- L94: \"A mount\" -> \"An amount\"?\n- L96: \"Observafions indicate\" Which observafions? (It's stated in the Methods, but please repeat it here\nfor the benefit of the reader.)\n- L114: THis comment on another paper should perhaps be in the discussions?\n- L143: \"0.1 vs 0.13\" These numbers seem idenfical, unless they have high precision. What is the\nuncertainty?\n- L246: \"see-saw like changes\": By only using the one period, 1979-2014, the authors do not really\ninvesfigate any see-saw like behaviour, I think? It underlies the results, but the trends shown would wash\nout any influence of a european emissions increase followed by a decrease?\n- L261: \"the problemtafic emissions inventory\" Please explain what is meant here, rather than just\nreferring to other studies.\n- L262: \"designing specific model experiments\" These now exist:\nhftps://gmd.copernicus.org/arficles/16/4451/2023/gmd-16-4451-2023.html"}, "author_response": "We are grateful to the reviewers for your precious time in providing helpful comments\nand guidance that have improved the manuscript. In this document, we describe how\nwe have addressed the reviewer’s comments. Referee comments are shown in black\nand author responses are shown in blue.\nReviewer #1 (Remarks to the Author):\nIn this study, the authors address two important topics: climate change over the Tibetan\nPlateau (TP) and aerosol forcing. ERA5 reanalysis data and several CMIP6 models\nwere used to identify the reason for the increased wetting of the TP in recent decades.\nIt was found that a significant decrease in the flux of moisture from the eastern\nboundary of the TP is the main contributor to the recent increase in total column water\nvapor over TP. This decrease, in turn, is related to a circulation anomaly induced by\nasymmetric aerosol forcings over Europe and East Asia. The conclusions are interesting.\nBut I hope they would have been more convincing. My foremost reason for worry is\nthe methodology. The study uses reanalysis and model-derived data sets at pressure\nlevel over a complex topography where the data itself could be erroneous due to\ninter/extrapolation, and the effect of vertical flux cannot be ignored. Secondly, it is not\nshown if the moisture budget is closed. Without a closed budget, it is difficult to\nconclude if a particular component is responsible for an observed change. My detailed\ncomments are below.\n⚫ Thank you very much for the positive and constructive comments and suggestions.\nAll the questions and concerns have been carefully answered and the paper has\nbeen revised accordingly. For your concerns about complex topography impacts\nand closed moisture budget, please refer to our detailed responses shown below:\nMajor comments\n1. The moisture budget at the time scale of a season of longer usually has very small\ncomponents of the storage term (TPW). This insignificant term is at the center of this\nstudy. An increase in TPW with year due to flux, that is a couple of orders of magnitude\nhigher than the d(TPW)/dt must be justified through a budget that is closed.\nResponse: Thanks for your comments! We are sorry that our Figure 1a (vertically\nintegrated precipitable water vapor) may mislead you. We would like to clarify that we\nfocused on moisture flux divergence term (− 〈𝛻(𝑞𝑽)〉) not storage term (−〈 𝜕𝑞/𝜕𝑡 〉) in this\nstudy. In the revised paper, we replaced Figure 1a with the difference (𝑃−𝐸) between\nprecipitation and evaporation representing the net water entering a region and improved\nthe descriptions of the methodology section as follows:\n“Calculation of WVB\nIn this study, we calculate the vertically integrated WVB at monthly time-scale based\non Trenberth and Guillemot 49:\n𝑃−𝐸 = −〈 𝜕𝑞/𝜕𝑡 〉 − 〈𝛻(𝑞𝑽)〉+𝑅 (1)\nWhere 𝑃, 𝐸, 𝑞, 𝑽 and 𝑅 represent precipitation, evaporation, specific humidity,\nhorizontal wind and residual term, respectively. 〈 〉 represents the vertical integration\nfrom the surface to top of atmosphere. The first term of −〈 𝜕𝑞/𝜕𝑡 〉 represents the time\nchange of atmospheric moisture storage, which can be negligible at monthly time-scale.\nThe second term of − 〈𝛻(𝑞𝑽)〉 represents the vertically integrated moisture flux\ndivergence, which can be calculated through water vapor fluxes at four boundaries of\nthe TP based on the 2D divergence theorem (Figure S10):\n∬− 〈𝛻(𝑞𝑽)〉𝑑𝐴 = ∫−〈𝑞𝑽〉𝑑𝑙 (2)\nWhere 𝐴 represents the area of the TP; 𝑙 represents the length of boundaries. Here,\nwe further calculate the net WVB across the boundaries of TP as follows:\n∫−〈𝑞𝑉〉𝑑𝑙 = 𝐵𝐸 +𝐵𝑊 −𝐵𝑆 −𝐵𝑁 (3)\nWhere 𝐵𝐸 , 𝐵𝑊 , 𝐵𝑆 and 𝐵𝑁 (kg s-1) present the water vapor flux in eastern, western,\nsouthern and northern boundaries, respectively, which can be calculated as follows:\n𝐵𝐸 = ∫𝑄𝑢 𝑑𝑙𝐸\n𝐵𝑊 = ∫𝑄𝑢 𝑑𝑙𝑊\n𝐵𝑆 = ∫𝑄𝑣 𝑑𝑙𝑆\n𝐵𝑁 = ∫𝑄𝑣 𝑑𝑙𝑁\n(4)\nWhere 𝑄𝑢 (kg m-1 s-1) and 𝑄𝑣 (kg m-1 s-1) represent the vertically integrated zonal\nand meridional water vapor flux, respectively, which can be calculated as follows:\n𝑄𝑢 = − 1/𝑔 ∫ 𝑝𝑡 𝑝𝑠 𝑞𝑢 𝑑𝑝\n𝑄𝑣 = − 1/𝑔 ∫ 𝑝𝑡 𝑝𝑠 𝑞𝑣 𝑑𝑝\n(5)\nWhere 𝑔 represents gravitational acceleration; 𝑞 represents specific humidity; 𝑢\nrepresents zonal wind vector; 𝑣 represents meridional wind vector; 𝑝 represents\natmospheric pressure; 𝑝𝑠 represents surface pressure; 𝑝𝑡 represents pressure at upper\nbound (setting to 50 hPa).” Lines 328-351\nFigure 1. Changes of summer water vapor budget. (a) Spatial pattern of linear trends\nin the difference (𝑃− 𝐸) between precipitation (𝑃) and evaporation (𝐸) (shade, unit:\nmm month-1 decade-1) and water vapor flux (vector, unit: kg m-1 s-1 decade-1) over the\nTibetan Plateau during 1979-2014. The black points and green vectors indicate that the\ntrends of precipitable water vapor and water vapor flux pass the test at 0.1 confidence\nlevel. (b) Time series (point) and linear fitting (line) of water vapor flux (unit: 106 kg s-\n1) in four boundaries of the Tibetan Plateau during 1979-2014. The asterisks represent\nstatistical significance at 90% confidence level. The numbers represent linear trends\n(unit: 106 kg s-1 decade-1) of water vapor flux in each boundary of the Tibetan Plateau.\n2. The study region is covered by complex topography. Over such regions. Moisture\nflux through mountain slopes is important. Can pressure level data capture such fluxes?\nCan it be shown that those fluxes are not important?\nResponse: Thanks for your comments! We agree that moisture flux through mountain\nslopes is important in complex topography. However, we would like to clarify that we\ndiagnosed vertically integrated moisture flux divergence over the TP through\ncomputing fluxes at the boundaries based on the 2D divergence theorem instead of\ncentered finite differences for the rectilinear latitude-longitude grids (Please refer to the\nabove reply for details). In four boundaries of the TP, the impacts of complex\ntopography are limited. In the revised paper, we also clarified as follows:\n“It is noted that previous studies revealed that complex topography can influence the\ndiagnosis of moisture transport over the TP 50. However, there are limited impacts of\ncomplex topography on WVB in this study, because we select a line integral method\nbased on the 2D divergence theorem to diagnose the net WVB over the TP.” Lines 352-\n355\n3. As a continuation of the previous comment, the role of evaporation is not described.\nDue to global warming, the retreat of glaciers and the increase of water bodies at that\naltitude could increase evaporation. What effect the surface evaporation have on the\nbudget of moisture and its trend?\nResponse: As responded to the first comment, our study mainly focused on the\ndifference (𝑃−𝐸) between precipitation and evaporation. According to the suggestion,\nwe have analyzed the trend of evaporation and added the descriptions as follows:\n“The increase of WVB has greatly contributed to the wetting trend of the TP in recent\ndecades 22, 30, 31 (Figure S2c), although there is a slight increase in evaporation due to\nglobal warming 32 (Figure S2b).” Lines 100-102\n(a)P (b) E\nmm month-1 decade-1 mm month-1 decade-1\n(c)P - E\nmm month-1 decade-1\nFigure S2 Observed linear trends of summer mean precipitation (𝑃, a), evaporation (𝐸,\nb) and its difference (𝑃−𝐸, c) during 1979-2014. The black points represent statistical\nsignificance at 90% confidence level.\n4. It is not obvious if the large-scale circulation anomaly (Fig 1a, 3) is due to east-west\nasymmetric aerosol forcing or due to a change in rainfall over the TP itself. A change\nin rainfall and related diabatic heating could change the circulation pattern and reduce\nthe outgoing flux of moisture through the eastern boundary on account of increased\nconvergence.\nResponse: Thanks a lot for this constructive comment. As you noted, the diabatic\nheating associated with the variations in precipitation can modulate large-scale\ncirculation anomalies. However, it is worth noting that in this study the circulation\nanomalies responsible for the precipitation anomalies belongs to a large-scale\nanomalous anticyclone centered to the northeast of the TP with distance more than 3000\nkm. It is difficult to attribute the formation of the anticyclone of this spatial scale to\nheating anomalies associated with TP precipitation anomalies. Furthermore, both the\nanticyclone and TP precipitation anomalies are atmospheric anomalies, which cannot\nmaintain longer than a month without external forcing factors. However, we also\nacknowledge that a change in rainfall over the TP caused by anomalous circulation and\nrelated diabatic heating release may in turn intensify the circulation anomaly. So, we\nadded the descriptions in the revised paper as follows:\n“Therefore, the decrease in water vapor export from the eastern boundary driven by\natmospheric circulation anomaly dominates the increase in summer precipitation over\nthe TP. The more latent heat release from increasing precipitation may strengthen\nconvection development, which may in turn intensify anomalous anticyclone in the\nupper troposphere of the northern TP 39, thus forming a local positive feedback effect.”\nLines 192-197\nMinor comments:\n1. It is not clear what the analysis is done for which season. Is it the annual average of\nthe boreal summer average? It’s not mentioned anywhere in the text (am I missing\nsomething?).\nResponse: Our study focused on summer mean water vapor budget over the TP. It’s\nmentioned in title, abstract and text in lines 82, 97, 105, 128, 135, 156, 167……\nAdditionally, we have added the “summer mean” in appropriate places of the text\nthrough careful checking.\n2. It is also not clear why the study is restricted from 1979 to 2014. The ERA5 data is\navailable from 1940 up to recent time and a longer data set should provide more robust\nresults, along with long CMIP6 runs.\nResponse: Thanks for your comments! This study focused on the period of 1979-2014\ndue to the following reason:\nThe accuracy of reanalysis before 1979 remains to be debated, especially in areas with\nrare observations such as the Tibet Plateau due to the lack of satellite data assimilation.\nThat's also why most reanalysis datasets (e.g., MERRA-2, NCEP-2 and JRA-55) started\nin 1979.\n3. In Fig 4 (as well as Supp Fig S6), it is mentioned that the data is averaged between\n90-120E. That means this dimension is reduced. As yet, the x-axis ticks are longitude.\nResponse: Sorry for our mistakes in Fig.4 and Fig.S6. The x-axis ticks should be\nlatitude and we have corrected it in the revised paper (new Figures 6 and S9).\nReviewer #2 (Remarks to the Author):\nThe authors present an investigation into the driving causes behind decadal scale\nchanges in the summertime water vapor budget over the Tibetan plateau. They use the\nstandard driver attribution method of analyzing the changes in relevant physical proxies\nin single forcer climate model simulations, and conclude that there is an outsized\ninfluence from the changes in anthropogenic aerosol emissions since 1979. Tibetan\nPlateau moisture budgets is clearly an important topic, and understanding the\nunderlying drivers of recent changes is crucial. Overall, the authors seem to have\nperformed a thorough analysis, and reached a number of insightful conclusions. I do\nhave a few methodological concerns, however. Also, importantly, the presentation of\nthe results makes it difficult to fully follow the logic of the paper. I therefore recommend\nmajor revisions before this paper can be considered for publication.\n⚫ Thank you very much for the positive and helpful evaluations. All the questions\nand concerns have been carefully answered and the paper has been revised\naccordingly.\nAlso, I note that a very similar paper appeared in Nature very recently. While the\nauthors were likely unaware of this work, it would be important to include, discuss and\ncompare to the results of Jiang et al. in the revised manuscript:\nhttps://www.nature.com/articles/s41586-023-06619-y\nResponse: Accepted. We have compared and discussed the results of Jiang et al. (2023)\nas follows in the revised paper:\n“Meanwhile, compared with a recent report by Jiang et al. 41, which showed that the\ndecrease in meridional troposphere temperature gradient over Eurasia caused by the\nuneven anthropogenic aerosol emissions led to the weakening of the Asian subtropical\nwesterly jet, which ultimately dominated the increase in summer precipitation in the\nnorthern TP since the 1950s, we suggests a new mechanism, i.e., a teleconnection\npathway across Eurasia induced by European aerosol forcing. We further reveal that\nthe above mechanism and anomalous temperature gradient in eastern China due to the\nuneven aerosol forcing in Eurasia collectively play a dominant role in the weakening\nof the westerly jet over the eastern TP. The weakened westerly winds decrease the water\nvapor output from the eastern boundary of the TP, rather than other boundaries, which\nmay be the primary cause of increase in precipitation in the northern TP.” Lines 239-\n- L125: I am very surprised that the authors chose to exclude the two CMIP6 models\nthat show a differing trend to observations. Since the objective of this study is to\ndisentangle the effects of various drivers, these models would be key for understanding\nwhether one particular driver has a strong influence. As an example, the high climate\nsensitivity of CanESM5 could be driving a strong moisture influence due to CO2\nchanges. Excluding models \"to improve model confidence level\" (L128) is not\ndefensible in this case.\nResponse: Thanks for your comments!\n(i) We agree that different models have different climate sensitivities to a specific\nforcing, e.g., greenhouse gases and aerosol. However, we first need to ensure\nconsistency between historical simulations and observation for a specific study, e.g.,\nthe increasing water vapor flux in the eastern boundary of the TP in the past four\ndecades. We usually consider the mechanism analysis to be reliable, only when the\nhistorical simulation for a specific model is consistent with observation. In this study,\ntwo climate models showed an opposite trend of water vapor flux compared to\nobservation, which may be attributed to a lower climate sensitivity of aerosol and a\nhigher climate sensitivity of greenhouse gases. Such climate model selection based on\nobservation and historical simulations were also common in many previous studies.\n(ii) We also analyzed the results of above two models and discussed in the revised paper\nas follows:\n“We select nine climate models to further analyze the possible mechanisms behind the\nrecent decadal change of summer WVB because CanESM5 and NorESM2-LM\nsimulate a trend of water vapor flux in the eastern boundary of the TP contrary to the\nobservation. However, it is noted that both CanESM5 and NorESM2-LM still show the\nlargest positive trend of summer WVB in the eastern boundary of the TP in the response\nto aerosol forcing (Figure S7).” Lines 257-262\nFigure S7 Attribution of changes in summer water vapor budget in the eastern\nboundary of the Tibetan Plateau from NorESM2-LM and CanESM5. The gray, red,\nblue and green points represent simulated linear trends (unit: 106 kg s-1 decade-1) of\nwater vapor budget in eastern boundary of the Tibetan Plateau during 1979-2014 from\nall (ALL), aerosol-only (AER), greenhouse gas-only (GHG) and natural-only (NAT)\nforcings simulations, respectively.\n- The authors use three ensemble members from each model in the single forcer cases.\nI know this is at times all that is available, but for weak changes like aerosol influences\non a regional scale it is clear that even three ensemble members may be too little. Is\nthere one or more models where there are more members available, where the authors\ncould test whether three members can be expected to give a reasonable representation\nof the forced signal?\nResponse: As suggested, we added the new Figure S6 to compare the results of ten\nensemble members from IPSL-CM6A-LR and MIROC6. The conclusion that aerosol\nforcing regulate the recent decadal change of summer water vapor budget over the\nTibetan Plateau remains the same in the results from ten ensemble members from two\nmodels. Also, we added these descriptions in the revised paper:\n“Although the only first three ensemble members in each model are used in this study,\nthe main conclusion that aerosol forcing regulate the recent decadal change of summer\nwater vapor budget over the Tibetan Plateau remains the same in ten ensemble members\nfrom IPSL-CM6A-LR and MIROC6 (Figure S6).” Lines 254-257\nFigure S6 Attribution of changes in summer water vapor budget in the eastern\nboundary of the Tibetan Plateau from ten ensemble members from IPSL-CM6A-LR\nand MIROC6. The gray, red, blue and green points represent simulated linear trends\n(unit: 106 kg s-1 decade-1) of water vapor budget in eastern boundary of the Tibetan\nPlateau during 1979-2014 from all (ALL), aerosol-only (AER), greenhouse gas-only\n(GHG) and natural-only (NAT) forcings simulations, respectively. The error bars\nrepresent one standard deviation of ten ensemble members for each model.\n- Would it be relevant to include a process schematic that illustrates what the different\ndrivers do to moisture transport into/out of the TP, and why? This would, I think, make\nthe dynamical explanations provided much easier to follow.\nResponse: As suggested, we added the schematic diagram Figure 7 to illustrating the\npaths of uneven aerosol forcing in Eurasia regulating the recent decadal change of\nsummer water vapor budget over the TP in the revised paper.\nFigure 7. Schematic diagram illustrating the paths of uneven aerosol forcing in Eurasia\nregulating the recent decadal change of summer water vapor budget over the TP.\nPresentation points:\n- My main concern here is that so much of the key material resides in supplementary\nfigures, that following the logic of the paper is virtually impossible without them. Also,\nthe figures are not really discussed in order, but are brought into the discussion at\nvarious times and without clear explanation. I strongly encourage the authors to expand\nthe number of figures in the main paper, and to slightly restructure the presentation so\nthat it links closer to the graphics. Examples are Figures S4 and S5, which are key to\nunderstanding the forcing that would be driving the changes described.\nResponse: As suggested, we removed some key figures from Supplementary to main\ntext, including Figure S2, Figure S4 and Figure S5. Additionally, we adjusted the order\nof some tables and figures to read logically.\n- The current main figures also need units on the axes, and improved clarity. Notably\nfigure 3 is more or less impossible to read, as the green arrows are hard to separate.\nResponse: As suggested, we added the units on the axes for all figures in main text and\nSupplementary. Additionally, we revised the figure 3 (new figure 4) to show arrows\nclearly.\nMinor points:\n- L94: \"A mount\" -> \"An amount\"?\nResponse: Corrected.\n- L96: \"Observations indicate\" Which observations? (It's stated in the Methods, but\nplease repeat it here for the benefit of the reader.)\nResponse: As suggested, we revised “observations” as “The ERA5 reanalysis” Line 96\n- L114: TH is comment on another paper should perhaps be in the discussions?\nResponse: As suggested, we removed this comment to discussion section. Lines 235-\n239\n- L143: \"0.1 vs 0.13\" These numbers seem identical, unless they have high precision.\nWhat is the uncertainty?\nResponse: As suggested, we revised the units from 106 kg s-1 year-1 to 106 kg s-1 decade-\n1. This comparison has changed to \"1.0 vs 1.3\" in the revised paper. Line 140\n- L246: \"see-saw like changes\": By only using the one period, 1979-2014, the authors\ndo not really investigate any see-saw like behaviour, I think? It underlies the results,\nbut the trends shown would wash out any influence of a european emissions increase\nfollowed by a decrease?\nResponse: In this study, “see-saw like” represents spatial pattern (decreasing\nanthropogenic emission in Europe but increasing anthropogenic emission in China) not\ntime variation. To avoid misleading readers, we have revised “see-saw like aerosol\nforcing” to “uneven or inhomogeneous aerosol forcing”.\n- L261: \"the problemtatic emissions inventory\" Please explain what is meant here,\nrather than just referring to other studies.\nResponse: As suggested, we clarified as follows in the revised paper:\n“However, the CMIP6 models fail to capture the trends in AOD and SSR in China\nduring this time due to an opposite trend of anthropogenic aerosol emissions in the\nCommunity Emissions Data System (CEDS) inventory 42, 43.” Lines 288-291\n- L262: \"designing specific model experiments\" These now exist:\nhttps://gmd.copernicus.org/articles/16/4451/2023/gmd-16-4451-2023.html\nResponse: Thanks very much for this useful information. We have added the\ndescriptions as follows in the revised paper:\n“It will be worthwhile to quantify the impact of reductions in aerosol emissions over\nChina in the future by using specific model experiments such as the recent Regional\nAerosol Model Intercomparison Project (RAMIP) 44.” Lines 291-294"}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors have made several changes to the manuscript in response to both reviewers. However, I sfill\nfeel that some of the major concerns raised by me are sfill unresolved.\n1. My first concern (Comment 1) was on the relafive magnitude (contribufion) of the tendency term of\nthe moisture budget equafion. In response, the authors have modified the analysis to include the P-E\nanalysis instead of the advecfion of moisture. However, this does not answer my query as (P-E) includes\n<dq/dt> in addifion to the advecfion term. If this tendency term is large, how can it be concluded that\nadvecfion is responsible for changes in P-E? In summary, without closing the moisture budget, it is a bit\ndifficult to make a profound statement on the causality of changes in P-E.\n2. In addifion, the equafion of moisture budget needs to be corrected (in both the main and reply to the\nreviewers). The term ( ) is missing the dot product. Please correct.\n〈𝛻 𝑞𝑽 〉\n3. I’m not sure how one can confirm that the effect of boundary slopes is limited in these calculafions of\ngradient (in response to my Comment 2). Can it be shown, at least quofing some previous study?", "reviewer_2": "I thank the authors for their thorough revision. The updated manuscript is now much clearer, and I\nappreciate the addifion of the schemafic.\nWhile I am broadly convinced by the arguments presented, I sfill have some concerns regarding the\nchoice of models and ensemble members. The new figures S6 and S7 help, but sfill do not provide\nsufficient informafion to assure that the aftribufion presented in Figure 2 - which is key for the\nconclusions of the paper - is fully supported.\nAn example of the issue is the CanESM5 result in Figure S7. Its net flux is very close to zero, however the\nauthors have chosen to exclude the model since its overall flux is negafive. The results for the individual\nensemble members are not shown, but it is very likely that at least one of them is net posifive - and\nshould therefore, by the authors methodology, have been included if the analysis was one on a per-\nmember basis.\nExcluding models based on the ensemble mean could sfill be defensible, however Figure S6, which\nshows 10-member results for two models, highlights the issue. Here, both models have a posifive net,\nbut with a 1 std.dev. spread that touches zero. Hence, there will be members in each sample that are\nnegafive. Consequently, it is fully possible to select a 3-member set that gives a negafive mean flux, even\nfor a model where the 10-member mean is posifive.\nSo, how can the authors be certain that this is not the case for the two excluded models? The answer is\nthat they can't, at least based on a 3-member mean.\nThe easiest solufion here is to include all models in the main result. Based on figure S6, this is unlikely to\naffect the overall conclusions, but it will give a less biased representafion of errors. At the same fime, I\nurge the authors to include the (3-member mean) results for all individual models in Figure 2 (or in a\ncorresponding supplement). This gives a much clearer indicafion of how robust the modelling/aftribufion\npart of the study is.\nIf the authors sfill wish to exclude two models, I believe that a more thorough analysis is needed to show\nthat these two really have robustly different behavior to the others.\nAnother suggesfion that would help make the argument clearer is to also show the results for all\nmembers in figure S7. Standard deviafions do not fully show the spread of the members, or how robust\nthe mean behavior is.\nI thank the authors again for their confinued efforts, and hope that these final comments are taken as\nconstrucfive."}, "author_response": "We are grateful to the reviewers for your precious time in providing further helpful\ncomments and guidance that have improved the manuscript. In this document, we\ndescribe how we have addressed the reviewer’s comments. Referee comments are\nshown in black and author responses are shown in blue.\nReviewer #1 (Remarks to the Author):\nThe authors have made several changes to the manuscript in response to both reviewers.\nHowever, I still feel that some of the major concerns raised by me are still unresolved.\n⚫ Thanks for your positive evaluations for our revised version.\n1. My first concern (Comment 1) was on the relative magnitude (contribution) of the\ntendency term of the moisture budget equation. In response, the authors have modified\nthe analysis to include the P-E analysis instead of the advection of moisture. However,\nthis does not answer my query as (P-E) includes <dq/dt> in addition to the advection\nterm. If this tendency term is large, how can it be concluded that advection is\nresponsible for changes in P-E? In summary, without closing the moisture budget, it is\na bit difficult to make a profound statement on the causality of changes in P-E.\nResponse: Thanks for your further comments! As suggested, we compared relative\nmagnitude of linear trends between atmospheric moisture storage and moisture flux\ndivergence. The related contents are added in revised paper as follows:\n“The increasing summer WVB over the TP can be attributed to the changes of\natmospheric moisture storage and moisture flux divergence (Details in Equation 1 in\nMethod section). During 1979-2014, the vertically integrated moisture flux divergence\nshows a significant increase (p<0.1), with a seasonal average rate of 4.7 kg s-1 decade-\n1. In contrast, there are limited changes of atmospheric moisture storage (−\n𝜕〈𝑞〉/𝜕𝑡\n), with\na seasonal average rate of -0.05 kg s-1 decade-1. This comparison demonstrates the\ndominated contribution of moisture flux divergence to increasing WVB over the TP.”\nLines 104-110\n2. In addition, the equation of moisture budget needs to be corrected (in both the main\nand reply to the reviewers). The term −∇∙〈𝑞𝑽〉 is missing the dot product. Please\ncorrect.\nResponse: Thanks! The term −∇〈𝑞𝑽〉 has been corrected as −∇∙〈𝑞𝑽〉 in the revised\npaper.\n3. I’m not sure how one can confirm that the effect of boundary slopes is limited in\nthese calculations of gradient (in response to my Comment 2). Can it be shown, at least\nquoting some previous study?\nResponse: Thanks for your further comments! In previous versions, we used the\nsimplified moisture flux divergence term −〈∇∙(𝑞𝑽)〉 instead of −∇∙〈𝑞𝑽〉. Based on\nthe Equation 5b in Trenberth and Guillemot (1995), −∇∙〈𝑞𝑽〉 can be written as −∇∙\n〈𝑞𝑽〉 = −〈∇∙(𝑞𝑽)〉−𝑞𝑠 𝑽𝑠 ∙∇𝑝𝑠 . The second term −𝑞𝑠 𝑽𝑠 ∙∇𝑝𝑠 represents the effects\nof complex topography. As suggested, we calculated the effects of complex topography\non the trend of moisture flux divergence (Figure R1). The results show that complex\ntopography has a small effect on increasing moisture flux divergence over the TP during\n1979-2014.\nFurthermore, we calculated moisture flux divergence term through water vapor fluxes\nat four boundaries of the TP based on 2D divergence theorem ∬𝑑𝑖𝑣 𝑭 𝑑𝐴 = ∫𝑭𝑑𝑙,\nthat is ∬−𝛻∙〈𝑞𝑽〉𝑑𝐴 = ∫−〈𝑞𝑽〉𝑑𝑙 (Equation 2 in Main text). Indeed, the effects of\ncomplex topography are included in calculation of moisture flux divergence term\nthrough water vapor fluxes at four boundaries.\nReference: Trenberth, K. E. & Guillemot C. J. Evaluation of the global atmospheric\nmoisture budget as seen from analyses. Journal of Climate 8, 2255-2272 (1995).\n(a) (b)\nmm/month/decade mm/month/decade\n(c)\nmm/month/decade\nFigure R1 The panels (a), (b) and (c) represent linear trends of −∇∙〈𝑞𝑽〉, −〈∇∙(𝑞𝑽)〉\nand −𝑞𝑠 𝑽𝑠 ∙∇𝑝𝑠 during 1979-2014, respectively.\nReviewer #2 (Remarks to the Author):\nI thank the authors for their thorough revision. The updated manuscript is now much\nclearer, and I appreciate the addition of the schematic.\n⚫ Thanks for your positive evaluations for our revised version.\n1. While I am broadly convinced by the arguments presented, I still have some concerns\nregarding the choice of models and ensemble members. The new figures S6 and S7\nhelp, but still do not provide sufficient information to assure that the attribution\npresented in Figure 2 - which is key for the conclusions of the paper - is fully supported.\nAn example of the issue is the CanESM5 result in Figure S7. Its net flux is very close\nto zero, however the authors have chosen to exclude the model since its overall flux is\nnegative. The results for the individual ensemble members are not shown, but it is very\nlikely that at least one of them is net positive - and should therefore, by the authors\nmethodology, have been included if the analysis was one on a per-member basis.\nExcluding models based on the ensemble mean could still be defensible, however\nFigure S6, which shows 10-member results for two models, highlights the issue. Here,\nboth models have a positive net, but with a 1 std.dev. spread that touches zero. Hence,\nthere will be members in each sample that are negative. Consequently, it is fully\npossible to select a 3-member set that gives a negative mean flux, even for a model\nwhere the 10-member mean is positive. So, how can the authors be certain that this is\nnot the case for the two excluded models? The answer is that they can't, at least based\non a 3-member mean. The easiest solution here is to include all models in the main\nresult. Based on figure S6, this is unlikely to affect the overall conclusions, but it will\ngive a less biased representation of errors. At the same time, I urge the authors to include\nthe (3-member mean) results for all individual models in Figure 2 (or in a corresponding\nsupplement). This gives a much clearer indication of how robust the\nmodelling/attribution part of the study is. If the authors still wish to exclude two models,\nI believe that a more thorough analysis is needed to show that these two really have\nrobustly different behavior to the others.\nResponse: We really agreed the further comments from Reviewer. As suggested, we\nused the simulations from all climate models to attribute and analyze the possible\nmechanisms behind increasing water vapor budget over the TP during 1979-2014. We\nhave carefully updated the figures and related texts in the revised paper. Furthermore,\nfollowing suggestions, we have revised the Figure 2 and included the results from all\nmodels in the revised paper as follows:\nFigure 2. Attribution of changes in summer water vapor flux in the eastern\nboundary of the Tibetan Plateau. The gray, red, blue, green and orange markers\nrepresent simulated linear trends (unit: 106 kg s-1 decade-1) of water vapor flux in\neastern boundary of the Tibetan Plateau during 1979-2014 from all (ALL), aerosol-only\n(AER), greenhouse gas-only (GHG), natural-only (NAT) forcings and piControl (PI)\nsimulations in eleven climate models, respectively. The points represent the ensemble\nmeans of eleven climate models. The error bars represent one standard deviation of\neleven climate models.\n2. Another suggestion that would help make the argument clearer is to also show the\nresults for all members in figure S7. Standard deviations do not fully show the spread\nof the members, or how robust the mean behavior is.\nResponse: As suggested, we added the results for all members in Figure S7 (New\nFigure S6) in the revised paper as follows:\n(a) IPSL-CM6A-LR (b) MIROC6\n1-e )\nd\na\nc\ne\nd\n1-s\ng\nk\n6\n0\n1\n(\nd\nn\ne\nTr\nFigure S6 Attribution of changes in summer water vapor flux in the eastern boundary\nof the Tibetan Plateau from ten ensemble members from IPSL-CM6A-LR and\nMIROC6. The gray, red, blue and green markers represent simulated linear trends (unit:\n106 kg s-1 decade-1) of water vapor flux in eastern boundary of the Tibetan Plateau\nduring 1979-2014 from all (ALL), aerosol-only (AER), greenhouse gas-only (GHG)\nand natural-only (NAT) forcings simulations in ten ensemble members, respectively.\nThe points and error bars represent ensemble means and one standard deviations of ten\nensemble members for each model.\nI thank the authors again for their continued efforts, and hope that these final comments\nare taken as constructive.\n⚫ Thanks for your constructive comments again."}, {"round": 3, "reviewer_comments": {"reviewer_1": "The authors have provided reasonable responses to the second round of comments by both reviewers.\nAt this point, I don't have any new queries and I congratulate the authors for a great job done.", "reviewer_2": "I thank the authors for this second revision. My concerns have been addressed in the revised figures, so I\nhave no further comments."}, "author_response": ""}]

Decision letter