Resolving glacial isostatic adjustment (GIA) in response to modern and future ice loss at marine grounding lines in West Antarctica

Abstract. Accurate glacial isostatic adjustment (GIA) modelling in the cryosphere is required for interpreting satellite, geophysical and geological records
and for assessing the feedbacks of Earth deformation and sea-level change on marine ice-sheet grounding lines. GIA modelling in areas of active ice loss
in West Antarctica is particularly challenging because the ice is underlain by laterally varying mantle viscosities that are up to several orders of
magnitude lower than the global average, leading to a faster and more localised response of the solid Earth to ongoing and future ice-sheet retreat
and necessitating GIA models that incorporate 3-D viscoelastic Earth structure. Improvements to GIA models allow for computation of the viscoelastic
response of the Earth to surface ice loading at sub-kilometre resolution, and ice-sheet models and observational products now provide the inputs to
GIA models at comparably unprecedented detail. However, the resolution required to accurately capture GIA in models remains poorly understood, and
high-resolution calculations come at heavy computational expense. We adopt a 3-D GIA model with a range of Earth structure models based on recent
seismic tomography and geodetic data to perform a comprehensive analysis of the influence of grid resolution on predictions of GIA in the Amundsen
Sea Embayment (ASE) in West Antarctica. Through idealised sensitivity testing down to sub-kilometre resolution with spatially isolated ice loading
changes, we find that a grid resolution of ∼ 13 of the radius of the load or higher is required to accurately capture the elastic
response of the Earth. However, when we consider more realistic, spatially coherent ice loss scenarios based on modern observational records and
future ice-sheet model projections and adopt a viscoelastic Earth, we find that predicted deformation and sea-level change along the grounding line
converge to within 5 % with grid resolutions of 7.5 km or higher, and to within 2 % for grid resolutions of 3.75 km and
higher, even when the input ice model is on a 1 km grid. Furthermore, we show that low mantle viscosities beneath the ASE lead to viscous
deformation that contributes to the instrumental record on decadal timescales and equals or dominates over elastic effects by the end of the 21st
century. Our findings suggest that for the range of resolutions of 1.9–15 km that we considered, the error due to adopting a coarser grid
in this region is negligible compared to the effect of neglecting viscous effects and the uncertainty in the adopted mantle viscosity structure.



S1 GIA modelling setup considerations
Here we briefly explore the influence of model setup factors that impact the predicted GIA, namely: input load resolution, GIA model resolution, and loading changes outside the region of interest. We use the ICE-RD model to explore these issues, which provides ice thickness at 10 km across Antarctica with a region at 1 km resolution 5 over the ASE. From this, we produce an ice model on the GIA model grid in a variety of ways summarised in Table A1. ANT_10km is the 10 km AIS-wide ice sheet model run and ASE_1km is the nested 1km ice sheet model run conducted under the same model forcing and receiving boundary conditions from the continental run. Figure S3.1 shows the ice models ASE_1km and ANT_10km, and corresponding modelled sea level change due to a purely elastic GIA response across the 150-years of ice loss from 1950 to 2100. Figure S3.2 provides an 10 overview of the errors due to different GIA model setup methods (Fig. S3.2d a, b,c) relative to the error from GIA model resolution (Fig. S3.2d).

15
Since ASE_1km only covers our region of interest (ROI) the Amundsen Sea Embayment, we first explore the importance of adding the loading pattern outside the ROI using ASE_ANT and ASE_1km ( Fig. S3.2a).
Comparing simulations with fixed and evolving ice outside the ASE indicate that deformation due to mass changes outside domain of interest result in a broad, superimposed signal of uplift or subsidence. Not 20 considering mass changes outside the region of interest (e.g. Kachuck et al., 2020) can result in a difference in predicted deformation of at least 6% (and up to 50% at the ROI edge) of the overall signal in the region. The implication of this result is that when modelling regional GIA, we must input the surrounding load changes beyond the ROI. The exact bounding region required is outside the scope of this study.

25
To isolate the effect of using different ice sheet model resolutions, we compared the results of ASE_10km and ASE_1km ( Fig. S3.2b). We also explored the effect of using the same load at different resolutions, by resampling the ASE_1km load grid by a factor of 0.2 to result in a 5 km resolution load grid (ASE_5km). We compute the effect of instantaneous removal of the ice load change from 1950 to 2100 and find that calculations of the resulting elastic GIA response are influenced by: ASE_1km -ASE_5km) we find up to 14 cm difference in SL predictions, with the largest error along the load edge (i.e. grounding line).
• GIA Model Grid Resolution (Fig. S3.2d): For the ice model ASE_1km, improving the GIA model resolution from 7.5 km to 1.9 km produces SL predictions with up to 16 cm difference.

40
Results from Figure S3.2 indicate that refining the GIA model grid resolution from 7.5 to 1.9 km has a similar effect as refining of the input ice load resolution. The effect of load and GIA model resolution both have a predictable pattern whereby the largest error occurs along the load edges. Accordingly, we recommend efforts in improvements in GIA model accuracy go towards constraining the wavelength of ice cover changes, and 45 improving the resolution of the ice model (i.e. ice sheet observations and models) accordingly. The load set up, including interpolation techniques and consideration of the load outside the ROI are also important.

50
This section aims to provide further detail on the scaling factor used to set viscosity variations in our 3-D GIA model and how our scaling factor values compare with results in past literature. We direct readers to the references herein for a complete understanding of the background theory and methodology. As discussed in the main text, 3-D mantle viscosity variability in our Earth models is derived from seismic velocity data by scaling shear wave velocity anomalies to viscosity variability using the method described in Ivins and Sammis (1995)  Adopting tabulated Earth properties of activation energy and background temperature . ! for the whole mantle convection from Kaufmann at el, 2005, it is possible to extend ! into an equivalent scaling function with the surface value of 0.02071 K -1 , slowly decreasing from the depth of 300 km towards the CMB to take the value of 0.01625 K -1 . A hybrid approach is then to write where F(r) is a normalized function of depth ranging from 1 (first 300 km) to 0.785 (CMB) and ! #$%& is a "free"     (Table S1). Each frame represents the difference in GIA predictions due to a) out of region ice loading, b) ice sheet model resolution, c) load resolution and d) GIA model resolution.