Role of the Tibetan plateau glaciers in the Asian summer monsoon

The Tibetan plateau (TP) plays an important role in the Asian summer monsoon (ASM) dynamics as a heat source during the pre-monsoon and monsoon seasons. A significant contribution to the pre-monsoon TP heating comes from the sensible heat flux (SHF), which depend on the surface properties. A glaciated surface would have a different SHF compared to a non-glaciated surface. Therefore, the TP glaciers potentially can also impact the hydrological cycle in the Asian continent by impacting the ASM rainfall via its contribution to the total plateau heating. However, there is no assessment of this putative link available. Here, we attempt to qualitatively study the role of TP glaciers on ASM by analyzing the sensitivity of an atmospheric model to the absence of TP glaciers. We find that the absence of the glaciers is most felt in climatologically less snowy regions (which are mostly located at the south-central boundary of the TP during the pre-monsoon season), which leads to positive SHF anomalies. The resulting positive diabatic heating leads to rising air in the eastern TP and sinking air in the western TP. This altered circulation in turn leads to a positive SHF memory in the western TP, which persists until the end of the monsoon season. The impact of SHF anomalies on diabatic heating results in a large-scale subsidence over the ASM domain. The net result is a reduced seasonal ASM rainfall. Given the relentless warming and the vulnerability of glaciers to warming, this is another flag in the ASM variability and change that needs further attention.


Introduction
The Tibetan plateau (TP), with an area of about 2.4 × 10 6 km 2 , occupies about 0.5% of the total area of the earth. Despite this relatively small area, TP imparts a profound dynamic and thermodynamic impact on the global climate mainly due to its high elevation (average 4000 m above sea level). The TP acts as a heat source during boreal summer (Yanai et al. 1992) which drives the upper-tropospheric Tibetan anticyclone (TAC). This anticyclone along with the low-level cross-equatorial southwesterly monsoon circulation provides the necessary sheared environment to drive the precipitation northward over the Asian landmass regulating the Asian summer monsoon (ASM). As the season progresses from winter to the pre-monsoon, the heating over TP is characterized by relatively little rain and thus a fairly small moisture sink. This renders TP a significant source of sensible heat to the atmosphere (He et al. 1987) facilitating the development of TAC. The TP heating plays an important role in determining the intensity of the ASM rainfall by modulating the upper-level circulations during the spring and summer seasons. Hsu and Liu (2003) found a close relationship between interannual variability of the TP diabatic heating and the dominant East Asian summer rainfall pattern. Furthermore, using station data, radiosondes, and ISCCP data,  showed that sensible heat flux (SHF) is a major component of the TP heat source, especially before the monsoon onset (Yanai and Li 1994). These findings imply that the surface conditions over TP play an important role in determining the intensity of the ASM via modulating the SHF and contributing significantly to the total diabatic heating which in turn regulates the TAC. A detailed review of the role of TP heating in shaping the ASM climate can be found in Yanai and Wu (2006). Following the seminal study by Boos and Kuang (2010) on the mechanical impact of TP, several studies have investigated the dynamical impact of TP on ASM. Park et al. (2012) found that monsoon onset downstream of TP is mechanically driven. Chen and Bordoni (2014) built on the work of Park et al. (2012) but suggested that both effects of TP, thermal and mechanical, are important. By analyzing the model response to the extent of the northern edge of TP, Kong and Chiang (2020) suggested that the termination of the Mei-Yu is dynamically linked to the westerlies impinging on TP. The Mei-Yu terminates when the core of the westerlies migrates beyond the northern edge of TP. The role of the TP dynamical forcing was further emphasized when Chiang et al. (2020) found that without TP, ASM rainfall evolution stages disappear in their simulations. The studies by Son et al. (2019Son et al. ( , 2020 and Seok and Seo (2021) argue that orographically forced Rossby waves by TP are a primary driver of the ASM.
It is evident from previous research that both the TP effects, thermodynamic and dynamic, are important for ASM. A considerable amount of studies are available investigating the impacts of TP on ASM. However, the role of the TP glaciers has not been studied in detail yet, to the best of our knowledge. Since land cover plays an important role in determining surface heating and SHF, TP glaciers may have direct consequences on the thermodynamic impact of TP via changing the surface albedo. In a modeling study, Hu and Boos (2017) showed that surface albedo plays an equally important role as the elevation of TP itself in creating the TP "elevated heating". One may argue that TP glaciers occupy only a small area compared to the total area of the plateau (~ 2.5% of the plateau) and hence may have a negligible contribution to the total plateau heating. However, the impact of TP on global climate itself is proof that we cannot a priori rule out the potential contribution of the TP glaciers to the total plateau heating. In this study, we have investigated the impact of the TP glaciers on ASM rainfall.
To distill the role of TP glaciers in ASM rainfall, we conducted a sensitivity experiment without the TP glaciers (noGLCR run) and contrasted it with a simulation where the glaciers are present (CTRL). The TP glacier location is shown in Fig. 1a (indicated by the blue color). We used an atmospheric general circulation model (AGCM) for these simulations, which is a common practice in process studies like many of the previously mentioned studies despite the limitations of such an approach. It is noteworthy that our experimental design is highly idealized and is capable of providing only a qualitative assessment of the linear response of the TP glaciers in the current climate which is essentially the goal of this study. Nonetheless, these forced simulations are not only computationally economic but also, by being simple, can alleviate the impact of biases in various coupled processes. For example, simulations of various air-sea coupled processes are not accurate even in the most advanced climate models which may contaminate the simulation of ASM (Long et al. 2020). Moreover, since our goal is to assess the atmospheric response to glacier change, we reduce the computational costs by prescribing the sea surface temperatures. Moreover, to avoid any model response due to ocean warming, we forced the atmosphere AGCM using monthly climatological SSTs.  25° data (1998-2019) and then interpolated to model grid for visualization), c CTRL and d noGLCR. e JJAS mean rainfall (mm/day) of the noGLCR simulation minus that of the CTRL simulation. The red (blue) contour indicates a difference of 0.5 mm/day (− 0.5 mm/day). The differences above 90% significance level (based on student's t test) are shown in colors. Part of TP with elevation > 3000 m is shown by the thick black contour in panels (b-e) Especially so because Indian Ocean warming is known to impact the monsoon strength (Roxy et al. 2015;Goswami et al. 2021;Goswami 2022). The choice of prescribing a stub glacier in our simulations is also driven by the same logic of adopting a simple approach to understanding a specific process. A non-evolving stub glacier means the model does not mimic any glacier-atmosphere feedback or glacier melting. The absence of several key feedback processes of the earth system, for example, air-sea interaction, glacier-atmosphere feedback, and glacier melting mean that the model response is relatively easily relatable to the forcing fields. A key element that we did not consider in our simulations is the debris cover. It is noteworthy that simulating the effects of debris cover is a challenging task and can produce several biases (Nicholson et al. 2021) which may complicate the interpretation of the model response. The use of vegetated land units to replace the glaciers in the noGLCR experiment is a computational constraint (details are provided in Sect. 2). Although this may be a crucial caveat for a quantitative assessment of the glacier impact, for a process study like ours, this is not a primary concern. In short, our experiment design is highly idealized but well-targeted to provide a qualitative assessment of the linear response of the TP glaciers in the current climate which is essentially the goal of this study. Such an assessment is necessary to formulate more sophisticated hypotheses which can be tested using state-of-the-art earth system models.
The simulations are made with the community atmosphere model (CAM) which simulates a fairly realistic basic state of the ASM climate. It is interesting to note that the "blue" glacier patch ( Fig. 1) is similar to the orientation of the "barrier" of Boos and Kuang (2010) and ). To investigate the role of TP on the South Asian summer monsoon (SASM), Boos and Kuang (2010) and  performed sensitivity experiments from which they concluded that the TP topography is more of a "barrier" between the warm and moist tropical air and the cold and dry extra-tropical air. This further emphasizes the possible role of the glaciers as they are located exactly over the "barrier" of  where they noted a "candle heating" ("candle heating" is a term used by ) to indicate a candle-flame shaped elevated diabatic heating along the southern edge of the TP between 26°N and 30°N; also see Fig. 5b of ). Following the argument of , the absence of glaciers over TP should impact the ASM via the modulation of "candle heating" in the troposphere. However, the response of the atmosphere to the absence of this glacier mass would depend on one pivotal process: is the SHF change (due to the change in albedo) in the noGLCR run strong enough to cause a significant difference in the tropospheric heating? The thermal adaptation theory (Hoskins 1991;Yanai and Wu 2006) suggests that the SHF over TP determines the vertical profile of diabetic heating which modulates the circulation anomalies (Yanai and Wu 2006). Analyzing model simulations, we shall demonstrate (in Sect. 3) that the atmospheric response is detectable in terms of the atmospheric dynamic response to SHF changes and the related thermal forcing resulting from the glacier removal. Since we aim to understand the role of TP glaciers in the ASM climate in this study, we have restricted our analysis to the premonsoon (MAM) and monsoon (JJAS) seasons.
This paper is structured in the following fashion. In Sect. 2, we describe the observational datasets used to force the model and also for validating our model simulations. This section also includes a brief description of the model we have used in this study and a description of the numerical experiments. In Sect. 3, we present the analysis of the mechanisms of ASM response to the absence of the TP glaciers. We present our conclusions in Sect. 4.

Data and methodology
The Tropical Rainfall Measuring Mission (TRMM) Multi-Satellite Precipitation Analysis (TMPA) 3B42 Version 7 (V7) at a resolution of 0.25° × 0.25° as daily averages (1998-2019) (Huffman et al. 2010) is used to assess the model simulated mean precipitation climatology. Atmospheric thermodynamic variables are obtained from the National Centers for Environmental Prediction (NCEP) reanalysis product, version 1 (Kalnay et al. 1996).
We have performed the simulations analyzed in this study using the Community Atmosphere Model, version 6 (CAM6) (Gettelman et al. 2019). The model runs are made via the compset "FHIST" of the Community Earth System Model, version 2.1.3 (CESM 2.1.3) (Danabasoglu et al. 2020) with a horizontal resolution of 1.25° × 0.9° in longitude and latitude, respectively and 32 vertical levels from the surface to 3.6 hPa. We prescribed the glaciers as a stub glacier model (SGLC) instead of using the community ice sheet model (CISM) (Lipscomb et al. 2019). We used HadISST1 monthly climatological SSTs (Hurrell et al. 2008) to force the model. We performed one simulation with the aforementioned details which we call the control (or CTRL) simulation. To isolate the effect of the TP glaciers, we performed an experimental run which we named noGLCR. In the noGLCR simulation, we replaced the TP glaciers with natural vegetation. Replacing the glacier with a vegetated land is necessarily a subjective choice. However, out of the five allowed land units in CLM5, namely glacier, lake, urban, vegetated, and crop units, we decided to opt for vegetated land units to replace glaciers (an explanation for this choice can be found in the Supplement). It should be noted that the land unit "Lake" has the same albedo as glaciers for a frozen state while for an unfrozen state; it is a function of the solar zenith angle only. On the other hand, for a vegetated land unit, albedo is a function of various plant types and soil moisture. While replacing the TP glaciers with vegetated land units, glaciers were simply replaced and the melting of the glaciers was not considered. In these simulations (with a horizontal resolution of 1.25° × 0.9° in longitude and latitude, respectively), one may question the representation of the complex topography of TP. A highresolution model is certainly desirable, particularly to realistically mimic the glacier distribution. However, before investing in a high-resolution model simulation, preliminary investigations with relatively simpler and computationally cheaper models are needed and instructive. Simpler models are advantageous in some situations; for example in this case, for being more suitable in process understanding owing to their simplicity (Held 2005). We have used idealized approaches like replacing the glacier land units with natural vegetation land units and a stub glacier model in a climatological SST forced relatively coarse-resolution atmospheric model for the ease of interpretation of the model responses at a minimal computational cost. Further highresolution model experiments with interactive coupling between different components of land, ocean, and atmosphere will be needed to verify the robustness of our results further.
The two simulations contrasted in this study are a total of 36 years long each starting from 01 Jan 1979 onwards. The first year of each simulation is discarded as a spin-up and the remaining 35 years are analyzed.

Changes in ASM rainfall
Before doing a detailed analysis of the simulated ASM rainfall, we evaluated the simulation for a few surface variables, namely snow cover, surface temperature, and SHF.
The seasonal mean values for these variables for the pre-monsoon (MAM) and monsoon (JJAS) months are provided in the supplementary document (Figs. S1, S2, and S3). The model captures the overall spatial distribution of these surface variables reasonably well. However, the model exhibits a few biases; for example, a pre-monsoon season cold temperature bias and an overestimation of snow cover. In a study, Chen et al. (2017) also noted similar biases in CMIP5 models. The SHF biases (Fig. S3) are consistent with the surface temperature biases. Further investigation of these biases will require a detailed analysis of the forcing and feedback variables associated with snow cover (Cohen and Rind 1991) which is beyond the scope and motivation of this study. Speculatively, the impact of no-glacier-TP thermal heating on the cold and dry Siberian air , the impact of the changed meridional temperature gradient over TP on the westerly jet, etc., are among the possible mechanisms that can affect the snowfall over TP. Nonetheless, since the overall mean pattern of these surface variables is reasonably well captured by the model, we proceed with the further analysis of the ASM rainfall response in the model.
The simulated seasonal mean JJAS rainfall in the CTRL run is fairly realistic (comparing Fig. 1 c and d with the observations in Fig. 1b). However, some biases are in fact long-standing issues in the state-of-the-art climate models (Goswami and Goswami 2016). At first glance, the simulated mean JJAS precipitation with and without the TP glaciers look nearly identical. However, on a closer look, some systematic impacts of TP glaciers are evident (Fig. 1e). For a clear depiction of the seasonal mean precipitation changes, the differences are shown only for those regions where it is at least ± 0.5 mm/day (red contour indicating + 0.5 mm/day and blue indicating − 0.5 mm/ day), and within these contours, changes above 90% significance levels are shown in colors. A significant increase is noted over the central-south Indian region (centered around 18°N and 80°E) and southeastern China (centered around 27°N and 120°E) and a significant decrease in rainfall is evident over southern China (centered around 20°N and 110°E) and subtropical western North Pacific. The regions of decreasing rainfall (blues shading) are all roughly located along 20°N. This zonal orientation indicates the possibility of a large-scale organized change in the atmosphere. It is noteworthy that the decrease in rainfall over the south-central boundary of the TP, the central Indian region, and over southern China have observational analogs in rainfall trends (Yao et al. 2012;Xu et al. 2015;Yadav and Roxy 2019) indicating that in addition to the widely accepted factors forcing ASM trends such as increased aerosols (Bollasina et al. 2011), ocean warming (Roxy et al. 2015;Goswami et al. 2021;Goswami 2022), etc., the retreating TP glaciers can also be a contributing factor; an issue that warrants further research. However, our motive in this study is to understand the role of the TP glaciers in the context of ASM dynamics. Considering the large, zonally oriented pattern of decrease in rainfall in the absence of the TP glaciers (in the noGLCR run) compared to a climate with the TP glaciers (in the CTRL run), we conjecture that a relatively large-scale impact of the TP glaciers is possible.

Changes in the TP thermal forcing
The prime candidate for explaining contrasts in large-scale ASM dynamics in the CTRL and noGLCR runs is the change in the thermal forcing caused by changes in surface properties, predominantly the albedo resulting from the removal of glaciers (Takeuchi et al. 2002). However, it should be noted that the surface albedo is dependent on the total snow cover/ thickness and not just the glaciers. The presence of snow cover and its melting affect the surface temperature because of its high reflectivity (higher albedo), higher thermal emissivity (more outgoing infrared radiation), and low thermal conductivity (less sensible heat flux from the ground to the atmosphere) (Cohen and Rind 1991). When we simulated the climate without the TP glaciers, the absence of glaciers was most felt in regions that were not snow-covered otherwise (that is, regions that are climatologically not snowy or relatively less snowy). The northwestern TP is a climatologically heavy snowfall region (Wang et al. 2018). Considering the glaciers and the climatologically snowy regions (Fig. S1), it is evident that the south-central boundary region of TP experienced the most contrasting surface conditions in the noGLCR experiment compared to the CTRL run during the premonsoon season. Since the contribution of the SHF and its impact are at their peak during the pre-monsoon season (MAM), we investigated the MAM difference in albedo, surface temperature, and SHF between the noGLCR and CTRL simulations. Removing the TP glaciers results in a decreased albedo (within the magenta-colored contours which indicate the location of the glaciers) (Fig. 2a) which means less insolation is reflected. Figure 2a shows two positive albedo anomalies located at about (86.5°E, 32.5°N) and (98°E, 32°N) (these locations are marked by X1 and X2, respectively, in Fig. 2a). However, these responses are not significant because of the high climatological mean and variability of albedo in these locations (Fig. S4). Moreover, as we did not explicitly alter the surface conditions in our sensitivity experiment anywhere except over the glacier regions, we start our investigation of the albedo anomalies noted within the magenta-colored contours in Fig. 2a. The negative albedo anomalies along the south-central boundary of TP in Fig. 2a lead to warmer Fig. 2 MAM mean a surface albedo (the green contours indicate the differences above 90% significance level (based on student's t test)), b surface temperature (no value is found to be significant above 90% confidence level) and c SHF in the noGLCR simulation compared to the CTRL simulation. d JJAS mean SHF in the noGLCR simulation compared to the CTRL simulation. In c, d, only differences above 90% significance level are shown. Part of TP with elevation > 3000 m is shown by the thick black contour. The purple contour indicates the TP glacier occupied region. In a, X1 and X2 mark the locations (86.5°E, 32.5°N) and (98°E, 32°N), respectively surface temperatures (Fig. 2b) in that region. This results in an increase in SHF (Fig. 2c). One intriguing point here is that the strong increase in SHF is also seen during the monsoon season in the noGLCR climate over the western TP (Fig. 2d). This is consistent with the fact that during the JJAS season, the western TP, which otherwise would have been glaciated, is a vegetated surface in the noGLCR run. And hence it has a lower albedo, a warmer surface, and higher SHF. In addition, these positive SHF anomalies in the western TP are a case of retaining the SHF memory into the monsoon season because of the circulation response to the locally rising air (Wang et al. 2014). We shall discuss this increase in SHF in the western TP persisting into JJAS in Sect. 3.3. As for MAM, we re-emphasize that the western TP is a climatologically snowy region that remains snowy in the noGLCR run as well (figure not shown). Hence, the impact of glacier removal is not very strong (Fig. 2c) in that region.
An increase in the SHF in the absence of the TP glaciers over the south-central boundary of TP during the pre-monsoon season is seen in Fig. 2c. It occurs during the seasonal transition of TP from a heat sink to a source (He et al. 1987;Yanai et al. 1992). The SHF over the plateau determines the vertical profile of diabatic heating which modulates the circulation (Yanai and Wu 2006). Furthermore, the contribution of SHF becomes prominent during the pre-monsoon season (He et al. 1987; also see Fig. S5). To investigate the possible ASM circulation response to these SHF differences, we analyzed the simulated diabatic heating, the associated meridional circulation, and the resultant tropospheric temperatures. These latitude-height cross-sections in Fig. 3a-d are computed over the longitude band 80°E-100°E, where maximum changes in SHF are evident in Fig. 2c. In our analysis presented in Sect. 3.3, we shall demonstrate that the positive SHF anomalies over the western TP (Fig. 2d) are driven by the positive SHF anomalies over the south-central TP (Fig. 2c). We chose the longitude band, 80°E-100°E, in Fig. 3a-d, to depict the impact of the positive SHF anomalies over the south-central TP without considering SHF anomalies over the western TP. During the pre-monsoon season (Fig. 3a), prominent anomalous heating is evident in the noGLCR run relative to the CTRL simulation over the southern boundary of TP; exactly over where the glaciers are located in the CTRL run and are absent in the noGLCR run. Interestingly, this anomalous heating looks somewhat similar to the "candle heating" of . The anomalous heating causes an updraft at 30°N and subsidence to the south and north of it as shown in Fig. 3c. This is consistent with the results of Sato and Kimura (2007), who performed a series of sensitivity experiments to elicit the mechanical and thermal impacts of the TP plateau. They showed that the TP thermal forcing causes subsidence over the northern Indian region during the pre-monsoon season which weakens as it moves poleward with the season and the onset of the Indian summer monsoon occurs. The associated changes in the tropospheric temperature field and the subsequent subsidence around 50°N latitude cause warming below 400 hPa, leading to a marginal increase in the low-level clouds (figure not shown) (Bony et al. 2004). In Fig. 3 a and c, during pre-monsoon (MAM) when the monsoonal latent heating is not strong, the overturning circulation can be seen being driven by the strong heating on the southern slope of the TP. However, for the monsoon season, from the heating and meridional circulation shown in Fig. 3 b and d, respectively, it is not immediately obvious that the JJAS rainfall changes are directly attributable to the absence of TP glaciers. Nonetheless, a decreased diabatic heating and increased subsidence in the latitude band 10°N-20°N is consistent with the decreased rainfall (Fig. 1e) over the central Indian region. Possibly, the pre-monsoon anomalous heating and the associated circulation and temperature responses which are attributable to the absence of TP glaciers persist into the monsoon season and drive a seasonal change during JJAS. The mechanisms for this persistence are discussed below.

Persistence of spring SHF anomaly into summer
To find a causal link between the pre-monsoon anomalous heating over TP and the monsoon circulation over the ASM region (Wang et al. 2014), we examined the longitudeheight cross-sections of the zonal circulation (averaged over 30°N-40°N) during the premonsoon season (Fig. 3e) and diabatic heating during the monsoon season (Fig. 3f). The motive behind depicting the pre-monsoon zonal circulation over TP is to demonstrate that when the TP heating interacts with the westerlies, the tropospheric warm air rises on the eastern side and descends gradually on the western side (Rodwell and Hoskins 1996;Wang et al. 2008). Our Fig. 3e is reminiscent of Fig. 2c of Wang et al. (2008; hereafter W08) whose altered albedo experiments can be logically considered parallel to our CTRL (increased albedo, or, "C" runs in W08) and noGLCR (reduced albedo, or, "W" runs in W08) runs. Furthermore, based on the same theoretical premise, Wang et al. (2014) showed that a positive SHF anomaly over the eastern TP during the pre-monsoon season drives a positive SHF anomaly over the western TP during the monsoon season thereby retaining the SHF memory into the following season. In our simulations, the positive SHF anomalies over the south-central TP during MAM (Fig. 2c) drive an anomalous pre-monsoon zonal circulation (Fig. 3e) leading to positive SHF anomalies over the western TP during JJAS (Fig. 2d) that results in an anomalous positive diabatic heating over the western TP region during the monsoon season (Fig. 3f). What is the response of ASM circulations to this heating? Wu et al. (2007) performed idealized experiments with an aquaplanet setup but with an idealized topography prescribing elevated land in the TP region. They found that sensible heating on the sloping lateral surfaces of the elevated land is a major driver of the ASM circulations. Using a series of sensitivity experiments, Wu et al. (2012) demonstrated that surface sensible heating in the Himalayan slope of TP is essential for , and e, f over the latitude band 30°N-40°N. Green contours indicate differences above 90% significance level (based on student's t test). Black shading indicates the TP region topography averaged, a-d over the longitude band 80°E-100°E, and e, f over the latitude band 30°N-40°N ASM rainfall.  also reported similar findings emphasizing the importance of the TP thermal heating in the southern fringe of TP. To mimic the TP warming due to sensible and latent heating, W08 imposed a heat source over TP centered at 90°E, 33°N, with maximum heating in the boundary layer, decreasing exponentially with height. W08 studied the response of the monsoon circulations using an atmospheric general circulation model and found two Rossby wave responses. One along the upper-level westerly jet stream to enhance the anticyclonic circulation to the east of Japan and another wave train along the low-level southwesterly monsoon into the South China Sea. Since we found a considerable resemblance between the heating imposed in W08's experiment and the heating resulting from the absence of TP glaciers in our study, we investigated whether similar circulation responses are seen as suggested by the experiments of W08.

Circulation response to spring SHF anomaly
The circulation responses at 200 hPa and 850 hPa in the absence of TP glaciers are depicted in Fig. 4. Two anticyclonic circulation anomalies are seen at 200 hPa (Fig. 4a); one over the eastern side of TP (just east of the 3000 m height contour) centered at 111°E, 29°N, and the other located further eastward centered at 145°E, 25°N. This response is different from what W08 had observed. An anticyclonic circulation response centered over TP was noted by W08. Moreover, the lower-level anomalous southwesterlies and the upper-level anomalous anticyclonic circulation (marked by the letter A in Fig. 4a) indicate a baroclinic structure suggesting that the anticyclonic response A is driven by the latent heating from the anomalous rainfall seen in Fig. 1e. This difference in the results of W08 and our experiment can be explained as follows: in W08, the TP heating is stronger because of the 3 times albedo difference between their "W" and "C" runs and also because the albedo difference is forced over the entire TP. On the other hand, in our experiment, the heating is generated only by the absence of the "blue" strip of glaciers (Fig. 1a). Hence, the non-linear effect and the eddies due to the larger latent heating obscure the impact of the TP heating anomalies in the noGLCR run compared to the CTRL in Fig. 4a. The 850 hPa circulation response is depicted in Fig. 4b. Consistent with the results of W08, we also note a wave train response at 850 hPa over the SASM domain south of TP. However, the circulation response east of Japan at 850 hPa is cyclonic, unlike W08, which is related to the anomalous increase in rainfall in that region. From the upper and lower-level circulation responses seen in our experiment in the absence of TP glaciers, we posit that the responses that we note are fundamentally similar to those of W08 except for some details which are attributable to the differences in the experimental set ups. The differences are due to the different strengths of heating imposed in the two studies. While W08 analyzed the response of ASM to TP warming as a whole, our results elicit the role of the glaciers in the total TP thermal forcing in terms of ASM circulation responses.
When we analyzed the divergence (Fig. 4c, d) associated with the circulation response, we get a better insight into the precipitation anomalies seen in Fig. 1e. Compared to the CTRL run, there is a large-scale convergence at the upper levels (Fig. 4c) and a low-level divergence (Fig. 4d), indicating large-scale subsidence over the ASM region at about the 20°N latitude. This is consistent with the reduced vertical velocity at 500 hPa (Fig. S6). However, the reduction in vertical velocities at 500 hPa (and the rainfall response in Fig. 1e) in the absence of TP glaciers are not as smooth as the large-scale circulation response. Rainfall at a given location depends on several factors involving local and remote drivers. The divergence patterns (Fig. 4c, d), which are suggestive of large-scale subsidence over the Asian monsoon region, does not necessarily guarantee a smooth and homogeneous reduction in rainfall. It only indicates more inhibition to convection. Nonetheless, a careful visual inspection reveals that the negative rainfall anomalies (Fig. 1e) are almost collocated with the centers of strongest subsidence (Fig. 4c, d). While the regional inhomogeneity in the rainfall response may depend on several factors, both local and remote, it is clear from Fig. 4c, d that the absence of TP glaciers would cause an overall inhibition to ASM rainfall, especially along ~ 20°N.

Discussion and conclusion
The TP glaciers are the "water tower" of Asia. However, glaciers do not just provide water, they serve as a highly reflective land surface as well. Therefore, their retreat or disappearance can potentially impact the surface properties significantly. The resulting albedo changes can potentially manifest in atmospheric heating and hence the ASM rainfall. We investigated the role of TP glaciers in ASM by analyzing the CAM model sensitivity to TP glaciers by replacing them with natural vegetation. Although the TP glaciers may not disappear completely soon, the recent rapid retreat of the glaciers cannot be ignored (Yao et al. 2012). The impact of TP glaciers on ASM, some of whose (the ones on the southern and eastern TP) accumulation depends on the ASM precipitation itself (Fujita 2008), determines the net feedback of glacier changes. This underscores the need to understand the role of the TP glaciers in the ASM dynamics. In our experiment, by replacing the TP glaciers with natural vegetation, we mainly altered the surface albedo which in turn affected the SHF. Because SHF contributes significantly to the total diabatic heating when the TP transitions from a heat sink to a source during the pre-monsoon to monsoon season, the atmosphere responds to these SHF changes. Response of the atmosphere to changed surface albedo has also been studied before (for ex., Wang et al. 2008). However, in this study, we demonstrated that even though the glaciers only occupy the west-south-east boundary regions which may seem relatively small compared to the whole TP domain, changes in surface properties in that glacier strip still have a substantial role in driving atmospheric responses. This dynamic link may have crucial implications for the future of ASM. To wit, in the absence of TP glaciers, we observed a decrease in the seasonal mean ASM rainfall (mostly along the 20°N latitude) (Fig. 1e). Based on the investigation of the model-simulated dynamic and thermodynamic responses, we summarize our results as follows: 1. The absence of glaciers is felt most along the south-central boundary of the TP during the spring season because of the decrease in albedo and the resulting increase in SHF. 2. The increased SHF results in increased diabatic heating along the south-central boundary of the TP. 3. The spring season heating anomalies over TP drive a zonal circulation anomaly with rising air in the east and subsiding air in the west which persists in the increase in SHF into summer. 4. This leads to a positive diabatic heating anomaly over the western TP which in turn drives an atmospheric circulation that is zonally oriented (along 20°N) and results in large-scale subsidence over the ASM domain. The net result is an inhibition of convection and a decreased ASM rainfall.
Despite the limitations (briefly described in supplement Sect. 4) that our experimental design did not consider glacier melt, debris, and the potential feedback but instead abruptly replaced the glacier with natural vegetation, which removes some feedback processes involving glaciers, land, ocean, and atmosphere, we posit that the results are still insightful. Since these key processes are missing both in the control and no-glacier simulation, the results are indicative of the model response to the absence of glaciers. The most intriguing result is that TP glaciers play a demonstrable role in the thermal heating over the entire TP. Arguably, the inclusion of the additional feedbacks may alter the model response. It is noteworthy that the differences between W08 and this study emphasize the fact that the response of the atmosphere and hence the ASM climate are sensitive to the location of the warming within the TP with the caveat that the experiments in the two studies are not identical. Furthermore, the TP warming may not occur homogeneously in the future (Shi et al. 2019). A more detailed analysis will thus be needed to affirm the robustness of these results. Another major implication is that this study emphasizes the need for a better representation of the glaciers in climate models. In a recent study, analyzing the summertime warming trend, Shi et al. (2019) reported that the CMIP5 models overestimate the southern TP warming rate in historical simulations. Based on our study, the possibility of this overestimated warming could be a contributor to CMIP5-simulated systematic dry bias over the ASM domain (Narapusetty et al. 2016).
A remaining question is the response of the albedo in the TP domain over the glacierfree regions. Figure 2a shows two positive anomalies centered at about (86.5°E, 32.5°N) and (98°E, 32°N). Since we did not alter the surface conditions anywhere else in our sensitivity experiment other than over the TP glacier regions, these albedo responses are intriguing. A detailed investigation of these responses is needed to bring forth the processes that generate these albedo anomalies. Nonetheless, the change in snow cover in the noGLCR run compared to the CTRL suggests that these albedo responses are strongly related to snow cover changes (Fig. S8). Snow cover over the complex terrain of TP is a result of interactions of dominant atmospheric circulations transporting moisture over to TP (You et al. 2020). The fact that these two positive albedo responses in the noGLCR run are located to the southeast of the two major elevated regions in central TP, the Chang Tang plateau and the Tanggula Mountains (Fig. S7), suggests a possible role of the TP terrain in causing these responses. Interaction of the upper level East Asian jet with the orography can cause snow plumes off the highest mountain peaks over TP (Moore 2004) and redistribute snow cover (Pu et al. 2007). At this point, it is worth mentioning that a marginally stronger and slightly equator-ward shifted westerly jet is indeed seen in the noGLCR run (figure not shown) which must be a consequence of the altered meridional temperature gradient (figure not shown). Meridional displacement of the westerly jet was reported to be strongly associated with the Silk Road Pattern (SRP) teleconnection, which is a stationary Rossby wave train along the Asian jet in the upper troposphere (Enomoto et al. 2003). Downstream Rossby waves can be generated also by the variability of the ASM precipitation (via the resultant diabatic heating variability) (Lin et al. 2016;Stephan et al. 2019;Son et al. 2021). In such a case, the absence of the TP glaciers can kick off a chain of feedback by reducing the ASM rainfall. Robust evidence would be needed to comment conclusively, for which further detailed analysis is required.
The fact that TP is experiencing rapid warming and this warming further increases with altitude (Qin et al. 2009), makes the TP glaciers extremely vulnerable to a changing climate. The fact that the glacier-occupied regions in the southeast and eastern boundaries of TP are not 100% snow-covered (Fig. S1) underscores the importance of the glaciers in determining the local surface albedo. Hence glaciers and their shrinkage can potentially impact the SHF which in turn should affect the TP total heating and eventually the ASM rainfall itself. This is likely to have crucial implications because while over the western TP, the glaciers accumulate snow mostly in winter through westerly atmospheric circulation, the glaciers in the eastern TP receive most accumulation during summer from the Indian monsoon precipitation (Fujita 2008). Therefore, a change in the ASM rainfall arising from the TP glacier shrinkage, if negative, may form a positive feedback loop. A summer snowfall decrease may accelerate the glacier loss further via a decrease in accumulation due to a decreasing ASM rainfall (Duan et al. 2004). A detailed study of the role of the TP glaciers in the ASM climate is therefore extremely important. Our results are very qualitative. To quantify our findings under realistic conditions, further research is warranted with coupled models with realistic representations of debris, transient evolution of the glaciers, etc.