South Asian black carbon is threatening the water sustainability of the Asian Water Tower

Long-range transport of black carbon from South Asia to the Tibetan plateau and its deposition on glaciers directly enhances glacier melt. Here we find South Asian black carbon also has an indirect effect on the plateau’s glaciers shrinkage by acting to reduce the water supply over the southern Tibetan plateau. Black carbon enhances vertical convection and cloud condensation, which results in water vapor depletion over the Indian subcontinent that is the main moisture flux source for the southern Tibetan plateau. Increasing concentrations of black carbon causes a decrease in summer precipitation over the southern Tibetan plateau, resulting in 11.0% glacier deficit mass balance on average from 2007 to 2016; this loss rises to 22.1% in the Himalayas. The direct (accelerated melt) and indirect (mass supply decrease) effects of black carbon are driving the glacial mass decline of the so-called “Asian Water Tower”.


Introduction
We describe the methods in more details in Section 1, including the accumulative anomaly methodology, model evaluation statistics, WRF-Chem model set up with an in-depth simulation evaluation and uncertainty analysis, as well as a monthly-scale glacial mass balance model. Next, we discuss precipitation and incoming moisture variations over the Tibetan Plateau, including precipitation trends over the Tibetan Plateau and changes in incoming moisture from South Asia to the southern Tibetan Plateau (Section 2). South Asian black carbon emissions, loadings, and atmospheric heating, are described in Section 3. To better explain mechanistically the relationship between black carbon and the change in the water vapor transport to the Tibetan Plateau, we have added a calculation of moisture flux and its divergence from 1989 to 201 (Section 4). Lastly, in Section 5, we have added additional analyses to support the key conclusions with respect to glacial retreat over the Tibetan Plateau; this includes (i) a reconstruction time series of reference-surface mass balance from 1979 to 2014 using an empirical model, (ii) calculated contributions of summer precipitation and temperature change to glacier mass change, as well as (iii) its spatial heterogeneity.
An itemized summary of the datasets used in this study is provided in Section 6.

The accumulative anomaly method
The accumulative anomaly method was used to reveal the inflexion point (year) where the change in precipitation occurred "due to" black carbon emission accumulation. The accumulative anomaly is an index to distinguish the changing tendency of discrete data 1 . For a discrete series, the accumulative anomaly (Xt) for data point can be expressed as: where is the mean value of the series xi, and n is the number of discrete points.
A positive accumulative anomaly indicates that the corresponding data point is higher than the average, otherwise lower than the average. If the relation curve is composed of least two parts as noted, then the inflexion point/year can be identified. Here, the variable x represents summer precipitation, or black carbon emission.

Model evaluation statistics
We examined the skill of model performance utilizing in-situ observations, satellite data, and a reanalysis dataset vs. model output. Point-based propagation was assessed for precipitation and black carbon. The following statistical metrics 2,3 were used in the study to evaluate the suitability of different physics and chemical schemes for a wide range of scenarios and applications.
The bias was calculated as the difference between the WRF-Chem simulations, the observed precipitation, and black carbon. The mean error (ME, equation (2)) was derived using all sites in the study area.
Mean absolute error (MAE, equation (3)) was calculated to quantify the difference between the simulations and the observational data while negating the effect of cancelling positive and negative errors seen in the bias.
Root mean square error (RMSE, equation (4)), which is a measure of the combination of systematic and random error, was calculated from the simulated results.
Pearson's correlation coefficient (R, equation (5) TMP (equation (7)) is a combined index derived from MAE, RESM, and the slope to quantify the overall performance so that each simulation ensemble member could be compared individually.
= ( + ̅ + (1 + R)) + |1+ | 4 (7) Here, the metric was calculated, where n, yi and xi are the number of observations, simulation value, and observation value, respectively. is the mean of observed values and is the mean from the simulations.

WRF-Chem setup
The WRF-Chem model is a newly developed regional dynamic/chemical transport model 4 which simulates gas-phase chemical and aerosol microphysical processes, along with numerous meteorological fields 5

Evaluation of precipitation
To assess the ability of the WRF-Chem model to reproduce the summer precipitation spatial distribution, we compared simulated output with other gridded datasets including GPCP (Global Precipitation Climatology Project) 19 , CRU (Climatic Research Units) 20 , and ERA5 (ECMWF Reanalysis 5th Generation) 21   were compared with WRF-Chem precipitation estimates.
Supplementary Fig. 3 The terrain height and location of in-situ observations for summer precipitation observations (mark with black dots) in the study area. Next, we investigated the sensitivity of each physics scheme to scrutinize bias in the simulated results; this was done by grouping all ensemble members together with a common physical option and comparing their respective bias range using box and whisker plots ( Supplementary Fig. 4). The analysis was conducted in the study area using the summer precipitation results from the WRF-Chem simulation along with station observations. The ensemble members with the Kain-Fristch scheme gave the least bias over both South Asia ( Supplementary Fig. 4a) and the Tibetan Plateau

Evaluation of black carbon
In previous work, we had confirmed that the CBMZ gas-phase chemical mechanism and MOSAIC aerosol module were the most suitable in the reconstruction of black carbon (BC) concentrations over the Tibetan Plateau and adjacent regions 25,26 .
As is shown in Supplementary An additional reanalysis dataset (the Modern-Era Retrospective analysis for Research and Applications version 2 (MERRA-2)) was also used as a comparative for WRF-Chem simulation output; this can be seen in Supplementary Fig. 6. Both

Model uncertainty analysis
In addition to different physical and chemical schemes, variations in spatial resolutions can also result in simulation uncertainties. In this study, we conducted the WRF-Chem simulations at a 25 km spatial resolution. To analyze model uncertainty caused by the spatial resolution, we designed a higher resolution WRF-Chem simulation, i.e., at a resolution of 8 km -the result can be seen in Supplementary Fig.   12. A comparison of both model resolution runs of the 500 hPa wind field with the ERA-Interim reanalysis is plotted in Supplementary Fig. 12a

Glacier mass balance model
The monthly-scale mass balance model from Radić and Hock 56 was used to model glacier mass balance for the Tibetan Plateau for the period 1979-2014. The primary input data in this model includes monthly precipitation and air temperature; this was obtained from the CRU dataset. Glacier area-weighted specific mass balance (B) for the whole glacier in each mountain range was calculated as a sum of the specific mass balance (b) of each elevation band on a glacier (i): where bi and Si denote specific mass balance and glacier area, respectively, subscript (i) represents the number of elevation bands on a glacier (i = 1, 2, 3, …, n) with an elevation interval of 50 m.
Monthly specific mass balance (bi, mm water equivalent) for each elevation band on a glacier was calculated as: where ai is glacier surface ablation (negative), ci is glacier mass accumulation (positive), and Ri refers to snowmelt refreezing (positive) at each elevation band.
For debris-free glaciers, the glacier surface ablation (ai, mm water equivalent) was calculated based upon a positive degree-day model 57 in which snow/ice melt is considered linearly correlated to monthly air temperature. Here, monthly ai was computed as: where fsnow/ice denotes the degree-day factor for snow/ice (mm water equivalent day −1 °C −1 ), and Ti denotes monthly air temperature (°C) above the glacier surface.
For debris-covered glaciers, the relationship between debris thickness and glacier surface ablation 58 was used to calculate the effect of debris cover on glacier surface ablation. In detail, an averaged curve for debris thickness and surface ablation was used in the mass balance model to calculate the reduction of ice ablation due to debris cover. To calculate the ablation of debris-covered glaciers, the ablation factor ki was introduced to quantify the effect of the debris on glacier surface ablation. Thus, monthly glacier surface ablation due to debris cover (ai, debris) was calculated as: where ki denotes the scale factor for each elevation band (i) and depends on the debris thickness.
Monthly glacier mass accumulation for each elevation band ci (mm water equivalent) was calculated as: where δm is a constant, Ti denotes air temperature at each elevation band, and Tsnow represents the threshold temperature differential between snow and rainfall. If Ti is below Tsnow, Pi is assumed to be snow. Otherwise, Pi is regarded as liquid precipitation.
Based on the relationship between annual potential refreezing Ri,pot (cm) and the air temperature Ta (°C) at each elevation band 40 , Ri,pot was calculated by equation (13): ,pot = −0.69 × + 0.0096 (13) where the lower boundary of annual snowmelt refreezing over the whole glacier is Monthly air temperature at each elevation band (Ti) was calculated as: where TCRU denotes the CRU TS monthly air temperature for the period from 1979 to 2014, and h represents the altitude of the glacier elevation band.
To scale up CRU TS monthly precipitation to hmax, a precipitation correction factor (kp) was assigned, while a precipitation gradient (dpre) was used to interpolate precipitation to each elevation band (the percentage of precipitation decreases with every 50 m decrease in elevation) from the highest to lowest elevation of the glacier.
Thus, monthly precipitation for each elevation band was calculated as: where PCRU denotes CRU TS monthly grid precipitation from 1979-2014 at each cell location on the glacier.
The glacier mass balance sensitivities to temperature change (∆ /∆ ) and precipitation change (∆ /∆ ) were calculated as: where B is the simulated annual average mass balance (m water equivalent year −1 ) for the period from 1979-2014, B (∆ ) is the same B but in response to a step-wise temperature change in the range from −6 K to +6 K incremented by 0.5 K, and B (∆ ) is the same B, but in response to a step-wise precipitation change from −30% P to +30% P with incremented by 5% P: these were advanced in an uniform stepwise change to the original monthly temperature and precipitation time series throughout the 36-year period.

Precipitation trends over the Tibetan Plateau
We utilized the GPCP monthly precipitation dataset to investigate the precipitation trend over the Tibetan Plateau. The GPCP is a monthly mean rainfall dataset. GPCP is a merged dataset on a 2.5-degree global grid derived from rain gauge station data, satellite derived data, and sounding observations for the period 1979 to the present. As shown in Supplementary Fig. 15, the GPCP precipitation trend corroborates the increased summer precipitation trend (derived from the CRU dataset) over the southern Tibetan plateau and, more importantly confirms the trend reversal at the beginning of 21st century. Based upon our earlier analysis of CRU monthly gridded precipitation data, it had been established that, in the northern Tibetan plateau, there has been an increase in summer precipitation over the past 56 years.

Supplementary Fig. 15 Summer precipitation trend over the Tibetan Plateau (TP).
Linear trends in summer area-averaged precipitation over the northern TP using the Climatic  Fig. 16a). Subsequently, we calculated a time subset (i.e., the 2004-2016 period using CRU) for the summer precipitation trends over the southern Tibetan plateau and adjacent regions as plotted in Supplementary Fig. 16b. As expected, the CRU analysis revealed a decreasing trend in the summer precipitation was evident over the southern Tibetan plateau, while the opposite was the case for many regions within Indian subcontinent.

Plateau
In the southern Tibetan plateau, the impact of localized surface evaporation is not particularly pronounced 59 , thus the inter-annual variability of the summer precipitation is primarily controlled by "long-range" moisture transport. To investigate moisture transport and moisture change from South Asia to the southern Tibetan Plateau, we computed the moisture flux and its divergence using the European Center

South Asian black carbon emissions, loadings, and atmospheric heating
To investigate black carbon emissions over South Asia, we utilized monthly global black carbon emission inventories split out by sector (energy production, industry, transportation, residential and commercial, agriculture, as well as deforestation and wildfire), at 0.1°×0.1° resolution; these data were obtained from Peking University. As shown in Fig. 19, the temporal trend and variation of black dataset was analyzed. It was discovered that atmospheric black carbon concentrations have been ever-increasing with time ( Supplementary Fig. 20a to Supplementary Fig.   20c). We also investigated the wet removal of black carbon during the summer monsoon, which have been also ever-increasing from 1986 to 2016 ( Supplementary   Fig. 20d to Supplementary Fig. 20f).
To examine the impact of South Asian black carbon on summer air temperatures, the WRF-Chem regional model was tuned and optimized to run simulations for each year from 2007 to 2016 ( Supplementary Fig. 21). The results of the WRF-Chem simulations showed that increased atmospheric black carbon loading does indeed result in an increasing temperature change in the atmosphere.

South Asian black carbon linkage with to summer precipitation decrease in south Tibetan Plateau
In the main text, we identified a decreasing trend in the summer precipitation over the southern Tibetan Plateau at an inflexion point around 2004; this decline, it was discovered, is regulated by the long-range moisture transport into the region.
Further analysis was conducted to examine the change of long-range moisture transport into the southern Tibetan Plateau and, it was found that a reduction had

The lag correlation analysis between black carbon and summer precipitation
As was plotted in Fig. 1  To establish as to whether a delay / tipping point exists between the two, we conducted a lag correlation analysis between South Asian black carbon and summer precipitation over the southern Tibetan Plateau. As can be seen in Supplementary Fig.  22, there is an insignificant (P＞0.05) lag relationship between summer precipitation over the Tibetan Plateau and South Asian black carbon; this is especially so for a lag time of more than 10 years. Thus, we concluded that the negative correlation is more than likely due to black carbon reaching a critical level in the atmosphere to affect local convection and precipitation. Returning to Supplementary Fig. 19 once again, i.e., the difference in South Asian black carbon emission compared to the global average, we remarked upon the larger divergence between the two since the 21 th century. We also referenced Ramanathan et al 62 . to infer that the divergence may well induce significantly higher radiative perturbations than the globally mean estimates.

Mechanism of South Asian black carbon affecting on precipitation
Using the WRF-Chem model, a series of sensitivity simulations were undertaken to examine the climate dynamics and forcing due to atmospheric black carbon over Another central factor that our analysis uncovered is that of convective available potential energy (CAPE). Supplementary Fig. 4c in the main text, confirms that there has been a significant increase in CAPE in the eastern Indian subcontinent. Intensified black carbon concentrations which is an absorbing aerosol, induces positive radiative forcing and subsequent heating of the low and middle troposphere ( Supplementary   Fig. 26b). Moreover, air temperature and change (as plotted in Supplementary Fig. 26) reflects the fact that black carbon induces increased meridional temperature gradients in the eastern Indian subcontinent, i.e., northern-increased versus southern-decreased. The enhancement of local convection and moisture source over the Indian subcontinent resulted in higher summer cloud water mixing ratios over that region ( Supplementary Fig. 27). As a consequence, localized increased precipitation reduced water vapor availability elsewhere and, as noted earlier, was coupled with a weakening downstream flow in the direction of the southern Tibetan plateau; the change in the induced regional dynamics ends up as reduced precipitation. As discussed earlier, the WRF-Chem model was turned and optimized to study the effects of increasing black carbon. With respect to the precipitation event part i.e., heavy rain days, simulations indicated that black carbon lead to increased daily precipitation ( Supplementary Fig. 29a) for most of the Indian subcontinent and, the results were consistent with the APHRODITE maximum precipitation centers (black dots in Supplementary Fig. 29a). Moreover, the WRF-Chem breakdown into the two components of convective and large-scale precipitation revealed (a) enhancement of the convective precipitation regime over the Indian subcontinent ( Supplementary Fig.   29b), while at the same time (b) a decrease in large-scale precipitation ( Supplementary Fig. 29c). The results reinforce prior lines of inquiry as to the effects of black carbon and enhanced convective instability, increased summer heavy rain and depletion of atmospheric water that is subsequently transported northward to the southern Tibetan plateau.

Glacial mass balance change over the Tibetan Plateau
The To investigate the extent of glacial mass decline, we adopted the method proposed by Brun 68  in which the signal penetrates to an unknown depth up to many meters into the snow and ice 69 . Supplementary Fig. 30 shows the spatially distributed differences in glacier

Drivers of glacial mass balance change
In order to determine the drivers of glacial change, we applied the optimization procedure from Radić and Hock 57  Applications of the empirical model also uncovered significant spatial heterogeneity in the annual glacial mass balance for 1979 to 2014 period; this is shown in Supplementary Fig. 33 67 . Besides, it is worth reiterating that summer precipitation accounts for more than 60% of the total annual precipitation that feeds the Tibetan plateau 75,76 ; this rises to 90% over the southern Tibetan Plateau.  33b) and precipitation ( Supplementary Fig. 33c). Moreover, precipitation has resulted in substantive change in glacial mass over the southern Tibetan Plateau, especially in the Himalayas. As reported by a plethora of previous analyses, summer precipitation affects the glacial mass balance in the southern Tibetan Plateau more than that in the northern parts ( Supplementary Fig. 33d).
Supplementary Fig. 34a shows the analysis of time series in glacial mass balance and its climate drivers over the Tibetan Plateau. It would seem that glacial melting by temperature is the dominant driver of glacier mass loss. Meanwhile, over time precipitations' contribution to glacial mass accumulation is unchanged. However, things are somewhat different for precipitation in the central Himalayas (southern Tibetan Plateau) as Supplementary Fig. 34b indicates a declining trend with greater variability. Moreover, the contribution from summer precipitation accounts for a larger proportion of the contribution from the annual precipitation. Finally, considering the important role of precipitation in glacial mass balance over the southern Tibetan Plateau, we investigated whether there has been a seasonal and decadal shift in precipitation, as well as its contribution to glacier mass accumulation. To account for shifts in precipitation over the southern Tibetan Plateau, we conducted a PDF (probability density function) analysis and found that the precipitation exhibited a marked seasonal pattern with significantly greater precipitation in the summer (Supplementary Fig. 35a). The contribution of summer precipitation to glacial mass accumulation ( Supplementary Fig. 33d) actually accounts for a large proportion of the contribution from the annual precipitation ( Supplementary Fig. 33c). Also, there is a decadal shift in precipitation over the southern Tibetan Plateau ( Supplementary Fig. 35b), forcing a decadal shift of summer precipitation to contribute to glacier mass accumulation ( Supplementary Fig. 35c).
Supplementary Fig. 35 The probability density functions. a seasonal precipitation during 1979-2014, b decadal summer precipitation and c their contribution to glacial mass accumulation for the 1979-2014 period.

Supplementary information for datasets used in this work
Here, we supplemented more details for the datasets used in this study, as shown in Supplementary Table 7, including the data type, spatial and temporal resolutions, and the time/period used.