Glacier evolution in high-mountain Asia under stratospheric sulfate aerosol injection geoengineering

. Geoengineering by stratospheric sulfate aerosol injection may help preserve mountain glaciers by reducing summer temperatures. We examine this hypothesis for the glaciers in high-mountain Asia using a glacier mass balance model driven by climate simulations from the Geoengineering Model Intercomparison Project (GeoMIP). The G3 and G4 schemes specify use of stratospheric sulfate aerosols to reduce the radiative forcing under the Representative Concentration Pathway (RCP) 4.5 scenario for the 50 years between 2020 and 2069, and for a further 20 years after termination of geoengineering. We estimate and compare glacier volume loss for every glacier in the region using a glacier model based on surface mass balance parameterization under climate projections from three Earth system models under G3, ﬁve models under G4, and six models under RCP4.5 and RCP8.5. The ensemble projections

Abstract. Geoengineering by stratospheric sulfate aerosol injection may help preserve mountain glaciers by reducing summer temperatures. We examine this hypothesis for the glaciers in high-mountain Asia using a glacier mass balance model driven by climate simulations from the Geoengineering Model Intercomparison Project (GeoMIP). The G3 and G4 schemes specify use of stratospheric sulfate aerosols to reduce the radiative forcing under the Representative Concentration Pathway (RCP) 4.5 scenario for the 50 years between 2020 and 2069, and for a further 20 years after termination of geoengineering. We estimate and compare glacier volume loss for every glacier in the region using a glacier model based on surface mass balance parameterization under climate projections from three Earth system models under G3, five models under G4, and six models under RCP4.5 and RCP8.5. The ensemble projections suggest that glacier shrinkage over the period 2010-2069 is equivalent to sea-level rise of 9.0 ± 1.6 mm (G3), 9.8 ± 4.3 mm (G4), 15.5 ± 2.3 mm (RCP4.5), and 18.5 ± 1.7 mm (RCP8.5). Although G3 keeps the average temperature from increasing in the geoengineering period, G3 only slows glacier shrinkage by about 50 % relative to losses from RCP8.5. Approximately 72 % of glaciated area remains at 2069 under G3, as compared with about 30 % for RCP8.5. The widely reported reduction in mean precipitation expected for solar geoengineering is unlikely to be as important as the temperaturedriven shift from solid to liquid precipitation for forcing Himalayan glacier change. The termination of geoengineering at 2069 under G3 leads to temperature rise of about 1.3 • C over the period 2070-2089 relative to the period 2050-2069 and corresponding increase in annual mean glacier volume loss rate from 0.17 to 1.1 % yr −1 , which is higher than the 0.66 % yr −1 under RCP8.5 during 2070-2089.

Introduction
High-mountain Asia (HMA) contains the largest number of glaciers outside the polar regions. These glaciers provide water for many large and important rivers (e.g. Brahmaputra, Ganges, Yellow, Yangtze, Indus, and Mekong), and most, but not all, have shrunk over recent decades (Yao et al., 2012). The response of these glaciers to future climate change is a topic of concern especially to the many people who rely on glacier-fed rivers for purposes such as irrigation.
Glacier evolution is expected to be sensitive to climate change. Temperature and precipitation are the important climate factors affecting glaciers. Geoengineering is a method of offsetting the global temperature rise from greenhouse gases, albeit inevitably also altering other important climate parameters, such as precipitation and global atmosphere and ocean circulation teleconnection patterns Ricke et al., 2010). There have been various studies on mountain glacier change under future-climate scenarios, such as A1B and the various Representative Concentration Pathway (RCP) scenarios (Marzeion et al., 2012;Radić et al., 2014;Zhao et al., 2014). In contrast to glaciers in higher latitudes, many on the Tibetan Plateau are summer accumu-lation type (e.g. Fujita and Ageta, 2000); that is, both surface snow fall and melting occur overwhelmingly in the 3 summer months of June, July, and August, with little mass gain or loss throughout the remaining 9 months of the year. However some glaciers, especially in the northwestern parts of HMA, are winter accumulation type . Hence, the glaciers are affected by both the South Asian monsoon system and the westerly cyclonic systems, depending on specific location across the region; thus the region integrates the climate response to two important global circulation systems .
The impact of geoengineering scenarios on ice sheets and glaciers has been limited to studies on global responses based on semi-empirical models (Moore et al., 2010;Irvine et al., 2012) or on simplified ice sheet responses (Irvine et al., 2009;Applegate and Keller, 2015) or implications of a climate model (McCusker et al., 2015), with nothing to date on mountain glacier impacts.
In this paper, we predict glacier area and volume change for every glacier in HMA under projections from six Earth system model (ESM) simulations of climate under the Geoengineering Model Intercomparison Project (GeoMIP) G3 and G4 scenarios . These scenarios envisage use of stratospheric sulfate aerosols to reduce the radiative forcing under the RCP4.5 greenhouse gas scenario during a 50-year period from 2020 to 2069 followed by sudden cessation of geoengineering to determine the "termination effect"  but continued RCP4.5 greenhouse gas forcing for a further 20 years. We address two questions here: (1) would glacier shrinkage and loss in HMA be alleviated under geoengineering by stratospheric sulfate aerosol injection? (2) How would the glaciers respond to the termination of geoengineering?

Study region and glacier data
The Randolph Glacier Inventory (RGI) database contains outlines of almost all glaciers and ice caps outside the two ice sheets (Arendt et al., 2015). Our study region covers HMA (26-46 • N, 65-105 • E), which corresponds to the defined regions of Central Asia, South Asia West, and South Asia East in the RGI 5.0. According to the RGI 5.0, the study region contains a total of 94 000 glaciers and a glaciered area of about 110 000 km 2 . The RGI 5.0 data inside China are based on the second Chinese glacier inventory , which provides glacier outlines from 2006 to 2010, except for some older outlines from the first Chinese glacier inventory where suitable imagery could not be found -mainly in southern and eastern Tibet (the S and E Tibet RGI 5.0 subregion), most of which were made in the 1970s. The RGI 5.0 data outside China are from the "Glacier Area Mapping for Discharge from the Asian Mountains" (GAMDAM) inventory , and nearly all come from 1999-2003, with images selected as close to the year 2000 as possible. Because the data range from each data source is only a few years, we take three reference years -1980, 2009, and 2000 -as start dates for our model simulations of glaciers in S and E Tibet, elsewhere in China, and outside China, respectively.
Following previous authors Zhao et al., 2016), we use median altitude from RGI 5.0 for each glacier as a proxy for equilibrium line altitude (ELA) in the respective initial years: that is, the altitude on the glacier where the local net surface mass balance (SMB) is zero. We use the Shuttle Radar Topography Mission (SRTM) version 4.1 (void-filled version; Jarvis et al., 2008) digital elevation model with 90 m horizontal resolution to estimate the elevation range spanned by each glacier.
Field measurements on SMB are rare in the HMA due to difficulty of access to the glaciers. Following Zhao et al. (2014), we collate SMB versus altitude measurements from 13 glaciers (Table 1 and Fig. 1) to set up parameterizations of mass balance with altitude relative to the ELA for all glaciers.

Statistical model of glacier change
The statistical model for estimating glacier change is based on Zhao et al. (2014Zhao et al. ( , 2016. Briefly the algorithm can be described as follows. We start from known glacier outlines from RGI 5.0 and glacier elevation distribution from SRTM 4.1. In the start year, SRTM DEM data (90 m horizontal resolution) inside the glacier outline are interpolated onto a regular grid with a spatial resolution of 10 m covering the glacier surface. Vertical spacing of altitude bands depends on glacier size, taken as 10 m for glaciers with a total elevation difference from top to bottom larger than 100 m and 1/10 of the glacier altitude difference for glaciers with less altitude range.
We parameterize the annual SMB as a function of altitude relative to the ELA for each glacier. We calculate no more than three SMB gradients using in situ SMB measurements for every glacier in Fig. 1 and Table 1. Following Zhao et al. (2014), the SMB-altitude profile is constructed for every glacier by using its own ELA and these SMB gradients. Where SMB data exist in the sub-region, we use them to parameterize the SMB of all glaciers in that sub-region. Otherwise, we use glaciers from nearby sub-regions.
Integrating the SMB over each glacier gives the volume change rate, which is converted to an area change rate using volume-area scaling (Marzeion et al., 2012): where A(n) is glacier area in the nth year, V (n + 1) is glacier volume and dA(n + 1) is area change rate in the n + 1th year, Figure 1. The HMA region analysed. Sub-regions of the HMA in RGI 5.0 are listed and colour-coded in the legend. Glaciers with SMBversus-altitude measurements (Table 1) are marked with black triangles.
c A = 0.0380 km 3−2γ and γ = 1.290 , and τ A is the response timescale of glacier area and calculated as where L(n) and τ L are glacier length and the response time scale of glacier length in the nth year, respectively. τ L is calculated by following the scaling of Johannesson et al. (1989), where V (n) and P solid (n) denote glacier volume and the solid precipitation on the glacier in the nth year, respectively. The initial glacier length L start is estimated by area-length scaling A start = c L L q L start , where c L = 0.0180 km 3−q (Radic et al., 2008) and q = 2.2 (Bahr et al., 1997). The glacier length change is calculated using the area-length scaling We assume all the area changes take place in the lowest parts of the glacier. The set of glacier surface grid points is updated every year -the number of the grid points that need to be removed or added is calculated using the area change rate, while the elevation of the grid points is updated using SMB. For retreating glaciers, we remove grid cells near the glacier terminus from the glacier surface grids and get the new glacier terminus position and hence the new outline for the next year. For advancing glaciers, we add grid points to the glacier surface grid, whose elevations are all supposed to be the glacier elevation minimum in the n + 1th year, z min (n + 1), which is obtained as follows by assuming a con-stant glacier surface slope, where z max (n + 1) denotes the glacier elevation maximum in the n + 1th year. We also limited the maximal surface increase at any point on the glacier to 15 m above the initial elevation at the starting year. We chose to do this because the valley glacier is physically constrained from growing above the level of the surrounding mountain ridge and side walls. The SMB-altitude profile on each glacier is evolved annually as the ELA changes, and the ELA evolution is estimated by using its sensitivities with respect to temperature and precipitation as follows: where ELA n is the ELA in the nth year from the beginning year; T and P are the inter-annual change of summertime (June-July-August) mean air temperature and annual solid precipitation on the glacier, respectively; and the coefficients α (unit: m • C −1 ) and β (unit: m m −1 ) are the sensitivity of ELA shift to air temperature change ( • C) and precipitation change (m), respectively, which are zonal mean values from energy-balance modelling of glaciers in HMA by Rupper and Roe (2008) (see also Zhao et al., 2014).

Climate scenarios and downscaling of climate data
We run the simulations for glacier change from the relevant start years (Sect. 2) to the year 2089. From the start years to 2013, we use the relatively high resolution monthly-mean gridded 0.5 • × 0.5 • temperature data from the Climatic Research Unit Time-Series (CRU TS) 3.24 dataset (Harris et al., 2014) and 0.5 • × 0.5 • monthly total gridded precipitation data from the Global Precipitation Climatology Centre (GPCC) Total Full V7 dataset (Becker et al., 2013).  (1991( , 1993( , 1994( , ELA varies in 4050-4450 1996( , 1999( , 2001 Ts. Tuyuksuyskiy (43 • 03 N, 3400-4200 0, z > ELA + 100; 1987-2011 WGMS Glacier (13-03) 77 • 05 E) 0.0057, z < ELA + 100; (1991( , 1993( , 1994( , 1996( , ELA varies in 3600-4200 1999( , 2001 (1991( , 1993( , 1994( , 1996 0.01, z < ELA ;1999, 2001, ELA varies in 3950-4175 2007 Haxilegen For the years 2014 to 2089 we use four kinds of climate forcing: experiment RCP4.5 and RCP8.5, and results from two GeoMIP scenarios (G3 and G4; Kravitz et al., 2011) which use stratospheric aerosols to reduce the incoming shortwave while applying the RCP4.5 greenhouse gas forcing. In G3 and G4, stratospheric geoengineering of sulfate aerosol injection starts in the year 2020 and ends in the year 2069. In the 50 years of geoengineering, G3 is designed to achieve a balance between reduction of incoming shortwave radiation and the increase in greenhouse gas forcing, while G4 specifies continuous injection of SO 2 into the equatorial lower stratosphere at a rate of 5 Tg per year from 2020. The across-model spread of temperatures under G4 is larger than under, for example, RCP4.5 (Table 2; there are too few ensemble member models under G3 to see this) because of differences in how the aerosol forcing is handled, and each model has a different temperature response to the combined long-and shortwave forcing . Following the abrupt end of geoengineering, both G3 and G4 specify 20 years of further simulation from 2070 to 2089.
We derived climate forcing data from three climate models participating in G3, five models in G4, and six models in RCP4.5 and RCP8.5 (Table 3). We use the Coupled Model Intercomparison Project Phase 5 (CMIP5; Taylor et al., 2012) output of all models. Yu et al. (2015) noted there was no significant change in surface temperatures after sulfate aerosol was injected in the GISS-E2-R model, possibly due to the efficacy of SO 2 forcing being relatively small as compared to CO 2 forcing in the model. We do not find a termination ef-fect in GISS-E2-R under G3 either. Therefore, we not use any results from GISS-E2-R. We also exclude the model CISRO-Mk3L due to its very coarse spatial resolution of about 4 • and the absence of simulation results in the year 2020 under G4; the models used and their resolutions are listed in Table 3.
Compared with the size of most glaciers in HMA (typically kilometre scale), the CRU, GPCC, and climate model grids have rather coarse resolution (Table 3). The direct use of coarse grid points naturally results in a poor representation of the local climate. Hence we downscale the CRU gridded temperature data, the GPCC gridded precipitation data, and the climate model output to a grid based on a land surface topography having resolution of 0.1126 • × 0.1126 •  using an altitude temperature lapse rate of 0.65 • 100 m −1 , an altitude precipitation lapse rate of 3 % 100 m −1 (Marzeion et al., 2012), and elevation difference of the fine local grid point relative to the climate model grid.
We bias-correct the downscaled model temperatures and precipitation output by using CRU gridded temperature data and GPCC gridded precipitation data as a reference climate. Downscaled series were produced for each climate model for the period 2013 to 2089 under each climate scenario by averaged monthly differences over the baseline period 1980 to 2005 taken from the models' CMIP5 historical simulations. We only use summer (JJA) mean near-surface air temperature. Therefore, future temperature time series T i (t) on each grid point were calculated by where T i,c (t) is monthly-mean temperature for the ith month from the climate model output from t = 2013 to 2089, and T i,CRU and T i,c,history are mean temperature from the CRU TS V3.24 dataset and climate model output, respectively, for the ith month averaged over the baseline period 1980-2005 on each grid point. Future precipitation time series P i (t) on each grid point were calculated by where P i,c (t) is monthly precipitation for the ith month from the climate model output from t = 2013 to 2089, and P i,GPCC and P i,c,history are monthly precipitation from the GPCC dataset and climate model output, respectively, for the ith month averaged over the baseline period 1980-2005 on each grid point. The temperature and precipitation on each glacier were calculated by an altitude temperature lapse rate of 0.65 • 100 m −1 , precipitation lapse rate of 3 % 100 m −1 (Marzeion et al., 2012), and the elevation difference of the glacier surface elevation relative to the nearest fine grid point. Moreover, the solid precipitation on the glacier is calculated by the fraction of solid precipitation, f solid , based on the monthly-mean temperature T a on the glacier as (Fujita and Nuimura, 2011)

Validation of the glacier model and methodology
In this section we justify the selection of various parameter values used in the method here. In Sect. 5 we indicate how elements in the model and climate forcing affect the uncertainties of the results we produce in Sect. 4, and how those results compare with previous estimates of glacier evolution in HMA. A crucial parameterization concerns the SMB-altitude gradients. The field data (Table 1) include three more glaciers than those used in Zhao et al. (2014Zhao et al. ( , 2016 and include a benchmark glacier from almost every sub-region. With so few glacier observations available, there is an issue of how representative they are of the general population. For inner Tibet, there are three glaciers (Zhadang, Gurenhekou, and Xiao Dongkemadi glaciers) with SMB observations, and they have almost the same SMB-altitude gradients, 0.0041 m m −1 , over their common elevation range (5515-5750 m, Table 1); two glaciers (Naimona'nyi and Kangwure) in central Himalaya (C Himalaya) have SMB gradients of 0.0038 m m −1 in their common altitude range of 5700-6100 m. These similarities suggest that the measured glaciers share some important characteristics with the vast majority which are not surveyed.
Next we consider the choices for the initial value of ELA at the start year, different V-A scaling parameters, and different ELA sensitivities to summer mean temperature and annual precipitation.
In choosing the initial ELAs for each glacier, there are several reasonable alternatives (Zhao et al., 2016): (i) using ELAs interpolated from the first Chinese glacier inventory, (ii) using median elevations from the RGI dataset, and (iii) using the elevation of the 60th percentile of the cumulative area above the glacier terminus. Zhao et al. (2016) showed that these three choices lead to a range of about 2.5 mm of global sea level in glacier volume loss by 2050. In this study, we use median elevations from the RGI dataset, which corresponds to the median result. Zhao et al. (2014) showed that different volume-area scaling parameterizations can lead to a ±5 % range of glacier volume loss. The set of parameters we use in this study corresponds to the lower bound of estimated volume loss, but one that is best matched to the observational dataset of 230 separate glaciers .
For the ELA sensitivity to summer mean temperature and annual precipitation, we use the zonal mean values from energy-balance modelling of glaciers in HMA by Rupper and Roe (2008). Alternatively, it can be estimated using an empirical formula for ablation and a degree-day method (Zhao et al., 2016). Zhao et al. (2016) calculated the ELA for nine glaciers in China, India, and Kyrgyzstan, and compared them with the observed ELA time series by similarities of decadal trends and also annual variability. The Rupper and Roe ELA parameterization produced the best fits to observed ELA decadal trends on nine glaciers, with a correlation coefficient of 0.6, which is significant (p < 0.05; the values we give for p are single-tailed Pearson correlation tests).
Combining the above uncertainties would require a Monte Carlo simulation since the parameters combine non-linearly to produce glacier volume and area change; this is prohibitively expensive to perform given that a single simulation of all glaciers in HMA requires about 60 CPU hours on an eight cores computer with parallel computing in Matlab. We did estimate elevation changes for individual glaciers directly from simulated volume and area changes; then we calculated the average rate of elevation change for all the glaciers in each sub-region and compared them with remote-sensing estimates from 2003 to 2009 from Gardner et al. (2013) (Table 4). The correlation coefficient between the Gardner et al. (2013) estimates for the 6 RGI 5.0 sub-regions with regional data and our modelled regional averages is 0.7, which is marginally significant (p < 0.1).
In our simulations we have used constant lapse rates for temperature (0.65 • 100 m −1 ) and precipitation (3 % 100 m −1 ). To check how reliable this is, we chose five We construct the climate forcing by using CRU temperature data and GPCC precipitation data before 2013; climate models (Table 3) Figure 3 shows the time series of JJA mean temperature and annual precipitation forcing from the beginning years to 2089, with the across-model range from the ensemble members; ranges found are slightly smaller than the regional spread found by Yu et al. (2015) due to the grid-point-by-grid-point bias correction we apply here. The multi-model mean temperature under G4 is higher than that under G3 in the geoengineering period. In contrast with the ensemble mean temperature, the HMA mean temperature projected by HadGEM2-ES under G4 is cooler than that under G3, and its G4 temperature is lower than the ensemble mean, while its G3 is higher than the ensemble mean (Fig. 2). The across-model spread in temperature response to G4 is larger than that under G3. Temperatures projected by BNU-ESM are lower than the ensemble mean under both G3 and G4.
The temperature averages over the whole region and the glaciated parts have similar trends. Temperatures under RCP8.5, as expected, increase at the highest linear rate among all the scenarios. Temperature rise under RCP4.5 is the next highest, and its rate becomes smaller after about the year 2050 as specified greenhouse gas emissions decline. There is relative cooling of 1.05 • C under G3 and 0.76 • C under G4 as compared with RCP4.5 during 2020-2069 across the whole region (Fig. 3). Yu et al. (2015) noted that G3 produced relative cooling under G3 of 0.58 • C and under G4 of 0.53 • C in globally averaged temperature over the 2030-2069 period.
There is no trend in temperature under G3 over the geoengineering period (2020-2069). But after the termination in the year 2069, there is a temperature rise of about 1.3 • C over the period 2070-2089 relative to the period 2050-2069 under G3. There is a less obvious termination rise of temperature under G4 than under G3. This is due to G4 having a constant stratospheric aerosol injection rate of 5 Tg SO 2 yr −1 , while G3 gradual ramps up the aerosol so that about twice as much is needed by 2069, depending upon the sensitivity of the particular model to stratospheric sulfate aerosols. Hence, the radiative impact of terminating G3 is about twice as large as terminating G4, and the termination temperature signal is much more obvious in G3 than G4.
The annual precipitation averages over the whole region do not show obvious trends in any climate scenarios (Fig. 3c). However, the annual solid-precipitation averages over the glaciers show decreasing trends in all the scenarios (Fig. 3d) until 2070, which is due to increases in surface air temperature (Fig. 3b). Under RCP8.5, annual solid precipitation averaged over each glacier decreases fastest (2.2 mm yr −1 ). Decreases are similar (about 1.5 mm yr −1 ) under RCP4.5 and G4 and least (0.3 mm yr −1 ) under G3 during the geoengineering period (2020-2069). After the year 2070, there are no trends in annual solid precipitation under G3, G4, and RCP4.5 (Fig. 3d) due to stable temperatures (Fig. 3b).

Glacier changes across HMA
Glacier volume changes for all the glaciers in the study region computed using temperature and precipitation data from the four scenarios are shown in Fig. 4a. The uncertainty we plot is due only to the differences between climate forcing across the models; it does not reflect uncertainty of the glacier model parameters. Volume loss rates increase in the following order, from lowest to highest, for the period 2020-2069: G3, G4, RCP4.5, RCP8.5. The RCP4.5 and RCP8.5 scenarios produce similar continuous mass loss until approximately 2035 (Fig. 4a) mainly due to the similarity of temperatures projected by RCP4.5 and RCP8.5 in the period 2020-2035 (Fig. 3a), and both show relatively slower loss rates after about the year 2050 probably because the most sensitive glaciers have already retreated before 2050. The multi-model mean glacier volume loss in equivalent to sea-level rise for the whole study region from 2010 to the end of geoengineering in 2069 is 9.0 mm (G3), 9.8 mm (G4), 15.5 mm (RCP4.5), and 18.5 mm (RCP8.5), with 91.8, 96.0, 98.5, and 99.7 % of glaciers retreating under these scenarios (Table 5). Volume  (Table 3).  region (b, d). Note the different temperature ranges in (a, b) and precipitation ranges in (c, d). Precipitation in (d) is the average annual solid precipitation at the ELA of each glacier in the start year of simulations, which is taken here to be representative of each glacier. The solid curves and shadings from 2013 to 2089 are the ensemble mean and the across-model spread between ensemble members for each scenario, respectively, which are colour-coded in the legend. loss using the climate projected by HadGEM2-ES under G4 is far less than that by other models (Table 5), this produces a larger standard deviation for the results than for other scenarios in Table 5. The cause is the combination of small precip-itation decrease under RCP4.5 and the G4 anomaly, accompanied by only modest warming (Table 2). These numbers may also be compared with the simulations run using the ensemble mean climate forcing (last row in Table 5), which are Table 5. The volume loss in millimetre sea-level equivalent, projected using forcing from all the climate models in the period 2010-2069 and the post-geoengineering period 2070-2089 under G3, G4, RCP4.5, and RCP8.5. The means of volumes lost driven by individual model forcing and its standard deviation are shown in the penultimate row. The simulated volume loss using the climate model ensemble mean forcing of temperature and precipitation is shown in the last row. The volume loss is calculated by assuming ice density of 900 kg m −3 and ocean area of 362 × 10 12 m 2 .  all close to the means of the individual model-driven mass losses, as are the time-varying loss rates (Fig. 4). Therefore, the geoengineering schemes G3 and G4 help to reduce glacier mass loss in our simulations, and G3 reduces glacier loss more than G4, which is due to stronger temperature cooling effect under G3 (Sect. 4.1.1).
There is a clear increase in volume loss rate under G3 after 2069, when geoengineering is terminated. Comparing the last 15 years of geoengineering (2055-2069) with the first 15 years post-geoengineering (2070-2084) shows annual mean volume loss rate for all the glaciers of 0.17 % yr −1 (referenced to the volume in the year 2010) increases to 1.11 % yr −1 , which is higher than the annual mean volume loss rates of 0.54 % yr −1 for RCP4.5 and 0.66 % yr −1 for RCP8.5 in the period 2070-2084. However, the volume loss rate under G4 shows negligible termination effect; annual mean volume loss rates change from 0.73 before to 0.86 % yr −1 after the termination. The glacier volume loss over the post-geoengineering period of 2070-2089 for both G3 and G4 is higher than for either RCP4.5 or RCP8.5 (Table 5). However, by 2070 under both RCP scenarios there is much less glacier ice volume remaining than under G4, or especially G3. When comparing ice loss rates at comparable total volumes, loss rates with RCP8.5 are similar to those of post geoengineering G3.
As may be expected, glacier area change trends under each climate scenario are similar to the volume change trends (Fig. 4b). We project 53, 41, 27, and 14 % of the area in 2010 remaining in the year 2089 under the G3, G4, RCP4.5, and RCP8.5 scenarios, respectively.
L. Zhao et al.: Glacier evolution in HMA under geoengineering 4.2 Sub-regional climate and glacier changes 4.2.1 Sub-regional temperature and precipitation change There are three RGI 5.0 regions in HMA: Central Asia, South Asia West, and South Asia East. They are named as Region 13, 14, and 15 and sub-divided into smaller sub-regions in the RGI 5.0 dataset (Fig. 1). In this section we plot the averages of JJA mean temperatures (Fig. 5) and that of annual solid precipitation at the ELA of every glacier in the start year ( Fig. 6) in every sub-region under all the climate scenarios. Temperatures under RCP8.5 increase at the highest rates (0.053-0.087 • C yr −1 ) among all the scenarios, with temperature increases under RCP4.5 in the range of 0.030-0.059 • C yr −1 , with its rate decreasing after about the year 2050 as specified greenhouse gas emissions decline. The temperatures rise of 0.030-0.050 • C yr −1 occurs under G4 across the sub-regions. Under RCP4.5, RCP8.5, and G4, temperatures increase most slowly in the southeast of the study area (S and E Tibet, C Himalaya, E Himalaya, and Hengduan Shan) and fastest in the northwest (mainly Tien Shan, Hissar Alay, Karakoram, Pamir, and Hindu Kush).
There is no trend in temperature under G3 in the geoengineering period  in all the sub-regions. The temperature cooling projected by G3 as compared with RCP4.5 during 2020-2069 is about 1.0 • C in sub-regions of Central Asia, 1.2 • C in South Asia West, and 0.8 • C in South Asia East (Fig. 5). After the termination in the year 2069, there is temperature rise of about 1.07-1.65 • C over the period 2070-2089 relative to the period 2050-2069 under G3. The post-termination temperatures increase the least (about 0.020 • C yr −1 ) in Karakoram and the most (about 0.046 • C yr −1 ) in eastern Kunlun. The temperature cooling projected by G4 as compared with RCP4.5 during 2020-2069 is very similar across all the sub-regions: 0.68-0.86 • C.
The annual solid precipitation projected by RCP8.5, and to a lesser degree by RCP4.5 and G4, decreases in all the subregions, with the rates larger than 1.5 mm yr −1 in S and E Tibet, Hindu Kush, W Himalaya, and the whole of Region 15 (C Himalaya, E Himalaya, Hengduan). There is no obvious trend of solid precipitation projected by G3 in the geoengineering period (2020-2069) in most sub-regions. But after the geoengineering termination under G3 in the year 2069, there is a significant decrease of solid precipitation in S and E Tibet, Hindu Kush, and the whole of Region 15.

Sub-regional glacier changes
Glacier volume changes in the HMA sub-regions are shown in Fig. 7. Glacier volumes in all the sub-regions decrease during the period 2020-2089, with the highest rates under RCP8.5 and the second-highest rates under RCP4.5, as expected. Glacier volumes decrease with lower rates under G3 and G4 in all the sub-regions except S and E Tibet, inner Tibet, and Hengduan Shan, where glacier volumes increase from the year 2020 to about 2040 under G4 and to the end of geoengineering period under G3 (Fig. 7). The glacier volume triples in S and E Tibet and increases by about 56 % in inner Tibet, while increasing slightly in Hengduan Shan in the geoengineering period under G3. The "termination effect" of geoengineering under G3 is significant in most sub-regions.
There are some noticeable difference between means of individual climate-model-forced simulations and the results using multi-model ensemble mean climate forcing (Fig. 7): for example, S and E Tibet under all the scenarios, Karakoram under G3, and inner Tibet under G4. This could be because (i) individual model differences in temperature and precipitation forcings are large between ensemble members and their means (especially for the three-model ensemble in G3) in particular sub-regions; (ii) glacier hypsometry differences between regions lead to sensitivity under some combinations of forcing when the ELA change is located around large amounts of ice; and (iii) glacier data inside S and E Tibet were measured in the 1970s (Sect. 2) and contain outlines of glacier complexes rather than individual glaciers, which has an impact on the volume estimate because of the nonlinearity of volume-area scaling relationship. Furthermore, the observations offer some support to the model simulations. Liu et al. (2006) and Shi et al. (2006) found that over 40 % of the glaciers in the Gangrigabu Mountains of the S and E Tibetan Plateau have been advancing since the mid-1980s, which is a peculiar phenomena and due to the increase of high precipitation brought by the Indian monsoon. The across-model spread for each sub-region is not shown for clarity. Note the difference of glacier volume ranges in the panels.

Uncertainties in projections
Glacier model parameter selection was discussed in Sect. 3.3 and is discussed in more depth by Zhao et al. (2016). In this section we address, and try to estimate, how systematic errors in climate forcing or glacier model parameters cause errors in projections of HMA glacier contributions of sea-level rise.

Climate forcing
There are several uncertainties in climate model forcing used to drive the glacier model in this study. The models are also relatively coarsely gridded, certainly as compared with the vast majority of glaciers, and so differences may be expected between statistically downscaled forcings based on lapse rates that we use here and those produced from highresolution dynamic climate model forcing.
Firstly, only three ESMs participated in G3, while five participated in G4, simply because doing the G3 experiment is difficult and time-consuming to set up. So the ensemble climate projection by G3 is less robust than that by G4. In many Figure 5. JJA-mean bias-corrected surface air temperature time series from 2010 to 2089 in the sub-regions of Region 13 (left column panels), 14 (middle column panels), and 15 (right column panels) under scenarios by row: G3 (top panels), G4,RCP4.5,and RCP8.5 (bottom panels). Note the different temperature ranges in the panels. The curves and shadings are the ensemble mean and the spread between ensemble members for sub-regions, respectively, which are colour-coded in the legend.
cases it seems that the results from G3 and G4 are statistically similar enough to be combined Moore et al., 2015). We tested the differences between RCP8.5, RCP4.5, and G4 using the four models in common (Table 5) and found the glacier responses are significantly different (p < 0.05). Although there are too few models in common between G3 and G4, the dominant influence of summer melting to the mass balance across the region (Zhao et al., 2016), and the clear difference in temperature across HMA between G3 and G4 (Figs. 2 and 3) suggests the glacier response in HMA is different between G3 and G4.
Secondly, although the goal of geoengineering schemes is to mitigate temperature rise, it inevitably also alters other important climate parameters, such as precipitation. Simulating change in the Asian monsoon is difficult for climate models under geoengineering since the deep convection involved may also be influenced by chemistry changes in the stratosphere caused by the injected aerosols -most of the ESM models in our study do not have sophisticated aerosol chemistry schemes (though the MIROC-ESM-CHEM model does). Tilmes et al. (2013) showed that changes under the G1 scenario (which specified a much larger shortwave radiation reduction than G3 or G4) produced a weakening of the Asian monsoon and the hydrological cycle by about 5 %. The reductions in solid precipitation (Fig. 3) under RCP8.5 are about one-third of historical levels, and the regions most affected in Region 15 (Fig. 6) are some of those most influenced by monsoon precipitation patterns (Fig. 1). Hence the temperature-driven shift from solid to liquid precipitation is probably more significant than changes in monsoon precip- itation suggested by the G1 results discussed by .
Thirdly, we note that the distribution of meteorological stations in the study region is very sparse, especially in the northwest of this region (Liu and Chen, 2000). Therefore, both the CRU gridded data and data from model projections that we used in this study may have low accuracy for specific glacier regions. This has also implications for the use of very high resolution dynamic models; one such model simulated air temperatures and down-welling radiative fluxes well, but not wind speed and precipitation, producing unstable results when used with the CLM45 land model that simulated ground temperatures and snow cover (Luo et al., 2013). Explicit glacier atmospheric mass balance modelling , a technique based on very high spatial and temporal resolution climate data (hourly and 60 m), was used on Zhadang Glacier (Fig. 1, Table 1), with in situ observations available, but not across the general expanse of the glaciated region; this study also noted the importance of wind speed to glacier mass balance in the region influenced by the Indian monsoon. Maussion et al. (2014) demonstrate that 10 km resolution dynamic modelling of the region can be done successfully and potentially can improve the precipitation modelling over the statistical downscaling methodology we employ here, though to date this is a reanalysis dataset with no prognostic simulations. Zhao et al. (2014Zhao et al. ( , 2016) used a 25 km resolution regional climate model, RegCM3, to drive their simulations of glacier response to scenario A1B. By 2050 under A1B (which is intermediate between RCP4.5 and 8.5 in temperature rise), a sea-level rise equivalent to 9.2 mm was projected from HMA. In comparison, our estimates are 11.1 and 12.5 mm for RCP4.5 and RCP8.5, respectively (Fig. 4).

Glacier model
The model we use is not particularly sophisticated; it simply relies on statistical relationships between mass balance and ELA. Compared with the method used in our previous studies (Zhao et al., 2014(Zhao et al., , 2016, we improved our method here by considering the area response time in the volumearea scaling (Eq. 1), which is more physical. We also allow the glacier area to grow (Sect. 3.1), giving better estimates of glacier area for advancing glaciers. The motivation to use a relatively simple model must be that it simulates the glaciers well given the available data. As previously discussed, there is a shortage of observational data both on glaciers and from climate stations across HMA. In Sect. 3.3 we discussed how the model performs when tested against by the limited data available from satellites and ground measurements; in this section we compare the model against previous simulations of HMA glaciers under climate warming and examine how its weaknesses may affect the reliability of projected mass changes.
Perhaps a strong limitation on the glacier simulation under geoengineering in our model is the lack of response to the changes in shortwave forcing that would be produced under aerosol injection schemes. Van de Berg et al. (2011) showed that Greenland mass balance during the Eemian interglacial could not be explained purely by temperature rise but must also include losses due to changes in the shortwave radiation flux on the ice sheet.
Testing our results for the greenhouse gas scenarios against previous studies, we project glacier volume loss, in equivalent sea-level rise, for all the glaciers from 2010 to 2089 as 18.2 ± 2.5 and 22.4 ± 1.3 mm under the RCP4.5 and RCP8.5 scenarios, respectively. The volume change of all the glaciers in HMA over the 21th century estimated by Radić et al. (2014) is about 15 ± 5 mm under RCP4.5 and 22 ± 5 mm under RCP8.5. Marzeion et al. (2012) estimate about 15.4 ± 4.5 mm under RCP4.5 and 18.8 ± 4.0 mm under RCP8.5 using projected temperature and precipitation anomalies from an ensemble of 15 CMIP5 climate models. The results projected by our method have higher means than theirs but do not differ significantly.
The across-model uncertainties we plot here (Fig. 4) are smaller than glacier method uncertainties (Sect. 3.3; Zhao et al., 2016). Hence, more mass balance and meterological stations on glaciers across the region, or longer and higherspatial-resolution time series of glacier elevation changes, would better constrain the projected mass losses than simply increasing the number, or resolution, of climate models used in the simulations. That is, the range of mass projections given by the mass balance model with different, but reasonable, choices of data-limited quantities, such as the ELA sensitivity to temperature or the SMB-altitude gradients, is larger than the across-model range for each climate scenario.

Summary and conclusion
We estimate and compare glacier volume loss for glaciers in HMA using a statistical model based on glacier SMB parameterization to the year 2089. We construct temperature and precipitation forcing by using CRU temperature data and GPCC precipitation data before 2013, and projections from 6 Earth system models running RCP4.5 and RCP8.5 and the stratospheric sulfate aerosol injection geoengineering scenarios G3 and G4 with model bias correction and downscaling to a high-resolution spatial grid based on fixed altitudinal lapse rates for temperature and precipitation. In assessing how glaciers respond to geoengineering climates, we consider only across-climate model differences between the scenarios rather than uncertainties in glacier mass caused by errors in the glacier model we use. The projections suggest that glacier shrinkage at the end of the geoengineering period in 2069 is equivalent to sea-level rise of 9.0 ± 1.6 mm (G3), 9.8 ± 4.3 mm (G4), 15.5 ± 2.3 mm (RCP4.5), and 18.5 ± 1.7 mm (RCP8.5) relative to their volumes in 2010 (Table 5), with 91.8, 96.0, 98.5, and 99.7 % of glaciers retreating under these scenarios, respectively. There are clear increases in temperature and glacier volume loss rate under G3 after 2069, when geoengineering is terminated, which is higher than the rate under RCP8.5. But the termination effect under G4 is negligible. Glacier volumes decrease in most sub-regions under all the scenarios, while they increase in inner Tibet, S and E Tibet, and Hengduan Shan from the year 2020 to about 2040 under G4 and to the end of geoengineering period under G3.
Although G3 keeps the average temperature from increasing in the geoengineering period, G3 only slows glacier shrinkage by about 50 % relative to losses from RCP8.5. Approximately 72 % of glaciated area remains at 2069 under G3, as compared with about 30 % for RCP8.5. The reason for the G3 losses is likely to be that the glaciers in HMA are not in equilibrium with present-day climate, so simply stabilizing temperatures at early-21st-century levels does not preserve them. To do that would require significant cooling, perhaps back to early-20th-century levels. Achieving that cooling by sulfate aerosol injection may not be possible. The 5 Tg of SO 2 per year specified in G4 is about the same loading as a 1991 Mount Pinatubo volcanic eruption every 4 years (Bluth et al., 1992). G3 requires increasing rates of injection, to 9.8 Tg, for the BNU-ESM in 2069. As aerosol loading increases, its efficacy decreases as particles coalesce and fall out of the stratosphere faster, while also becoming radiatively less effective (Niemeier and Timmreck, 2015). This effect is so strong that it appears unfeasible to use sulfate aerosols to completely eliminate warming from scenarios such as RCP8.5. Greenhouse gas emissions would require very drastic reduction from present levels, and net negative emissions within the next few decades, to limit global temperature rise to 1.5 or 2 • C (Rogelj et al., 2015). If such targets are met, then it is conceivable that plausible quantities of sulfate aerosol geoengineering may be able to maintain 2020 temperatures throughout the 21st century. Our simulations suggest that, even if this politically very difficult combination of drastic emission cuts and quite aggressive sulfate aerosol geoengineering were done, the disappearance of about one-third of the glaciated area in HMA by 2069 still could not be avoided.