Insensitivity of mass loss of Icelandic Vatnajökull ice cap to solar geoengineering

. Geoengineering by stratospheric aerosol injection (SAI) may reduce the mass loss from Vatnajökull ice cap (VIC), Iceland, by slowing surface temperature rise, despite relative increases in ocean heat flux brought by the Atlantic Meridional Circulation (AMOC). Although surface mass balance (SMB) is affected by the local climate, the sea level contribution is also dependent on ice dynamics. We use the Parallel Ice Sheet Model (PISM) to estimate the VIC mass balance under the CMIP5 (Coupled Model Intercomparison Project Phase 5) RCP4.5, 8.5 and GeoMIP (Geoengineering Model Intercomparison Project) 15 G4 SAI scenarios during the period 1982-2089. The G4 scenario is based on the RCP4.5, but with additional 5 Tg yr -1 of SO 2 injection to the lower stratosphere. By 2089, G4 reduces VIC mass loss from 16 % lost under RCP4.5, to 12 %. Ice dynamics are important for ice cap loss rates, increasing mass loss for RCP4.5 and G4 by 1/4 to 1/3 compared with excluding ice dynamics, but making no difference to mass loss difference under the scenarios. We find that VIC dynamics are remarkably insensitive to climate forcing partly because of AMOC

throughout this century (Schmidt et al., 2017;Schmidt et al., 2020;Aðalgeirsdóttir et al., 2020). Although their contribution to global mean sea-level rise would be just 1 cm, even if all the ice melted (Björnsson and Pá lsson, 2008), the local impacts of rapid glacier loss will be obvious and deeply moving for Icelanders and cause profound changes in hydrology.

35
Vatnajökull ice cap (VIC), the largest nonpolar ice cap in Europe with a volume of 2870 km 3 and an area of 7700 km 2 (Aðalgeirsdóttir et al., 2020), is shrinking at an increasing rate (e.g., Björnsson et al., 2013). Surface mass balance (SMB, the sum of accumulation and ablation) significantly decreased from a slightly positive balance in the 1980s to -0.8 m yr -1 during 1995. Projections show that VIC volume will decrease 51-94 % by the year 2300 under the representative concentration pathway (RCP) 8.5 scenario (Schmidt et al., 2020). 40 Geoengineering by stratospheric aerosol injection (SAI) may help to meet IPCC (Intergovernmental Panel on Climate Change) targets of limiting global mean temperature rises to less than 2℃, since present greenhouse gas emissions trajectories suggest temperatures over much of the latter half of the 21st century will exceed this limit (MacMartin et al., 2018). Moreover, SAI may be a relatively cheap way to offset temperature rises on the global scale (Smith and Wagner, 2019). SAI does reduce 45 surface runoff over VIC compared with RCP4.5 and RCP8.5 (Yue et al., 2021), but by relatively little compared with the SAI induced reductions in mass loss from both Greenland (Moore et al., 2019) and smaller mountain glaciers (Zhao et al., 2017).
However, Yue et al. (2021) did not consider non-surface mass balance generated by changes in ice flow and discharge (e.g., calving of ice and basal melting) that are driven by changing climate or impacts due to the warming ocean in contact with the ice. It is this component that we tackle here. Although there are many other potential impacts that might be expected if SAI 50 were ever undertaken, here we focus only on the mass balance of a single ice cap in Iceland. While this is a very narrow focus, the study is topical and of wider relevance because of the atypical behaviour of the North Atlantic under both greenhouse gas forcing and SAI, primarily because of the compensatory effects the climate forcing has on the AMOC. The AMOC is simulated as slowing considerably under greenhouse gas forcing (Cheng et al., 2013), while slowing much less under geoengineering (Hong et al., 2017;Yue et al., 2021). Furthermore, the Arctic is currently warming at least twice as fast as the global mean 55 temperature, leading to concerns on the stability of the Arctic cryosphere and examination of possible roles for geoengineering methods its preservation (Lee et al., 2021). Thus, the North Atlantic Arctic sector is a key region for examination of unwelcome impacts from geoengineering.
We simulate the response of the VIC with the state-of-the-art Parallel Ice Sheet Model (PISM) driven by monthly surface mass 60 balance from 2006-2089 under SAI and greenhouse gas climates. The RCP4.5 emissions scenario (Thomson et al., 2011) is close to future emissions under the 2015 Paris climate agreement (Kitous and Keramidas, 2015), while the RCP8.5 scenario (Riahi et al., 2011) is a "business-as-usual" extreme failure to mitigate scenario. The SAI G4 scenario branches off the RCP4.5 scenario at 2020, specifying 5 Tg yr -1 of SO2 to be injected into the equatorial lower stratosphere until 2069, and then continues with RCP4.5 forcing to 2089(Kravitz et al., 2013.

Model and Verification
The PISM model (version 1.0; Bueler and Brown (2009); https://pism-docs.org/wiki/doku.php) is an open-source ice sheet thermo-dynamic model that has been used in numerous studies of a wide range of ice sheets and glaciers (e.g., Aschwanden et al., 2019;Yan et al., 2020). PISM includes a hybrid stress balance model (Bueler and Brown, 2009) with both Shallow Ice Approximation (SIA; Hutter, 1983) and Shallow Shelf Approximation (SSA; Morland, 1987) to solve ice vertical deformation 70 and longitudinal stretching allowing simulation of both slowly flowing ice cap interiors and fast flowing outlet glaciers. Schmidt et al. (2020) optimized the free parameters in PISM using observations over Vatnajökull, and we use these choices for all our simulations, with the SIA+SSA hybrid schemes for ice stress balance, the isothermal Glen flow law for ice flow, the pseudo-plastic flow law for ice basal sliding, and the Eigen scheme for ice calving.

75
To initialize PISM over the VIC, we need the boundary conditions of the surface elevation, bedrock altitude, ice thickness, upward geothermal flux, ice temperature, and surface mass balance (Fig. 1). We re-grided these to 500 m×500 m resolution.

90
We drive the PISM model with monthly SMB fields from four ESMs during 1982-1999, when the ice cap was close to steady-state, repeated for 2000 years (Fig. 2), achieving a near steady-state ice cap geometry. The final year of spin-up is then used as the initial condition in simulations.  During spin-up, VIC volume changes by between -1.3% and 0.8% for the four ESMs relative to present day, while area is reduced by around 16% (Fig. 2). Ice area loss is mainly over the outlet glaciers of Dyngjujökull, Brúarjökull and Sí ðujökull and are consistent across all the ESMs because their climate variables (e.g., surface air temperature, downward longwave and shortwave radiation) that drive SMB are bias-corrected with ERA5 reanalysis. This tends to impose a fixed spatial pattern to SMB while preserving the long-term trends in the ESM forcing fields, which are small during the 1982-1999 period. The largest discrepancies between the spin-up and present area for VIC, especially the reductions over the northern glaciers and increased thickness over Skeiðará rjökull, are likely caused by the quasi-periodic surging glaciers that are not parameterized 120 by PISM, but which cover about 75% of the VIC . VIC surges occur on timescales from several years to centuries, and rapidly move accumulated ice from the upper glacier towards the terminus, meaning that the upper and lower glacier are usually out of steady state. Thus, the spin-up is unlikely to achieve a present-day area coverage, although total volume is close to observed.

Ice cap SMB, MB and non-SMB from 2020-2089
In Fig. 6 we separate the SMB and non-SMB (ice dynamics and basal melting) components of overall mass balance. SMB increases the ice thickness over the interior of VIC, with the maximum of more than 400 m over the southern region of VIC, while decreasing the ice thickness over the margins, especially over the outlet glaciers of Skeiðará rjökull and Breiðamerkurjökull. The degree of surface thinning under RCP8.5 is stronger than for the RCP4.5 and G4 scenarios, resulting 165 in the smallest area of surface thickening. Non-SMB components display the opposite pattern to SMB. Positive non-SMB contributions are visible in all ESMs and scenarios over the margins, because as the negative SMB steepens the margins, it is compensated by increased ice flow from the interior and so dynamically thickens. Similarly, the drawdown of ice from the interior thins the higher elevation ice cap making the non-SMB component there negative. Basal melting is driven by nonclimate factors and so remains essentially unchanged under the scenarios. The pattern of non-SMB contributions for individual 170 ESM are all quite similar, the largest differences being mainly over the ablation zone, with ensemble standard deviations of more than 10 m (Fig. S2). Fig. 6b demonstrates that surface height differences (G4-RCP4.5 and G4-RCP8.5) by 2089 are mainly caused by SMB rather than dynamic effects, SMB under G4 increases VIC mean surface height by around 20 m. SAI dynamically thickens the ablation zone relative to the RCP scenarios, while thinning the accumulation area. The dynamic impact on surface height change differences between G4 and RCP4.5 are confined to ±10 m, which is much less than 175 differences from RCP8.5. Fig. 6c illustrates the across-ESM mean time series of modeled MB, SMB and non-SMB. MB is well-correlated (R 2 =0.98, p<0.01) with SMB, and generally exhibits a downward trend in all ESMs and scenarios, although with high annual variability.
The non-SMB contributions, however, remain nearly constant (around -0.25 m yr -1 ) over time and across all scenarios (Fig.  180   7). By 2069, ice dynamics are important for ice cap loss rates, increasing mass loss for RCP4.5 and G4 by 1/4 to 1/3 compared with excluding dynamics, but making no difference to mass loss difference under the scenarios.  N=4) compared with RCP4.5 and RCP8.5, although two MIROC models project almost no MB differences between G4 and RCP4.5. The mean SMB and MB under G4 have much larger differences between the ESMs than the two RCP 195 scenarios (Fig. 7). During the historical period, our simulations show the overall mass loss on VIC is about equally divided between SMB and non-SMB components, but as SMB becomes more negative, the non-SMB fraction becomes less important. Yue et al. (2021) simulate that SAI G4 (2020-2069) could reduce surface runoff by 6.2±6% (95% confidence intervals), relative to RCP4.5.
Including the dynamic and basal melt components leads to 12±2% mass loss under G4, and 16±4% under RCP4.5 of present-205 day mass by 2089. So, G4 saves 25% of mass lost relative to RCP4.5, while for High Mountain Asia (HMA) it saves 19% more than RCP4.5 (Zhao et al., 2017). The differences in efficacy are related to the degree of imbalance of the ice masses to present and recent climate, with most of HMA losing ice mass throughout the last century, so losses under RCP4.5 are 73%, and under G4 59%, of present-day glacier mass. Iceland has been closer to balance until recently.

210
The relative effectiveness of SAI on reducing surface runoff over VIC is much smaller than for the Greenland Ice Sheet (Yue et al., 2021), largely due to the compensating impact of AMOC changes. How will this change in future when non-SMB components are also considered? Greenland ice sheet MB has been simulated by PISM under forcing from 4 ESMs simulating the RCP2.6, and RCP8.5 scenarios (Goelzer et al., 2021). AMOC declines during 2020-2089 irrespective of ESMs and scenarios (Fig. 8a), associated with weakened deep water convection over the northern North Atlantic, but clearly, ocean 215 overturning near Iceland is more vigorous than Greenland (Fig. 8d, e). AMOC under G4 is significantly stronger than that under RCP4.5, producing about 10 W m -2 higher heat flux from ocean to atmosphere around Iceland than over Greenland (Yue et al., 2021). Fig. 8b-c shows that VIC MB is very significantly dependent on AMOC (R=0.91, p<0.01), while for Greenland its impact is not significant (R=0.42, p=0.35), consistent with the SMB behavior (Yue et al., 2021). Thus, AMOC plays a role through its impact on SMB and indirectly on dynamic change through the long-term changes in ice cap geometry. Because 220 VIC is much thinner than the Greenland ice sheet, and has higher accumulation and ablation rates, the mass turnover time in VIC is at least 10 times faster than in Greenland meaning that surface climate may induce larger dynamic effects earlier.   Instead, we find that non-SMB accounts for 0.25 m yr -1 mass loss for all future scenarios. We find only an extra 0-2% volume and area loss by 2089 because of the SMB-elevation feedback on VIC geometry, consistent with the 1-3% found by Schmidt et al. (2019) who considered feedbacks based on temperature and precipitation lapse rates. The SMB-elevation feedback will become increasingly important in longer simulations e.g., VIC volume and area will additionally reduce by 9-14% and 9-20% by 2300 than without feedback corrections (Schmidt et al., 2020). For a period as short as the 50-year SAI application specified 235 https: //doi.org/10.5194/tc-2021-318 Preprint. Discussion started: 19 October 2021 c Author(s) 2021. CC BY 4.0 License.
in G4, the change in ice elevation feedback and the induced ice dynamical effects might be expected to be rather too small to be seen in large ice caps with a mass turnover time measured in hundreds or thousands of years. However, over Vestfonna ice cap (Svalbard), scenario-dependent impacts from 10 to 50% (RCP2.6 to RCP8.5) appeared on century timescales (Schä fer et al., 2015), and even over the large Greenland ice sheet, effects due to changing elevation-SMB and ice dynamics become detectable (5-8%) by 2100 (Edwards et al., 2014;Goetlzer et al., 2020) and amount to 14-31% after 200 years (Goetlzer et al., 240 2013). Moreover, the extreme maritime environment of VIC would make it one of the glaciers most likely to exhibit a dynamical response to the SAI or RCP scenario, but we see no such effect. Perhaps the chief reason for the scenarioinsensitivity is that the proportion of calving in non-SMB is small; if the calving losses at the terminus are converted to an equivalent area-average ablation rate over the upper surface of the ice cap, they only amount to 0.06 m yr -1 (26%) over the recent few decades (Jóhannesson et al., 2020). The remaining 74% come from basal melt by geothermal heating, volcanic 245 eruptions, and dissipation of potential energy. Furthermore, retreat of the margins from the ocean limits ocean losses as the result of increased calving into rapidly growing terminus lake in Jökulsá rlón (location see Fig. 1). Hence it is perhaps not surprising we find no trend in dynamic losses in future scenarios.
(2019), (17% volume loss for RCP4.5), perhaps unsurprisingly as we use the same ice dynamic model but with different SMB forcing, which leads to local differences in various basin ice thicknesses by 2089. 255 The influence of volcanic deposits on snowpack albedo has not been incorporated into any snowpack simulations as yet, especially the relatively parameterized SEMIC model; in some (infrequent) years volcanic ash blankets some parts of the ice cap, drastically changing local ablation rates. The steep geometry of some outlet glaciers is still not perfectly captured by the 0.025°×0.025° grid although the bias-correction using satellite observations of albedo compensates for the resolution in the 260 higher resolution model grid.
Of the four ESMs we utilize here, the two MIROC models share many of the same model parameterizations except for the Although geoengineering by SAI is not particularly effective for VIC, and probably other Icelandic ice caps, it does still slow 270 the rate of ice loss. Given the unique geographical location of the ice cap, we may infer that SAI as specified by G4 will not lead to greater mass loss of any glacier of ice cap in the northern hemisphere than expected under any plausible greenhouse gas scenario. However, we are not advocating that SAI be done. Arctic SAI would inevitably modify the entire global climate.
If the population of Iceland wished to make local efforts to conserve their ice caps, locally targeted interventions may offer more palatable governance issues (Moore et al., 2020). 275

Conclusions
We simulate the VIC volume and area change during the period 1982-2089 using the ice sheet model PISM that forced with monthly SMB fields calculated by four ESMs. By 2089, the SAI we simulate reduces VIC mass loss by 4 percentage points, demonstrating that even in Iceland where the impact of AMOC is most felt, SAI could help preserve VIC from melting. The relative unimportance of calving losses due to lack of ocean terminating glaciers, the high basal melt due to active geothermal 280 areas, and the compensating changes in temperature and accumulation due to AMOC mean that VIC is relatively insensitive to climate scenario, with the ice dynamics especially so. We find that ice dynamics are almost constant over both time and scenario because they are relatively unaffected by changing air and ocean temperatures.

285
Code and data availability. All scripts used for simulations are available upon request from the authors. The modeled ice cap volume and area data by PISM are available at https://doi.org/10.5281/zenodo.5410852 Competing interests. The authors declare that they have no conflict of interest.
Acknowledgements. We thank Yimeng Xu for assistance with the draft writing. We thank the climate modeling groups for participating in the Geoengineering Model Intercomparison Project and the scientists managing the earth system grid data 295 nodes who have assisted with making GeoMIP output available.