Globally elevated chemical weathering rates beneath glaciers

Physical erosion and chemical weathering rates beneath glaciers are expected to increase in a warming climate with enhanced melting but are poorly constrained. We present a global dataset of cations in meltwaters of 77 glaciers, including new data from 19 Asian glaciers. Our study shows that contemporary cation denudation rates (CDRs) beneath glaciers (2174 ± 977 Σ*meq+ m−2 year−1) are ~3 times higher than two decades ago, up to 10 times higher than ice sheet catchments (~150-2000 Σ*meq+ m−2 year−1), up to 50 times higher than whole ice sheet means (~30-45 Σ*meq+ m−2 year−1) and ~4 times higher than major non-glacial riverine means (~500 Σ*meq+ m−2 year−1). Glacial CDRs are positively correlated with air temperature, suggesting glacial chemical weathering yields are likely to increase in future. Our findings highlight that chemical weathering beneath glaciers is more intense than many other terrestrial systems and may become increasingly important for regional biogeochemical cycles. Global glacial chemical denudation is one of the largest contributors to global elemental cycles and, amplified by climate warming, will significantly impact nutrient loads in downstream ecosystems.

Here we present concentrations of major cations in meltwaters of the 19 glaciers within Asia. The new data is combined with a global analysis of CDRs in glacial environments to help elucidate the role of future climate-driven glacier retreat and meltwater export in regional chemical weathering elemental cycles. We propose that the subglacial chemical weathering rates are directly correlated to changes in climate and glacial runoff, and subglacial chemical weathering will become increasingly important for regional elemental cycles and weathering-climate feedbacks in response to accelerated global glacier (and ice sheet) mass loss in a warming environment [47][48][49][50][51] .

Results and discussion
Glacial cation concentrations. We developed a global database of major cation concentrations from 77 glaciers (63 mountain/ valley glaciers and 14 ice sheet outlet glaciers) containing 5465 samples generated in this study and published literature ( Fig. 1A and Supplementary Tables 1-3). For our 19 Asian glaciers, the dominant cation is Ca 2+ for 14 of the 19 glaciers, followed by Mg 2+ , Na + , and K + (Fig. 1B). This is consistent with published cation data from glaciers in Asia ( Supplementary  Fig. 1). However, at SG2, SG3, and SG4 in the Qilian and at MKG and YZG in the Kunlun, meltwater compositions are dominated by Mg 2+ (Fig. 1B), suggesting the contribution of magnesium-rich minerals [52][53][54] . Globally, Ca 2+ is also the dominant cation in most glacial basins ( Supplementary Fig. 1b). However, the dominant cation is Mg 2+ in some basins within Asia, and Na + and K + in Scandinavia and Svalbard. The basinmean cation concentrations range between one and two orders of magnitude, with arithmetic mean of 159 ± 341 μeq L −1 for Na + , 42.7 ± 46.0 μeq L −1 for K + , 393 ± 655 μeq L −1 for Mg 2+ , and 581 ± 225 μeq L −1 for Ca 2+ for the 19 Asian glaciers ( Fig. 2a and Supplementary Table 4).
The range of cation concentrations for the 19 Asian glaciers (where data from 4 glaciers extended existing records) overlaps the majority of previously published data from 28 glaciers in Asia and from other glaciers worldwide ( Supplementary Fig. 2). However, mean cation concentrations for the 19 Asian glaciers are generally higher than previously published cation means for Asian glaciers and worldwide glaciers. For the 77 glaciers in the current global data set mean cation concentrations vary over two to three orders of magnitude (3.38 to 1323 μeq L −1 for Na + , 2.31 to 177 μeq L −1 for K + , 6.50 to 2605 μeq L −1 for Mg 2+ , and 11.5 to 1701 μeq L −1 for Ca 2+ ; Fig. 2b). Mean concentrations from glaciers (93.3 ± 177 μeq L −1 for Na + , 34.6 ± 29.0 μeq L −1 for K + , 191 ± 321 μeq L −1 for Mg 2+ , and 472 ± 366 μeq L −1 for Ca 2+ ; Supplementary Table 5) are comparable to those from global non-glacial rivers (265 μeq L −1 for Na + , 59.0 μeq L −1 for K + , 342 μeq L −1 for Mg 2+ , and 900 μeq L −1 for Ca 2+ ) after accounting for the uncertainty 55,56 .
Diurnal variations in cation concentration. For the 19 Asian glaciers, hourly cation concentrations differ between one to three orders of magnitude (Fig. 3a, d, g, j) and display strong diurnal trends, with higher values in the morning hours (22:00-12:00 h the following day) and lower values in the afternoon and evening hours (13:00-21:00 h; Fig. 3b, e, h, k). Compared with seasonal change of glacial cations 30,45,46,[57][58][59] , the diurnal changes reflect the interaction time between water and rock/sediment coupled with meltwater dilution 26,60,61 . Field observation showed that during the morning hours lower air temperature reduces the supraglacial melt, leading to less meltwater, lower proglacial runoff, and higher cation concentrations 26,[60][61][62] . Conversely, cation concentrations are diluted by increasing meltwater discharge and reduced transit times in the subglacial hydrological system later in the day (Fig. 3b, e, h, k). Our data are consistent with the diurnal changes for glacial runoff cations and/or electrical conductivity from Haut Glacier d'Arolla in Central Europe 61 , and from Qiyi Glacier and Dongkemadi Glacier in Asia 26,60,62 . To identify the "most representative" daily sampling time, the ratios of hourly to daily mean cation concentrations were calculated (Fig. 3c, f, 2+ during the afternoon and evening hours with low cation concentration and high runoff. The daily amplitude of these ratios provides a template that sampling at 13:00 h or~21:00 h can capture the "most representative" meltwater ( Fig. 3c, f, i, l). This can be used to help estimate the daily mean concentrations and therefore the fluxes of glacial cations in future work on similar glacial systems.
Glacial chemical weathering and cation sources. Glacier meltwater chemistry is related to the principle weathering reactions, rock mineralogy, and drainage system configuration 19,21,41,57,63 . The Gibbs plots and elemental stoichiometry were used to identify the sources of cations 13,63,64 . Cation weathering products from glacial basins and regions confirm the dominance of rock weathering as the major source (Fig. 4a, b). Notably samples from MKG and YZG plot away from the other glaciers, which may be related to differences in the underlying geology. Na-normalized ratios for Ca 2+ and Mg 2+ are a useful tool for exploring the relative importance of different weathering processes 26,33,65 . Clusters from the 19 Asian glaciers are widely distributed between silicate and carbonate end-members, but MKG and YZG are close to the silicate end-member (Fig. 4c). This suggests that glacial chemical weathering in the Kunlun has a higher relative proportion of silicate to carbonate weathering than in other mountain ranges within Asia. In detail, carbonate weathering dominates the weathering regime for glaciers in Asia, Low Latitudes, Arctic Canada and Europe; however, in Svalbard, Iceland, Scandinavia, Greenland, and Antarctic and Subantarctic silicate weathering appears to be comparatively more important (Fig. 4d and Supplementary Table 6).
The chemical composition of meltwaters has been shown to be typically controlled by the presence of carbonate weathering in most glacial environments 6,13,26,27,32,34,41,45,46,60,66 . Carbonate minerals (e.g., calcite, dolomite, and siderite) are preferentially weathered from glacial flour even if bulk bedrock is dominated by silicate minerals (e.g., potash feldspar, plagioclase, and illite) 13 . Previous studies also indicate that Ca 2+ and Mg 2+ are primarily derived from carbonate weathering and Na + and K + from silicate weathering in glacial environments 19,27,[67][68][69] . Our data set is consistent with this, reflected by higher contributions of Ca 2+ and Mg 2+ to total cation load from the 19 Asian glaciers (71-97%) and from previously published glaciers (51-98%; except for Scott Turnerbreen with 49%; Supplementary Tables 3-5) 31 . In detail, dolomite may be a source of Ca 2+ and Mg 2+ at SG2, SG3, and SG4 which is supported by their similar concentrations (Supplementary Table 3 and Table 7). However, at MKG and YZG silicate and evaporite may be a key source of Mg 2+ and Na + which fits with higher concentrations of Mg 2+ and Na + than Ca 2+ in glacial meltwaters (Fig. 4c, d).
Export of cations from glaciers. Estimates of cation fluxes should ideally be based on the discharge-weighted or midsummer normalized cation data (to correct for variable flow and associated concentration effects; Figs. 2 and 3) 26,45,46,57,61 . However, our global compilation of glacial cations (Supplementary Table 3) shows that discharge data are sparse and available for only a few glaciers. Future studies need to address this NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-28032-1 ARTICLE omission which will seriously hamper attempts to compare the regional weathering rates and generate accurate flux (and CDRs) estimates. Here we proceed to use our mean cation concentration data in combination with modeled meltwater flux data to produce a first-order estimate of regional cation fluxes from glaciers in Asia and other glacial regions and then generate a global estimate (Table 1). This approach is tested against estimates based on the methods of regional discharge-weighted cation concentration and regional discharge extrapolation using all data in nine non-ice sheet glacial regions (apart from Greenland Periphery, and Antarctic and Subantarctic) and the mid-summer data in six non-ice sheet glacial regions (apart from Alaska, Arctic Canada, and Svalbard), respectively, in the current global data set ( Fig. 1A and Supplementary Table 8).
Within Asia, the highest cation flux of 337 Gg year −1 is found in the Kunlun, which has the largest glacier cover (11,524 km 2 ) and mean cation concentrations (apart from K + ), despite the lowest specific discharge (0.43 m year −1 ) of mountain ranges ( Fig. 5a and Table 1). The other two mountain ranges with large glacier cover, the Tianshan (7180 km 2 ) and Himalayas (6821 km 2 ), also have high cation fluxes of 156 Gg year −1 and 271 Gg year −1 , respectively, compared to all other Asian glaciers. The smaller glacier-covered mountain ranges of Hengduan (1395 km 2 ) and Altai (179 km 2 ) have the lowest cation fluxes of 14.9 Gg year −1 and 11.8 Gg year −1 , respectively (Fig. 5a). For the Asian glacial region (121,694 km 2 ) we estimate a regional cation flux of 6305 ± 5843 Gg year −1 from Asian glaciers ( Fig. 5b and Table 1). The glacial region with the highest glacier cover, Arctic Canada (145,767 km 2 ), is also the region with the highest cation flux (7796 ± 321 Gg year −1 ). The lowest cation flux of 37.0 ± 36.4 Gg year −1 is found in Scandinavia, one of glacial regions with least glacier cover (2851 km 2 ; Fig. 5b and Table 1).
For eight mountain ranges in Asia, glacial cation fluxes are linearly related to glacial area (R 2 = 0.90, p < 0.01; Supplementary  Fig. 3a). This indicates a relationship between glacial area and runoff for mountain ranges (R 2 = 0.65, p < 0.02; Supplementary  Fig. 3b). There are some exceptions to these trends. The Kunlun has a higher cation flux (337 Gg year −1 ) than the Himalayas (271 Gg year −1 ) even though glacial runoff is smaller (by around 5 km 3 year −1 ; Table 1). For ten glacial regions (including Greenland Periphery), cation fluxes are also closely related to glacial area (R 2 = 0.74, p < 0.01; Supplementary Fig. 3c), suggesting a close relationship between glacial area and melt (R 2 = 0.95, p < 0.01; Supplementary Fig. 3d). Notably, cation export from Asia only accounts for 25% of that from nine non-ice sheet glacial regions (see above) despite its relatively large percentage (32%) of total runoff (Table 1). Conversely, glacial runoff from Arctic Canada accounts for 19% of total runoff yet appears to account  Table 4); b 77 glaciers in eleven glacial regions (Supplementary Tables 3 and 5). The boxes and horizon lines above indicate the mean cation concentration and standard deviations, respectively.
for 31% of cation export. This indicates that using a single glacier or even a few glaciers to represent a whole glacial region may introduce large uncertainties into the calculations of cation export owing to large spatial variations of cation yield and/or weathering processes.
Given that glacial runoff from these nine non-ice sheet glacial regions (1130 km 3 year −1 ; Table 1) accounts for 79.0% of total glacial runoff globally (excluding ice sheets) 70 , we can crudely estimate worldwide glacial cation export as 32,035 ± 14,446 Gg year −1 based on the discharge extrapolation (Table 1). This equals that calculated by the regional discharge-weighted mean cation concentrations from all glaciers in these nine non-ice sheet glacial regions multiplied by total glacial runoff (Supplementary Table 8). Based on mean cation concentrations from nine outlet glaciers in the Greenland Periphery 22,54,[71][72][73] (Supplementary  Table 5) multiplied by total runoff (542-652 km 3 year −1 ) 9,74,75 , total cation export from the Greenland Ice Sheet (GIS) was estimated as ranging from 2158 to 2595 Gg year −1 (Table 1), consistent with another recent estimate of cation release from the GIS (2386 Gg year −1 ) 76 . Cation release from the Antarctic Ice Sheet (AIS) was estimated as 1209 to 5584 Gg year −1 based on mean cation concentrations from five outlet glaciers in the Antarctic and Subantarctic multiplied by total runoff (55-254 km 3 year −1 ) 77,78 (Table 1 and Supplementary Table 5).
Cation export from glaciers wordwide is~3% of that from global non-glacial rivers (1,172,786 Gg year −1 ; Table 1). However, we estimate it is currently higher (3.9-9.5 times) than the combined GIS and AIS, even though total runoff 70 is only 1.6 to 2.4 times higher ( Fig. 5c and Table 1).
Our estimates contain substantial uncertainties, especially in the absence of catchment discharge-weighted mean cation concentration data from glaciers worldwide.  Table 1).

Fig. 3
Hourly variations in cation concentration as well as the anomalies and ratios for the 19 Asian glaciers generated in this study over a three-day period. a, d, g, j Hourly concentrations; b, e, h, k Anomalies of hourly concentration (i.e., the difference between hourly and mean concentrations divided by mean concentrations); c, f, i, l Ratios of hourly to mean concentrations. Note that the diurnal changes were indicated by gray shadows (M and AE indicate the morning and afternoon/evening, respectively), the mean anomalies and ratios were indicated by green solid lines, and the pink dots indicate the time at 13:00 and 21:00 when hourly concentrations equal to the daily mean.  Fig. 4) suggesting that the proglacial weathering will have relatively minor effects on cation exports in most cases. (iv) If only the mid-summer data (n = 593) is used from 24 non-ice sheet glaciers worldwide, the estimated glacier cation flux is~64% of that using all data (n = 4641) from all non-ice sheet glaciers, which is lower but does not significantly change the conclusions of this study (Fig. 1A  and Supplementary Tables 8, 9).
(v) Geochemical data is currently not available for the large glaciers in Patagonia and Southern Andes (Fig. 1A). Although total runoff is comparatively small compared to global glacier runoff, glaciers in these regions are among the most rapidly retreating on Earth, and are likely to be regionally important in cation denudation.
Despite these reservations, our findings provide evidence that glaciers play a key role in catchment geochemical weathering and associated elemental mobilization.
Glacial chemical weathering rates. Chemical weathering rates for glaciers within Asia appear to be controlled by several variables. Glacial weathering in the Altai is impacted by a very high discharge (3.63 m year −1 ) 82,83 , and the high observed Ca 2+ concentration (750 μeq L −1 ) is the main component of the very high CDR (3425 Σ*meq + m −2 year −1 ) ( Fig. 5d and Table 1). The lowest CDR of 486 ± 105 Σ*meq + m −2 year −1 is found in the Hengduan, where the specific discharge (0.65 m year −1 ) and cation concentration are relatively low. Higher CDR from the Himalayas (2000 ± 1586 Σ*meq + m −2 year −1 ; Fig. 5d) is also consistent with a positive relationship between CDRs and specific discharge (Table 1). Lower CDRs and specific discharge along with slow glacier mass loss for mountain ranges on the Tibetan Plateau 84,85 is further evidence of a link between meltwater generation and chemical weathering rates. On a regional scale, glacial CDR for Asia (2655 ± 2565 Σ*meq + m −2 year −1 ) is similar to Arctic Canada, Svalbard and Jan Mayen, Alaska,  Table 6). Note that data within the blue dashed circles denote Meikuang Glacier (MKG) and Yuzhufeng Glacier (YZG) in the Kunlun (Fig. 1B), and cations are not corrected for the sea-salt and aerosol contributions.     Fig. 5e and Table 1).
Our data indicate that regional variability in glacial chemical weathering rates is related to the hydrological regime. We find a positive linear relationship between glacial CDRs and specific discharge on both mountain range (for eight mountain ranges in Asia; R 2 = 0.77, p < 0.01) and regional (for ten glacial regions worldwide; R 2 = 0.46, p < 0.04) scales (Fig. 6a, b). The highest CDR of 9683 Σ*meq + m −2 year −1 is estimated for glaciers in the low latitude region, where the specific discharge is very high (6.39 m year −1 ). The lowest regional glacial CDR of 285 ± 233 Σ*meq + m −2 year −1 we estimated occurs in Greenland Periphery (with a specific discharge of 1.66 m year −1 ; Fig. 5e and Table 1). The disproportionately high denudation rates from low latitude glaciers reflect high rates of mass loss 50,70,86 as well as potentially higher rates of physical erosion compared to glaciers at high latitudes 14,87 , creating an abundance of reactive mineral surfaces.
CDRs from glacial basins worldwide also show a linear relationship to specific discharge on a basin scale (for >29 glacial basins globally; R 2 = 0.25, p < 0.01; Fig. 6c and Supplementary  Table 10). Glacial CDRs are positively related to mean annual temperature (0.18 < R 2 < 0.51, p < 0.03) and mean annual precipitation (0.26 < R 2 < 0.42, p < 0.05) on both basin and regional scales (Fig. 6e, f, h, i). This indicates that climate warming driven meltwater runoff increase 50,51 is a key control of glacier chemical denudation. The inverse CDR response to increasing latitude on both basin and regional scales (0.40 < R 2 < 0.69; p < 0.01; Fig. 6k, l and Supplementary Table 10) is further evidence of the link between climatic warming, meltwater generation (specific discharge), and chemical weathering rates. This suggests that cation export will increase as glacier and ice sheet melt accelerates, while it will decrease as glaciers and/or ice sheets are removed from the landscape (or reaches a critical threshold). Geographically widespread studies are needed to further constrain estimates  Table 1). Note that L and H in the bottom plots (c, and f) denote the high-and low-end estimates, respectively, and the horizontal lines above the bars indicate the standard deviations. and help further address these relationships. Over the next few decades, non-ice sheet glaciers are likely to become increasingly important in chemical weathering budgets and related regional elemental cycles. This study highlights the role of glaciers on the catchment scale chemical weathering rates and emphasizes the poorly understood biogeochemical consequences of waning glacier cover with predicted warming of the Earth's climate.

Methods
Study area. Fieldwork was conducted at the 19 glaciers in Asia, including four glaciers (KOG, DGG, BGG, and HGG) in the Tianshan, six glaciers (SG2, SG3, SG4, YG1, YG5, and LHG) in the Qilian, two glaciers (MKG, and YZG) in the Kunlun, one glacier (DKG) in the Tanggula, three glaciers (GG1, GG2, and GG3) in the Pamir, two glaciers (HLG, and HG2) in the Hengduan, and one glacier (RBG) in the Himalayan mountain ranges (Figs. 1B, and Supplementary Fig. 5 and Tables 1, 2). In this area, the regional climate is influenced by Indian monsoon and westerlies, with limited influence from East Asian monsoon on the Tibetan Plateau. In addition, mineral compositions of proglacial deposits at KOG are dominated by potash feldspar, quartz, and plagioclase (93.4%), followed by illite, calcite, dolomite, kaolinite, and siderite (6.4%), and chlorite, pyrite, and gypsum (0.2%) (Supplementary Table 6). This is similar to the dominance of mineral compositions at DKG in the Tanggula 26 , at HLG in the Hengduan 46 , at Urumqi Glacier No.1 in the Tianshan 45 , and at Qiyi Glacier in the Qilian 60 within Asia.
Sampling and laboratory analysis. Glacier outflows from the 19 Asian glaciers were sampled manually at locations as close to ice margin as possible during Fig. 6 Relationships between cation denudation rate (CDR) and specific discharge (SQ), mean annual air temperature (MAT), mean annual precipitation (MAP), and the latitude (LAT) for glaciers globally. a, d, g, j Glaciers in eight mountain ranges within Asia; b, e, h, k Glaciers in ten glacial regions (excluding ANG; Table 1); c, f, i, l Glaciers at more than 29 glacial basins worldwide (Supplementary Table 10). Note that the horizontal lines above/below the symbols indicate the standard deviation. NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-28032-1 ARTICLE fieldwork that is supported by the PRC Ministry of Science and Technology ( Supplementary Fig. 5). The distance of sampling sites from ice margin ranged from 10 to 100 s of meters, depending on the confluence conditions and the accessibility of meltwater river channels. Bulk meltwater samples (n = 591) were taken every hour over a 1-4 day period in July-October of 2007 or April-August of 2008 to observe the diurnal process of cation concentrations and relevant drainage system in glacial environments (Supplementary Tables 2 and 3). Before sampling, precleaned bottles were rinsed three times in situ by river water. In 2013, glacial deposits (n = 10; diameter is below 5 mm) were sampled as close to the ice margin and river channels as possible in the proglacial area at KOG to characterize their mineral compositions (Supplementary Table 7). After collection samples were immediately frozen and transported to the State Key Laboratory of Cryospheric Science of the Chinese Academy of Sciences for chemical analysis. Samples were immediately filtered through 0.45 μm cellulose nitrate membranes after melting at room temperature (Supplementary Table 1). Concentrations of total dissolved solids were determined with a precision of ±1% using a CON410 meter. Major dissolved cations (Na + , K + , Mg 2+ , and Ca 2+ ) were determined by ion chromatography using a Dionex-600 fitted with an IonPac CS-12A-HC column, using 20 mM MSA eluent and a CSRS-ULTRA-II suppressor. The detection limits for all ions were below 0.85 μeq L −1 and the precision was better than ±1% for all measured ions. The comparative experiments on some samples indicated that the influence of the frozen-melt process after sampling on cation concentrations is negligible. The mineral compositions were determined by X-ray diffraction (D/ Max-IIIB) using a copper butt, step-continuous scanning with 4°/min of speed, and tube 40 kV and 25 mA 26,60 .  Tables 1 and 3). This data set is combined with new data from the 19 glaciers (see above) in eight mountain ranges (ATG, TSG, QLG, KLG, TGG, PAG, HDG, and HMG) within Asia generated in this study to produce a global data set (where cations were not corrected for the sea-salt and aerosol contributions) containing 5465 samples from 77 glaciers (including 14 ice sheet outlet glaciers) worldwide ( Fig. 1A and Supplementary Table 3).
Modeled annual glacial runoff for eight mountain ranges within Asia during 2000-2014 was taken from Wang 82 , and the modeled runoff for other ten glacial regions (excluding ANG; see above) and global glaciers (Table 1) during 2003-2022 were taken from Bliss and others 70 . In brief, annual glacial runoff from Asian mountain ranges was calculated based on a glacier mass balance model and a glacier runoff model. The mass balance model was run on a monthly timescale 89,90 , driven by monthly precipitation and air temperature. For the debris-free glaciers, ice surface ablation was calculated based on a degree-day model. But for the debriscovered glaciers, an ablation influencing factor was introduced to the degree-day model based on a relationship between the debris thickness and sub-ice ice melt rates 91 . The runoff model was proposed by Bliss and others 70 , in which the evaporation, sublimation, and any storage changes (groundwater, englacial water, and subglacial water) were assumed to be negligible 92 . Notably, glacier mass balance is the basis of the runoff model, and all variables involved are integrated over an entire glacier-covered area. The same input data, calibrated parameters for the mass balance model were utilized to drive the glacier runoff model, and the observations from 45 monitored glaciers were used to calibrate related parameters 93 , following the method of Radic and Hock 89,90 . To improve the accuracy of modeling results, the ASTER DEMs-derived geodetic glacier mass balance from Brun and others 94 were used to further constrain the modeled results. In detail, the glacier mass balance and glacier runoff models are briefly illustrated below.
Glacier mass balance model. The monthly-scale mass balance model 89,90 was used in this study to model glacier mass balance for mountain ranges in High Mountain Asia during 1952-2014. In detail, the primary input data involved in this model includes monthly precipitation and air temperature. Glacier area-weighted specific mass balance (B) for the whole glacier in each mountain range was computed as a sum of specific mass balance of each elevation band on a glacier.
where b i and S i denote specific mass balance and glacier area, respectively, and subscripts (i) represent the number of elevation band on a glacier (i = 1, 2, 3, …, n). Monthly specific mass balance (b i ) of each elevation band on a glacier (with an interval of 50 m) was calculated as where a i denotes glacier surface ablation (negative), c i represents glacier mass accumulation (positive), and R i refers to snowmelt refreezing (positive) at each elevation band.
For debris-free glaciers, the glacier surface ablation (a i ) was calculated based on a degree-day model, in which snow/ice melt is considered as a linear correlation to monthly air temperature. Here, monthly a i (mm w.e.) was expressed as where f snow/ice denotes degree-day factor for snow/ice (mm w.e. d −1°C−1 ), and T i denotes monthly air temperature (°C) above the glacier surface. For debris-covered glaciers, the relationship between debris thickness and glacier surface ablation 91 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 k i was introduced to quantify the effect of debris on glacier surface ablation. Thus, monthly glacier surface ablation due to debris cover (a i,debris ) was calculated as where k i denotes the scale factor of each elevation band (i) which depends on the debris thickness.
Monthly glacier mass accumulation for each elevation band c i (mm w.e.) was calculated as where δ m is a constant, T i denotes air temperature at each elevation band, and T snow represents the threshold temperature. If T i is below T snow , P i is assumed to be snow; if not, P i is regarded as liquid precipitation. Based on the relationship between annual potential refreezing R i,pot (cm) and air temperature T a (°C) at each elevation band 95 , R i,pot was calculated as where the lower boundary of annual snowmelt refreezing over the whole glacier is zero, while an upper boundary is applied in the ablation zone and assumed equal to the accumulated snow. Monthly snow meltwater frozen on glacier surface (R i ) does not flow away until accumulated melt in a mass balance year exceeds the annual potential refreezing (R i,pot ).
To interpolate the Climatic Research Unit Time-Series (CRU TS) monthly temperature at each elevation band, two temperature lapse rates (T lap , and G lap ) were utilized. T lap is a parameter similar to a 'statistical lapse rate' between the CRU TS altitude of grid cell with glacier location (h CRU ) and the highest elevation of glacier (h max ), while G lap denotes glacier surface temperature mainly taking into account glacier surface difference, such as orientation and glacier surface climate environment. Monthly air temperature at each elevation band (T i ) was calculated as where T CRU denotes the CRU TS monthly air temperature during the period of 1952-2014, and h represents the altitude of glacier elevation band.
To interpolate the CRU TS monthly precipitation to h max , a precipitation correction factor (k p ) was assigned, while a precipitation gradient (d pre ) was used to interpolate precipitation to each elevation band (the percentage of precipitation decreases with every 50 m decreases in elevation) from the highest to lowest elevation of the glacier. Thus, monthly precipitation at each elevation band was calculated as where P CRU denotes CRU TS monthly grid precipitation for 1952-2014 at the cell of the glacier located.
Glacier runoff model. According to the method proposed by Bliss and others 70 , glacier runoff was calculated using a water balance approach. For individual glaciers in mountain range within High Mountain Asia during 1952-2014, annual glacier runoff (Q i ) at each elevation band on a glacier was calculated as: where a i denotes the melt of ice, firn and snow, as in Eq. (3), which is meltwater from glacier net mass loss. P liq,i represents liquid precipitation as in Eq. (8).
Notably, if T m is above the threshold temperature (T snow ), P i is regarded as liquid precipitation. R i refers to potential refreezing at each elevation band as in Eq. (6). S i indicates glacier area of each elevation band i.
Cation flux estimates. To calculate the annual cation fluxes from eight mountain ranges and nine non-ice sheet glacial regions (Fig. 1) we used mean annual runoff 70,82 multiplied by data set mean Na + , K + , Mg 2+ , and Ca 2+ concentrations to produce annual Na + , K + , Mg 2+ , and Ca 2+ fluxes (Table 1 and Supplementary  Table 5). These were used to generate cation fluxes for each mountain range or glacial region (Table 1). Annual cation flux from all glaciers worldwide was estimated by two methods: (i) the discharge extrapolation of annual Na + , K + , Mg 2+ , and Ca 2+ fluxes from nine non-ice sheet glacial regions based on the percentage of annual glacial runoff from these nine glacial regions to annual glacial runoff globally 70 and then the extrapolated annual Na + , K + , Mg 2+ , and Ca 2+ fluxes were summed, and (ii) the multiplication of regional discharge-weighted mean Na + , K + , Mg 2+ , and Ca 2+ concentrations (as the mean cation concentrations for glaciers worldwide) by annual glacial runoff globally 70 and then annual Na + , K + , Mg 2+ , and Ca 2+ fluxes were summed (Table 1 and Supplementary Table 8). In order to evaluate the effect of discharge-concentration seasonality on estimated cation flux from global glaciers, annual cation flux from global glaciers was also estimated by the discharge-weighted method as described above based on the mid-summer (July-August) data (n = 593) from 24 glaciers worldwide in six non-ice sheet glacial regions (ASG, CEG, ICG, LLG, SCG, and CUG; Supplementary Tables 8 and  9). Similar extrapolations have been adopted in recent estimates of element flux from glaciers and/or ice sheets [10][11][12]26,66 . The approaches are simplistic but provide a first-order estimate of cation fluxes on the mountain range and regional scales. The approaches taken assume that glaciers with available data are representative of entire mountain ranges or glacial regions, which is unlikely given the differences in geology, glacial hydrology, and temporal variation in sampling. Error estimates of cation fluxes were calculated from the standard deviation of data set mean cation concentrations multiplied by annual glacial runoff (Table 1).
Cation denudation rate estimates. To calculate the crustal CDRs we use the annual *Na + , *K + , *Mg 2+ , and *Ca 2+ fluxes divided by glacial area 11,12 (asterisks indicate crustal-derived solutes; Table 1). The annual *Na + , *K + , *Mg 2+ and *Ca 2+ fluxes were estimated by the annual Na + , K + , Mg 2+ , and Ca 2+ fluxes as described above multiplied by the mean contribution (%) of the annual *Na + , *K + , *Mg 2+ , and *Ca 2+ fluxes to the annual total Na + , K + , Mg 2+ , and Ca 2+ fluxes which is calculated based on published data from 2 glaciers in Asia, more than 8 glaciers in Svalbard, 1 glacier in Central Europe and 1 outlet glacier in Greenland for each mountain range or glacial region (Table 1). In detail, for eight mountain ranges in Asia and the Asian glacial region, the mean contributions of the annual *Na + , *K + , *Mg 2+ , and *Ca 2+ fluxes to the annual total Na + , K + , Mg 2+ , and Ca 2+ fluxes were calculated based on the contributions of crustal-derived cations to total cations at Batura Glacier 27  For the GIS and AIS, the mean contributions were calculated based on the crustalderived cation contributions at Kuannersuit Glacier 31 in the Greenland. For the other six glacial regions (ALG, ACG, ICG, LLG, SCG, and CUG) and global glaciers, the mean contributions were calculated based on the crustal-derived cation contributions at all non-ice sheet glaciers (apart from Kuannersuit Glacier) because data are not available in these six regions (Table 1).

Data availability
All data used in this study are present in supplementary materials and have been deposited in the National Tibetan Plateau/Third Pole Environment Data Center (http:// data.tpdc.ac.cn/en/disallow/5a1eb9ee-462a-4840-af81-98caf3738665/).