Resolving GIA in response to modern and future ice loss at marine grounding lines in West Antarctica

Accurate glacial isostatic adjustment (GIA) modeling in the cryosphere is required for interpreting satellite, 10 geophysical and geological records and to assess the feedbacks of Earth deformation and sea level change on marine ice-sheet grounding lines. Assessing GIA 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 localized 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 15 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 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 20 (ASE) in West Antarctica. Through idealized sensitivity testing down to sub-kilometer resolution with spatially isolated ice loading changes, we find that a grid resolution of ~3 times the radius of the load 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 errors of less than 5% along the grounding line can be achieved with a 7.5 km grid, and less than 2% with a 3.75 km grid, even when the input ice 25 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. 30 https://doi.org/10.5194/tc-2021-232 Preprint. Discussion started: 11 August 2021 c © Author(s) 2021. CC BY 4.0 License.

geophysical and geological records and to assess the feedbacks of Earth deformation and sea level change on marine ice-sheet grounding lines. Assessing GIA 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 localized 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 15 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 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 20 (ASE) in West Antarctica. Through idealized sensitivity testing down to sub-kilometer resolution with spatially isolated ice loading changes, we find that a grid resolution of ~3 times the radius of the load 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 errors of less than 5% along the grounding line can be achieved with a 7.5 km grid, and less than 2% with a 3.75 km grid, even when the input ice 25 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 21 st 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. 30

Introduction
Changes in sea level in response to ice-mass loss are spatially variable because of glacial isostatic adjustment (GIA), which is the deformational, gravitational, and rotational response of the viscoelastic solid Earth to changes in surface ice and water distribution. The response of the bedrock consists of an instantaneous elastic response of the solid planet's lithosphere and 35 mantle, and a longer timescale viscous relaxation of the mantle towards isostatic equilibrium. GIA models produce predictions of changes in the height of the sea surface equipotential and solid Earth surface (i.e. sea-level changes) in response to surface ice loading changes, which are in turn used to interpret satellite, geophysical and geological records and serve as input to models of ice-sheet dynamics.

40
Accurate GIA modelling is required to constrain the sea level and solid earth feedbacks on ice dynamics in the coming centuries, especially along unstable marine-grounded ice fronts in Antarctica where bedrock uplift and gravitational drawdown of the sea surface due to ice loss act as a negative feedback to stabilise the retreat of marine-grounded ice-sheet grounding lines (e.g. Gomez et al., 2010;De Boer et al., 2014;Konrad et al., 2015;Larour et al., 2019). Furthermore, the GIA response to past and modern ice cover changes makes a significant contribution to satellite records of modern mass 45 changes in marine sectors of the West Antarctica that are actively experiencing ice loss (e.g. King et al., 2012;the IMBIE team, 2018).
GIA modelling in Antarctica is complicated by the fact that the continent is characterised by strong lateral variability in lithospheric thickness and upper mantle viscosity, with low viscosities in the west and high viscosities in the east (e.g. An et 50 al., 2015a;Heeszel et al., 2016;Lloyd et al., 2020). In particular, the low-viscosity mantle and thin lithosphere observed under the West Antarctic Ice Sheet (WAIS) identified from increasingly resolved seismic tomography and geodetic and geologic constraints (Ritzwoller et al., 2001;Morelli and Danesi, 2004;Kaufmann et al., 2005;Nield et al., 2014;Heeszel et al., 2016;Barletta et al., 2018;Shen et al., 2018;Lloyd et al., 2019) leads to a more spatially localised (short wavelength) and faster viscoelastic response to surface loading than cratonic regions covered by Late Pleistocene ice sheets (e.g. Hay et al., 2017;55 Powell et al., 2020). Over West Antarctica, upper mantle viscosities are thought to vary by several orders of magnitude over short spatial scales reaching as low as 10 18 Pa s in the Amundsen Sea Embayment (ASE) beneath areas of active marine ice loss (e.g. Nield et al., 2014;Barletta et al., 2018). This implies that viscous effects due to 20 th century and more recent ice loss will become significant on annual to decadal timescales and accelerate during the timeframe of instrumental records (Barletta et al., 2018;Powell et al., 2020). Viscous effects due to ongoing ice loss also have the potential to influence ice sheet grounding 60 lines in the coming centuries (Gomez et al., 2015) but have not been included in recent high resolution coupled projections (Larour et al., 2019).
To accurately capture the timing and wavelength of GIA effects across Antarctica, models must be capable of both accounting for 3-D earth structure (i.e., 3-D GIA models such as Latychev et al., 2005or van der Wal et al., 2015, and be of sufficient 65 spatiotemporal resolution to capture the geometry of grounded ice cover. Commonly, GIA, ice-sheet and coupled ice-sheet -GIA modelling (e.g. Gomez et al., 2015;Konrad et al. 2015) studies of the sea level change in response to modern and future ice loss have been performed with only 1-D (radially varying) Earth structure models (e.g., Kendall et al., 2005;Spada and Stocchi, 2007;Adhikhari et al., 2016), or with coarse spatial resolutions of > 20 km due to the computational expense. GIA models capable of km-to sub-km-resolution have also been developed (e.g. the 1-D GIA model by Adhikhari et al., 2016;the 70 3-D GIA model by Latychev, et al. 2005 with updates described in supplementary materials of Gomez et al., 2018). For computational efficiency, some of these models implement regional grid refinement techniques which allow for a higher resolution along ice retreat margins. Alongside this, improvements in ice-sheet models (e.g. Nowicki et al., 2016;Seroussi et al., 2020;DeConto et al., 2021) and observational products (e.g. Studinger, 2014;Bamber et al., 2020;Smith et al., 2020;Morlinghem et al., 2020) allow similarly high-resolution km-to sub-km-ice thickness and bedrock topography datasets that 75 serve as input to GIA models. These advancements allow for GIA models to capture short-wavelength bedrock deformation and input ice loading changes at unprecedented detail, but at heavy computational expense, particularly for global 3-D GIA models.
It is well-established that dynamic ice-sheet models are sensitive to the chosen grid resolution (e.g. Durand et al., 2009;Van 80 den Berg & Van de Wal, 2006), requiring km to sub-km resolution to accurately represent ice dynamics and grounding line migration in some applications (e.g. Gladstone et al., 2012;Pattyn et al., 2013;Cornford et al., 2016). It has also been suggested that as fine as 1-km resolution bedrock topography may be needed to capture the influence of fine-scale topographic features on the ice-sheet evolution (Durand et al., 2009) and high resolution may also be needed to represent some embayment walls and pinning points that act to slow down retreat (e.g. Favier et al., 2012;Joughin et al., 2014;Berger et al., 2016). 85 While topographic features themselves can be very fine scale, the changes in bedrock elevation and sea level in response to ice cover changes tend to be longer wavelength, and the corresponding spatial resolution required to accurately resolve these changes in GIA models and their influence on ice dynamics remains poorly understood. Larour et al. (2019) suggested that kilometer-scale resolution may be required to capture the elastic component of deformation in response to ice loss. However, 90 the idealized tests they performed considered an isolated, and increasingly localized load, and their conclusion may not hold for more realistic, spatially coherent ice loss scenarios. Furthermore, their model did not include viscous deformation in response to ongoing ice loss during the simulation, or account for lateral variations in Earth structure. There have been limited studies assessing the length scale of realistic viscoelastic bedrock response beneath the structurally complex WAIS, though a recent high-resolution bedrock deformation modelling study by Zwinger et al. (2020) suggests a convergence in modelled 95 deformation at resolutions finer than 5 km. The broad spatial nature of bedrock deformation and the spatially coherent nature https://doi.org/10.5194/tc-2021-232 Preprint. Discussion started: 11 August 2021 c Author(s) 2021. CC BY 4.0 License. of ice-sheet retreat, which becomes less localised over longer timescales, suggest that sub-km to km grid resolution, which comes at great computational cost, may not be necessary for accurate GIA model calculations.
The aim of this study is to assess the sensitivity of predictions of GIA in response to modern and future ice loss to spatial 100 resolution in the rheologically complex, marine sectors of the WAIS. We build a 3-D viscosity model based on the most recent Antarctic-wide seismic tomography model (Lloyd et al., 2019) to serve as input to a 3-D finite volume, global sea level model (Latychev et al., 2005) to assess the performance of 3-D GIA model predictions across surface grid resolutions of 1.9-15 km.
We repeat calculations with a range of Earth models, considering the contribution from elastic and viscous deformation separately. We focus on the response to observed modern ice loss over the last two decades  and 105 projected future ice sheet retreat in the coming century DeConto et al., 2021) in the Amundsen Sea Embayment of West Antarctica. Our study is motivated by the following questions: What 3-D grid resolution is necessary to adequately capture the elastic and viscous uplift and associated gravitational and rotational effects in response ice loading changes? How significant is the effect of grid resolution compared to sources of uncertainty and simplifications made in some previous modelling, in particular the neglect of viscous deformation. 110

Methods
To investigate the influence of GIA model grid resolution, we first conduct idealized load sensitivity tests over a range of surface grid resolutions from 7.5 to 0.5 km, for the instantaneous removal of cylindrical loads from 0.5 to 16 km in radius (Section 3). We then widen our "aperture" to assess the model grid resolution required to accurately capture GIA due to modern 115 ice-sheet cover changes from satellite observations  and future ice loss from ice-sheet model projections DeConto et al., 2021) in the rapidly evolving Amundsen Sea Embayment of West Antarctica (Section 4). We choose to locate our study region (light blue square in Fig. 1a) on the ASE both because of the ongoing ice loss and vulnerability to large-scale future marine ice-sheet retreat in the region, and because the region is characterized by low upper mantle viscosities and a thinned lithosphere (e.g. Barletta et al., 2018) making the ice there sensitive to solid-earth ice-sheet 120 feedbacks. In the ASE, the Pine Island and Thwaites Glaciers together contributed 95 Gt/year of the 159 ± 8 Gt/year total WAIS mass flux in 2017 (Rignot et al., 2019), with studies estimating that collapse of Thwaites Glacier is already underway (Joughin et al., 2014). Accurate GIA predictions are critical to assess rates of future ice-sheet retreat and associated sea level changes making it an ideal location to study the effects of grid resolution on modelled GIA. To represent the radially and laterally variable Earth rheology in this region, we use a viscoelastic Earth rheology and a range of 3-D viscosity structure 125 models in Antarctica derived from seismic tomography (An et al., 2015a;Heeszel et al., 2016;Lloyd et al., 2019). We adopt a range of 3-D Earth model grids with surface resolutions from 1.9-15 km and compare results to first assess the resolution required to capture the elastic deformation associated with ice loss. Lastly, we repeat these experiments with viscoelastic https://doi.org/10.5194/tc-2021-232 Preprint. Discussion started: 11 August 2021 c Author(s) 2021. CC BY 4.0 License. effects and compare results to the elastic calculations to assess the contribution of viscous effects to modern and future sea level, and the model resolution required to capture these effects. In the sections that follow, we describe the adopted 3-D GIA 130 model, computational grids, Earth rheological models, and adopted modern and future ice loss scenarios.

3-D GIA Model
We compute GIA predictions with a global 3-D finite volume sea level and Earth deformation model (Latychev et al., 2005) capable of regional grid refinement (Gomez et al., 2018). The model solves the sea level equation (Kendall et al., 2005) over a 3-D global spherical tetrahedral grid defined from the surface to the core-mantle boundary (CMB) that allows us to resolve 135 the laterally and radially varying Earth structure. which are a strong feature in Antarctica. We also adopt this 3-D GIA model because it is capable of regional grid refinement to achieve regional resolution at sub-km scale in regions of interest within a lower resolution globe (Section 2.2). The model computes gravitationally self-consistent solutions for the sea level equation, incorporating effects of time-varying shorelines, Earth rotation changes and viscoelastic deformation of the Earth assuming an elastically compressible Maxwell viscoelastic rheology. The GIA model requires two main inputs, a 3-D Earth model of 140 viscoelastic rheological properties and a time series of ice thickness changes. These components are described in the following sections. The model also requires global topography as an initial boundary condition, including the elevation of the bedrock beneath the ice. Topography globally outside of Antarctica is set to etopo2 (NOAA National Geophysical Data Centre, 2006) in all experiments with a realistic loading scenario, and the Antarctic bedrock elevation for each of the experiments is described below. 145

Computational model grid: regional grid refinement
To compute GIA model predictions, we construct model grids of various surface resolutions ( Fig. 1) with the regional grid refinement process detailed in the supplementary material of Gomez et al. (2018). Grid refinement is achieved by incrementally bisecting grid edges in the selected region to achieve a desired resolution, and a final smoothing operation along the region boundary to ensure a well-behaved transition. We perform calculations on a base grid with a global surface resolution of 15-150 km, which consists of 20 million nodes and 70 radial layers between the core-mantle boundary and the Earth's surface. The radial layers are defined to respect the unconformities in the radially varying (1-D) seismic reference model STW105 (Kustowski et al., 2008), with the shallowest layers at 12, 25 and 43 km depth. Eight regionally refined grids are constructed from this base grid: five for the idealized load sensitivity tests at 7.5, 3.75, 1.75, 1 to 0.5 km surface grid resolution over a minimum 0.3-degree radius region around the test load, and three for the more realistic calculations from observed modern 155 and future model projected ice loss, at 7.5, 3.75 and 1.9 km surface resolution in incrementally smaller regions converging over the ASE (Fig. 1). The highest-resolution 1.9 km grid over the ASE has ~ 29 million nodes. As our study focusses on surface grid resolution, the grid refinement is limited to the surface few layers down to ~10 km. Testing with deeper grid refinement during experiment design process indicated that this was sufficient, and our results indicate that km-scale model resolution is only needed at the surface to accurately capture the geometry of surface loading. 160

Earth rheological model
The spatial pattern and amplitude of surface deformation in response to ice loading is dependent on the underlying Earth structure. For the idealized sensitivity tests in section 3, we adopt a purely elastic Earth model with a 1-D elastic and density structure. In section 4, we move to a set of realistic simulations using observed or modelled AIS ice loading, adopting 3-D 165 viscoelastic Earth models with a range of viscosity structures constrained by seismic tomography (An et al., 2015a,b;Lloyd et al., 2019) and informed by GNSS-inferences of local mantle structure. The elastic and density structure for these models are based on seismic reference model STW105 (Kustowski et al., 2008). Laterally varying lithospheric thickness ( Fig. S1(d)) in all simulations is a composite of a regional lithospheric thickness model by An et al. (2015a) over Antarctica, and a global lithospheric thickness model by Conrad and Lithgow-Bertelloni (2006) everywhere else. Over Antarctica, lithospheric 170 thickness is scaled to have an average of 96 km, resulting in a minimum of 40 km, as was done in Hay et al. (2017).
3-D mantle viscosity variability is estimated from seismic velocity models by scaling isotropic seismic shear wave velocity anomalies to a viscosity variability term using the method developed by Ivins and Sammis (1995). This term defines the variability of mantle viscosity with reference to a chosen 1-D viscosity profile. We follow the same procedure as described in 175 Latychev et al. (2005), Austermann et al. (2013) and Gomez et al. (2018), whereby lateral variations in mantle viscosity are established by sequentially converting the field of relative variations in isotropic shear wave velocity to density, temperature and eventually viscosity anomalies. The final conversion adopts a factor that scales the dependence of viscosity to temperature variations, determining the peak-to-peak lateral variation in viscosity. A different scaling factor is applied to regional and global seismic velocity models to account for the amplitude differences between the different seismic velocity models. 180 To address the substantial uncertainty in Earth structure, we repeat our simulations with three different 3-D viscosity models: EM1_L (Figs. 1c, d, S1), EM1_M, and EM2. EM1_M and EM1_L adopt the latest high-resolution 25 km Antarctic seismic tomography model by Lloyd et al. (2019)   in the region south of 45°S and extending from the surface down to the transition zone, and S362ANI (Kustowski et al., 2008) for the rest of the globe. The two variations EM1_M and EM1_L were 185 scaled to represent a moderate range of viscosities across Antarctica that match regional averages, and a lower-viscosity endmember adjusted to match GPS-derived inferences of minimum viscosity beneath WAIS (Nield et al., 2014;Zhao et al., 2017;Barletta et al., 2018).
The viscosity variations in EM1_M are more moderate with scaling factors of 0.0263 for ANT-20 and 0.035 for S362ANI,190 both close to the preferred value in Kaufmann et al. (2005). These viscosity variations are superimposed on a 1-D viscosity profile of 5 x 10 20 Pa s in the upper mantle and 5 x 10 21 Pa s in the lower mantle typical in most GIA-based inferences of mantle viscosity (e.g., Mitrovica and Forte, 2004). As noted, the viscosity variations in EM1_L (Figs. 1c, d, S1) are of higher amplitude https://doi.org/10.5194/tc-2021-232 Preprint. Discussion started: 11 August 2021 c Author(s) 2021. CC BY 4.0 License. with lower viscosities under WAIS and the ASE. EM1_L adopts larger scaling factors of 0.033 for ANT-20 and 0.04 for S362ANI and an accompanying 1-D viscosity profile of 1 x 10 20 Pa s in the upper mantle, which is aligned with Lambeck et 195 al., (2014) estimates of 1-D upper mantle viscosity using far-field sea-level proxy records, and 5 x 10 21 Pa s in the lower mantle.
The larger scaling factors applied in EM1_L were calibrated to best reflect the absolute upper mantle viscosity estimates from dynamically derived GPS bedrock uplift rates in the WAIS at three locations: ~ 6 x 10 17 to 2 x 10 18 Pa s at the northern Antarctic Peninsula (Nield et al., 2014), ~2 x 10 19 to 2 x 10 20 Pa s at the Fleming Glacier in central Antarctic Peninsula (Zhao et al., 2017), and ~2.5 x 10 18 to 4 x 10 19 Pa s at the Amundsen Sea Embayment (Barletta et al., 2018). Figures 1c, d and S1 shows 200 the resulting low viscosity earth model structure (EM1_L), which has the lowest viscosity at Marie Byrd land of ~ 9 x 10 17 to 7 x 10 18 Pa s in the upper mantle.

in the East Antarctica and Antarctic 205
Peninsula, and the model by Heeszel et al (2016) for West and Central Antarctica. The full construction of EM2 is detailed in Hay et al. (2017).

Ice model and topography
We consider three ice melt scenarios with resolution ranging from 1-to 5 km in the ASE. The total ice thickness change from start to end of each scenario is shown in  For observations of modern ice loss, we adoted a reconstruction we term ICE-SH from Shepherd et al. (2019) of surface elevation change (∆h) from 25 years (1992-2017) of multi-mission satellite altimetry data resolved over a 5-km grid at 5-year intervals. We treat ∆h as a proxy for ice thickness change (Carrivick at al., 2019), apply the Bedmap2 (Fretwell et al., 2013) grounded ice mask and saturate ice thickness change > 20 m/yr to control against spurious data points. Initial ice thickness is 215 given by Bedmap2. Observations of ice surface elevation changes in Antarctica are continuously improving in resolution, and currently range from metre-scale resolution over short observational tracks (e.g. Studinger et al., 2014), to sub-km to km-scale resolution at the regional scale (e.g. Bamber et al., 2020), to ~ 5 to 35 km from radar and laser satellite altimetry derived records over the whole Antarctic (e.g. Martin Español et al., 2016;Schröder et al., 2019;Shepherd et al., 2019;Smith et al., 2020). ICE-SH was selected from the available observational datasets due to the 5-km resolution being the highest of its class 220 of available satellite altimetry derived-records providing decadal time span surface elevation change records covering the whole Antarctic.
For Antarctic evolution over the next century, we apply modelled ice thickness changes from two Antarctic-wide ice-sheet model projections: (1) ICE-GOL which predicts AIS evolution under RCP 8.5 and including meltwater feedbacks from 2000 225 to 2100 at 5-km resolution over 5-year intervals , and (2) ICE-RD which predicts AIS evolution from 1950 to 2100 at 10-km resolution over the whole AIS with a nested 1-km resolution simulation over ASE at annual resolution (Extended Data Fig. 5 from DeConto et al., 2021). For ICE-RD, we ran simulations at yearly intervals from 1950-2100 but the interval from 1950-2000 is a period of ice model initialization and we therefore focus on the period between 2000 and 2100 in our results. We also take initial Antarctic bedrock topography from the ice models. Further information of each model can be 230 found in the corresponding references. In selecting these scenarios, the goal is to provide a representative sample of spatially coherent ice-sheet retreat scenarios at high resolution from the literature, rather than to capture all possible projected future ice loss scenarios.

Idealized experiments: sensitivity of elastic uplift predictions to grid resolution
Our main goal in this analysis is to assess the relationship between grid resolution and elastic GIA predictions, and to identify, 235 for a given load dimension, the grid resolution required to accurately model the associated GIA response. Realistic ice retreat has complex geometry, making it difficult to pinpoint the cause of GIA inaccuracies due to resolution, which may be due to poor representation of the ice load, or numerical errors in representation of the response to Earth loading. To isolate the effect of changing grid resolution on GIA predictions, we first perform a suite of idealized sensitivity tests modelling the instantaneous elastic deformation from unloading of an isolated cylindrical ice load with differing surface grid resolutions that 240 are iteratively bisected from 7.5 km down to 0.5 km. We chose to perform the test with short wavelength, spatially isolated ice loading changes because these would be most poorly represented by a coarse grid compared to coherent ice loss over a broader area. Furthermore, these tests with idealized loads are less computationally costly and enable a systematic assessment reaching higher spatial resolution. In total, 85 GIA simulations were run using 17 ice cylinders of height 100 m and radii ranging from 0.5 to 16 km (0.5, 1, 2,…,15, 16), across five different grid surface resolutions: 0.5, 1, 1.75, 3.75 and 7.5 km. 245 The simulations are performed with the purely elastic 1-D earth model (Sect. 2.3), and an idealized topography of 3800 m south of 24.5 °S and -835 m everywhere else to reflect the 30:70 land to sea ratio on earth. A radially symmetric ice sheet with steady-state Antarctic ice dome profile (Paterson and Colbeck, 1980) sits on top of this topography extending from the south pole to 69 o S, with a maximum height of 3500 m. We also consistently place the centre of the cylindrical load on an 250 arbitrary model grid node in the ASE (76°S 150°W).

Idealized experiment results
Figure 3 summarises the error in predicted elastic GIA response with varying load radius and grid resolution. Figure 3a shows the bedrock deformation along a transect from the centre of the load for ice cylinders of 2, 5 and 10 km radius, with maximum bedrock uplift predicted on the finest 0.5 km grid of 48, 120 and 185 mm, respectively. Figure 3b indicates whether the grid 255 over-(blue) or under-(red) represents the mass of the load within the model. Errors in the GIA prediction compared to the result for the finest resolution 0.5 km grid (yellow lines) are typically the largest at the load centre where they underestimate the magnitude of peak displacement. Although a coarser grid may either under-or over-estimate the mass of loading represented in our model (Fig. 3b), it will always dampen the magnitude (effective height) of the load by spreading the load area over a larger grid region. For example, a 5 km radius load will be represented by 3 grid points on a 7.5 km grid, resulting 260 in an overestimated 11.25 km radius loading footprint. For certain radius and grid combinations, the wider load footprint on a coarser grid leads to another zone of peak error occurring outside the load edge (e.g., compare dashed black line to yellow line in Fig. 3a). Even further from the load, the magnitude of deformation decreases and the results from various grid resolutions begin to converge.

265
These sensitivity test results highlight that the accuracy of predictions depends on the placement of the edge of the load relative to grid nodes and find that the load will be best captured if its edge lies sufficiently close to a grid node (e.g., in Fig. 3b the 1.75 km grid more closely captures the 2 km radius load than the 3 km radius load). Finally, the grid is unable to resolve the load when the grid resolution is more than approximately three times the radius of the load. This is illustrated, e.g., in the black solid line in Fig. 3a, where unloading a 2 km radius load on a 7.5 km grid resulted in no deformation, whereas a 3 km radius 270 load is captured with a 7.5 km grid.
To quantify grid-related error, we calculate the difference between a given simulation and the corresponding simulation with the finest (0.5 km) resolution. We plot the root mean square error (RMSE; Fig. 3c) as an absolute measure of error and the average and standard deviation of the absolute percentage error (Figs. 3d, e) as a relative measure of error, beneath and within 275 2-km of the loaded region. Fig. 3c shows that the magnitude of RMSE remains relatively constant. This RMSE remained between ~10-20 mm for a 3.75 km grid, and ~ 20-40 mm for a 7.5 km grid, for example. As the load radius increases, the magnitude of deformation increases as well. However, the error due to grid resolution becomes less significant relative to the total deformation (i.e. the percentage error in Figs. 3d, e decreases).

280
While the dependence of grid performance on load position relative to grid nodes complicates matters, in order to arrive at an approximate relationship between grid resolution and load size, we assume a linear relationship between the two, which allows us to estimate, for this GIA model, a threshold beyond which grid-related error becomes sufficiently low to no longer merit a further refinement in grid resolution (Figs. 3d, e). For example, in the cross-sectional view of deformation in Fig. 3a, the 10 km radius load deformation is equally well represented by grid resolutions between 0.5 to 1.75 km. Considering the average 285 absolute percentage error (Figs. 3d, e), we found that a 3:1 ratio (represented visually on Figs. 3b-e in the form of a black dashed line) between load radius and grid resolution (6:1 ratio between load diameter and grid resolution) brings the error to < 7 ±3 % (where 3% represents one standard deviation of the absolute percentage error calculated within 2 km of the load region). These results provide a rough estimate of the magnitude of error one can expect from a given model resolution and loading scenario, and can serve as a guide for selecting the appropriate grid resolution for a given problem. For example, for an input load with significant isolated locations of ice loss ~ 3 km in radius (or a ~ 6 km in diameter), a grid resolution of 1 km should be adopted. However, these idealized cylindrical load experiments are unlikely to capture the sensitivity of GIA predictions to grid resolution when realistic ice loss geometries are adopted. Such geometries are rarely characterized by spatially localized 295 loads. Furthermore, these experiments capture only elastic deformation and neglect viscous effects, which can be significant on short timescales in low viscosity zones of the West Antarctic. In the following sections we explore how the dependence on grid resolution of GIA model predictions identified here changes when more realistic ice loss geometries and 3-D viscoelastic Earth structure are adopted.

Results with realistic modern and future ice loss in the Amundsen Sea 300
In this section, we consider the importance of grid resolution error for more realistic, spatially coherent modern and future ice loss scenarios. We begin with a consideration of the influence of grid resolution on sea level change in simulations adopting a purely elastic Earth model in section 4.1. Following this, we adopt 3-D viscoelastic Earth models to consider the contribution to sea level change from viscous deformation.  (Fig. 2d), while for future ice loss projections, the sea level fall reaches up to 9.06 m and up to 12.8 m from 2000 -2100 for ICE-GOL and ICE-RD respectively. Note that in addition to signal coming from local ice loss in the ASE, ice outside this region of interest also contributes a broad signal of smaller magnitude (see Supplementary Section S1). 315

Influence of grid resolution on predictions of elastic deformation 305
To explore the resolution dependence of sea level predictions that adopt an elastic Earth model, we repeat the calculations in Figs. 2d-f with a surface grid resolution ranging between 1.9 and 15 km. Fig. 4a-i, shows the difference between results for simulations performed at 1.9 km grid resolution relative to coarser resolutions. The coarser grid simulations tend to underestimate the magnitude of ice unloading and associated sea level fall in most of the domain (red regions in Fig. 4). 320 The highest grid resolution error occurs at the periphery of the load within a few kilometres of the final grounding line position rather than at the location of maximum deformation (compare Figs. 3d-f to Fig.4) suggesting that the representation of the load is the main source of error rather than the representation of the response of the Earth to the loading. For example, for the ICE-GOL ice loss scenario, the greatest difference between 1.9 and 15 km grid simulations ranges occurs along the entire final 325 grounding line (Fig. 4b), but the maximum sea level fall of over 9 meters occurs ~ 2 km away from the grounding line.
The error decreases with increasing resolution, with minimal differences between the 1.9 km and 3.75 km grid resolutions.
The maximum absolute error in the case of a 15 km grid (i.e. the maximum difference between the 15 km and 1.9 km resolution cases) is 44 cm at 2100 in ICE-RD (Fig. 4c), 47 cm at 2100 in ICE-GOL (Fig. 4b) and 9.1 cm after 25 years of modern melt 330 in ICE-SH (Fig. 4a). That is 3.4%, 5.2% and 13% of the peak elastic sea level fall predicted at that time for each respective scenario. The errors are approximately an order of magnitude smaller when a grid resolution of 3.75 km is adopted: 9 cm for ICE-RD, 5 cm for ICE-GOL and 0.7 cm for ICE-SH, or 0.7%, 0.4% and 1.0% of the peak elastic sea level fall respectively.
Since maximum grid resolution error is concentrated along the grounding line for elastic runs, in Fig. 5 we explore how the 335 error evolves during the ICE-RD simulation along a 10-km region bounding the grounding line. The error increases in absolute magnitude with increasing ice loss ("Error" in Fig. 5a; Fig. 5b) but the relative error decreases across the same runs ("Percentage Error" in Fig. 5a). In the case of 15 / 7.5 / 3.75 km grid resolutions, the peak error is ~10 / 5 / 1 cm at the 25 year mark of the simulations, and ~50 / 15 / 5 cm at the 150 year mark (whiskers in Fig. 5a top row). In contrast, the percentage, peaks at 20 / 6 / 1.5% of the signal at 25 years and drops to < 5 / < 2 / < 0.3% after 150 years (Fig. 5a). This decrease in 340 percent error with time reflects that the ice geometry changes become broader in wavelength and can therefore be resolved by a coarser grid compared to the more spatially isolated changes occurring earlier in the simulation. Given the uncertainty in modelled and observed ice loss and bedrock elevation in Antarctica (e.g. Morlinghem et al., 2020), we suggest that for most applications, errors of less than 5% can be achieved with a 7.5 km grid, and errors of less than 2% with a 3.75 km grid.

Contribution of viscous earth deformation to sea level predictions 345
So far we have focussed on the resolution dependence of the contribution to sea level change from elastic Earth deformation, as this has been a focus of recent literature (Larour et al., 2019). However, the Antarctic Ice Sheet is underlain by strongly laterally varying viscosity structure, and the Amundsen Sea region in particular is underlain by a low viscosity zone and thinned lithosphere (e.g. Barletta et al., 2018;Lloyd et al., 2020). Viscous deformation associated with ongoing ice loss is neglected in Larour et al. (2019) but is expected to be significant on decadal to centennial timescales in this region. In Figs. 2g-i the 350 calculations of sea level change associated with the three ice loss scenarios shown in Figs. 2d-f are repeated with the 3-D viscoelastic Earth EM1_L described earlier. As with the elastic case, sea level falls beneath regions that experience ice loss in all three cases, but the magnitude of the sea level fall is significantly larger than predictions based on an elastic Earth model (compare Figs. 2g-i to Figs. 2d-f). In particular, peak sea level fall in this case reaches -0.79 m over 25 years in the ICE-SH ice loss scenario, and -14.9 m and -29.1 m from 2000 to 2100 for ice loss scenarios ICE-GOL and ICE-RD, respectively. The 355 latter (Fig. 2i) is more than double the sea level calculated with the elastic Earth model (Fig. 2f). viscous deformation to occur (Figs. 6a-c), but as we highlight below, substantial viscous deformation still occurs along the current grounding line in the simulation. This is less evident in the modern because migration of the location of maximum ice mass loss is minimal. Faint blue areas further from the region of ice retreat in Fig. 6b indicate a sea level rise due to peripheral bulge subsidence, a viscous process.

Influence of resolution on viscoelastic sea level predictions 370
In Fig. 7, we repeat the assessment of grid resolution error in Fig. 4, but with a viscoelastic rheology based on the 3-D earth model EM1_L. With the inclusion of viscous behaviour, the magnitude of the grid resolution error is similar to the elastic case (compare the range of errors on Figs. 7 and 4) but the spatial pattern of the error becomes more complex. The maximum error is no longer solely concentrated along the current grounding line since the solid earth continues to respond viscously to the poorly resolved loading changes along previous locations of the grounding line. This is particularly evident in the ICE-RD 375 simulations where the grounding line retreats across a large area. In this case, the grid resolution error over the region of past ice loss and grounding line migration is equal to or larger than the error along the active grounding line (Figs. 7b,c). The error increases during the simulation as viscous deformation builds, and thus it also has a dependence on the 3D viscosity structure.
Note that in the blue region of Fig. 7c, the sign of the error is not the same as in Fig. 7f because ice retreat is not consistent 380 within this particular region in ICE-RD. Specifically, between the years ~2020 to 2050 in ICE-RD, the grounding line in this blue region stays relatively fixed, experiencing multiple episodes of localized ice retreat and re-advance (unloading and loading). Situations of re-advance tend to occur at lateral scales < 15 km, such that these sub-grid scale movements were not adequately captured with a coarser, 15 km grid. Along the final grounding line for ICE-RD, the coarser grids consistently result in a lower magnitude sea level fall as in the elastic case.

Time evolution of influence of grid resolution and viscous effects
To compare the relative contributions of grid resolution and Earth model differences over time, we extract predicted sea level 400 time series from all simulations with the ICE-RD ice loading scenario and elastic and viscoelastic Earth models at two locations in Fig. 8:  lines) because the upper mantle viscosity is set to 10 20 Pa s and thus a significantly longer timescale is required for viscous effects to become significant. Differences between these two simulations and any of the simulations adopting 3-D Earth structure is larger. For example, the differences between simulations using EM1_L and the 1-D viscoelastic earth model, reach up to 16.5 m and 12.7 m at the sites shown in Figs. 8b and c, respectively, by 2100.

410
Starting early in the century, the influence of grid resolution becomes smaller than the effect of adopting different Earth models (compare the differences between the dashed and solid lines to differences between different colored lines in Fig. 8). Using the 1-D average earth model (Fig. 8 red lines), viscous effects start to emerge after 50 years and reach only 4% of the peak signal by the end of the run at year 2100. Nevertheless this signal is more significant than the error incurred by using 15 km grid resolution versus of 1.9 km grid resolution by 2040. With a 3-D earth rheology and low viscosities beneath the ASE, 415 viscous effects are pronounced within decades in the simulation (blue and green lines in Fig. 8) and become larger than the difference in predictions based on 15 km and 1.9 km grid resolutions within 25-30 years and before substantial ice loss has occurred in the simulation. Figure 9 provides a more detailed picture along the grounding line of the contributions to differences in predicted sea level 420 described in the preceding sections. We consider the impact of each factor on the predicted sea level signal at the end of the simulations, plotting the distribution of differences between simulations across all grid points within 10-km of the final grounding line. In interpreting Fig. 9, note once again that in viscoelastic runs, the region of maximum grid-resolution error does not necessarily occur near the grounding line (e.g. Figs. 7c, f). To visualise the distribution, we plot a classic box-whisker diagram where the boxes represent (from left to right) the 25 th percentile, median and 75 th percentile of the distribution, whilst 425 the whiskers represent the "minimum" and "maximum" (25 th and 75 th percentiles minus 1.5 times the interquartile range), overlain by a density curve. A box-whisker plot was chosen to mitigate against the effect of outlier points, which we plot as hollow diamonds.
Note that the error due to grid resolution consistently has a unimodal distribution that peaks at ~0-1% near the grounding line 430 across a range of durations from 25-years of ICE-SH, to 100 and 150 years (starting from 1950) of ICE-RD. The hollow diamonds show that a significant number of points are statistical outliers, which is likely due to the fact that predicted sealevel change is low in magnitude at some regions along the grounding line that experience less ice-loss, causing even a small magnitude of error to contribute a large percentage error. Inclusive of these outliers, the range of percent error due to a 15-km grid peaks at ~20% after 25-years of the ICE-SH simulation, but decreases to less than ± 8% by 100-years of the ICE-RD 435 simulation as the magnitude of predicted sea level fall becomes larger across the entire region and the spread of outlier points diminishes. We conclude that the range denoted by the box-whisker plot likely provides a more accurate assessment of the error contributed by each factor in Fig. 9 in zones of active grounding line migration. However, in the following paragraphs we continue to describe the range of errors inclusive of outlier points so as to not under-estimate the possible spread.

440
Our results indicate that the differences due to the choice of adopted Earth model equals and, in most cases, exceeds the size of the error due to grid resolution near the grounding line by the end of all our simulations. Over the 25-year modern observed ice loss scenario, the difference in predictions associated with earth model configuration lies between ~2-10% within 10-km of the grounding line, which is within the range of error due to insufficient grid resolution in a viscoelastic run, which ranges from ± ~20% with a 15 km grid, ± 6% with a 7.5 km grid, and < 3% with a 3.75 km grid at all grid points (range of the 445 diamonds outside the shaded distributions in Fig. 9a, though noting the discussion above, the percent error is substantially smaller than these end member values at most points). However, with more ice loss over longer timescales, the difference due to adopted Earth model far exceeds the grid resolution error (Figs 9b, c).
If we look beyond the grounding line and consider the difference in predicted sea level between different adopted 3-D viscosity 450 models in the entire study region (difference in results with EM1_L and EM1_M, plotted in Fig. 6 d-f at 2100), the lower viscosity model results in additional viscous deformation that is up to 10.2 %, 21.4% and 20.9 % of the total signal after 25years of ICE-SH, 100-years of ICE-RD and 150-years of ICE-RD respectively. In all cases, the error due to neglecting viscous effects far surpasses the error due to grid resolution with a 15-km grid (compare the bottom rows to top 3 rows of Figs 9b and c)., with the lower viscosity EM1_L model resulting in an additional viscous deformation that is up to 23.8 %, 58.9 %, 62.4 % 455 of the elastic signal near the grounding line after 25-years of ICE-SH, and 100-years of ICE-RD and 150-years of ICE-RD respectively.

Discussion
Our study provides an assessment of the model grid resolution needed to capture decadal to centennial-scale GIA in the vicinity of active ice loss. We targeted the ASE in West Antarctica as our study location as it is a region with ongoing and projected 460 marine ice sheet retreat and where low mantle viscosity and thin lithosphere result in a rapid and localised GIA response to ice loading. We adopted a 3-D GIA model to accurately capture the viscoelastic response at high resolution, including the complexity introduced by laterally varying earth rheology in the region. Accurate assessments of solid earth deformation from past and present ice evolution are important for constraining the negative sea level -solid earth feedback on ice sheet retreat, and more accurate interpretation of geophysical observables. For the former, our study focusses on the sensitivity of GIA 465 predictions along the ice sheet grounding line where this feedback occurs.

Influence of grid resolution
For our suite of simulations with elastic and viscoelastic Earth models, modern and 21 st century ice loss scenarios, and surface grid resolution ranging between 15 and 1.9 km, we found that improvements in the accuracy of model predictions with increasing grid resolution was limited, remaining within centimetres to decimetres at most at the grounding line. Furthermore, 470 our results converged at higher resolutions, with errors from a 3.75 km grid resolution reaching at most 6 cm within 10 km of the grounding line in all simulations, even when the input ice sheet model results were available at 1 km resolution in the case of ICE-RD. The error introduced in assuming an elastic Earth model and neglecting viscous deformation in the ASE builds to an order of magnitude or more larger than the grid resolution error within three to four decades, and up to tens of meters by the end of the century. In addition, predictions adopting different 3-D Earth models that reflecting the uncertainty in 475 viscoelastic Earth structure in the region diverge by up to 1 meter within 50 years and upwards of 2-3 meters after 100 years in the simulation.
For coupled ice-sheet -GIA model applications, our results suggest that adopting high resolution in the ice sheet model does not require a similarly high-resolution GIA model. In our simulations, a 3.75 km grid was sufficient to bring errors to < 2% 480 along the grounding line for all scenarios (Fig. 9). Furthermore, this percentage decreased over time our simulations, and would continue to decrease multi-century and millennial simulations as the magnitude of viscous deformation and the scale of the ice loss continue to grow. While bedrock topography has smaller scale features (Morlinghem et al., 2020), our results suggest that the GIA signal is less localized and may be computed at lower resolution relative to the ice sheet dynamics and then interpolated and added to the initial topography on the higher resolution ice sheet model grid, as is done in, e.g., Gomez et al. 485 (2015) and DeConto et al. (2021).
Our results showing that the location of maximum error consistently lies along the load edge for elastic model runs (Fig. 4ac) suggest that the error due to coarse model resolution is predominantly a result of poor representation of surface ice cover changes rather than representation of the solid Earth response. For the latter, we would expect the error to occur instead at the 490 location of maximum GIA response (compare differing spatial patterns in Fig. 2e, f to 4b, c). When the viscous response is incorporated, the time-evolving nature of viscous deformation leads to an additional peak in grid resolution error at locations past the grounding line due to inaccurate representation of past loading (Fig. 7). This additional zone of error will not affect active ice sheet grounding lines (though it may be important for interpretation of modern records), and while the spatial pattern of the error differs, the magnitude of the error due to grid resolution was similar across both elastic and viscoelastic simulations. 495 Our findings on the size and source of resolution error are in contrast to recent work by Larour et al., (2019) who suggested that kilometre resolution was required to capture elastic deformation. This discrepancy may be due in part to Larour et al. (2019) considering only point loads in their idealized resolution experiment, while our conclusions are based on more realistic, spatially coherent ice loss scenarios. Differences may also arise due to the nature of the computational grid and processing of 500 inputs (see section 5.3) which should be explored in more detail in future GIA model inter comparison efforts. Nevertheless, our predictions based on an elastic Earth model converge to theirs for more spatially broad loads.
One possible limitation in this study is we do not reach sub-km grid resolution in our GIA model, and our highest resolution ice model is 1 km. In sensitivity tests with idealized loading scenarios in Section 3 we adopted a grid resolution as low as 0.5 505 km grid and found that a minimum 1:3 ratio between grid resolution and load radius was required for results to incur an error of < 5±3 (s) % along the load edge, suggesting that a 3.75km grid would be unable to capture a spatially isolated, < 1 km radius ice unloading event. That we did not see significant error at this resolution in the realistic simulations indicates that the ice cover changes are spatially coherent and there are no significant spatially isolated ice unloading events (i.e. no ice thickness changes occurring over only a few grid points) predicted in the 1-km resolution ICE-RD ice model simulation (Fig. 5). 510 To further investigate if short-wavelength, spatially isolated ice loss scenarios exist over Antarctica, we assessed the surface elevation change observables from 40 and 25 years of multi-mission satellite altimetry data by Schröder et al., (2019) and Shepherd et al. (2019) respectively, and 15 years of airborne laser altimetry from Operation IceBridge (OIB ATM L4; Studinger, 2014(Updated 2018). While a more detailed investigation is merited, in our initial analysis of these datasets we 515 observed that spatially isolated ice loss events have a lower magnitude, only persist over short timescales, and found no evidence of high magnitude, short-wavelength ice loss occurring with spatial scales < 5 km. Thus, we expect that spatially isolated ice unloading occurs rarely and will not have a significant impact on the overall accuracy of GIA model results in a given region. Nonetheless, with improving observational products and ice sheet model resolutions, we expect to obtain regional-scale ice loading grids of sub-km resolution that may warrant further study with a sub-km GIA model grid (e.g. Durkin 520 et al., 2020).

Influence of viscous deformation and earth model uncertainty
Within decades in the ASE, viscous deformation is a significant contributor to the GIA signal (Hay et al., 2017;Barletta et al., 2018;Powell et al., 2020;Kachuck et al., 2020). The GIA response can be decomposed into the following: perturbation of the sea surface equipotential, elastic deformation and viscous deformation of the solid earth surface. Previous studies have isolated 525 the GIA signals from each of these components to assess the relative importance of each factor on ice sheet dynamics. Over decadal to centennial timescales, Larour et al. (2019) show that purely elastic deformation was more significant than the sea surface perturbation on continental-scales, while Kachuck et al. (2020) found that viscous deformation is more significant than either elastic deformation or the perturbation to the sea surface. In this study, we have confirmed that viscous deformation effects are significant within decades, particularly in the low viscosity region of the ASE, where the viscous component of 530 deformation can reach multi-metre scales by the end of the century. This body of work highlights the importance of incorporating viscous behaviour in GIA modelling applications in regions of low mantle viscosity.
Complicating efforts to accurately characterise the viscous deformation behaviour is the uncertainty in Earth's viscosity structure. The timescale of viscous solid-earth deformation on ice sheet dynamics is strongly dependent on the assumed earth 535 rheology. The global average mantle viscosity of ~ 10 21 Pa s (Forte and Mitrovica, 1996) corresponds to response times from centuries to millennia, whilst recent seismic (Lloyd et al., 2020) and GPS observations suggest an upper mantle viscosity under the ASE as low as ~ 10 18 Pa s (Barletta et al., 2018). Rapid viscous uplift response was similarly identified in Kachuck et al. (2020) which used a 2-D GIA model of mantle viscoelastic deformation and found that sea level fall associated with viscoelastic mantle deformation led to a 30% reduction in modelled ice sheet volume loss by 2150. Our study compares results 540 generated with three 3-D earth rheology models; EM1_M and EM2 have a comparable viscosity range, while EM1_L has the lowest viscosity values under the ASE. These results demonstrate that uncertainties regarding the effective mantle viscosity are significant and can contribute up to multiple metres of uncertainty by the end of a 100-year simulation (Fig. 6). Furthermore, additional uncertainty arises from the model of viscoelastic behaviour. We adopt a viscous Maxwell rheology, but studies suggest that incorporating a short-term transient component of deformation may result in even faster viscous deformation (e.g. 545 Pollitz, 2019).

On GIA model setup
In choosing a method for representing a finer resolution load grid onto a coarser model grid, we found that it is important to consider how the model itself discretizes the load, and the input load interpolation schemes. Here, it is worth noting that our GIA model grid is a tetrahedral grid (triangular grid on surface), and these findings may not translate perfectly to other model 550 grid compositions. Our GIA model grid consists of a uniform global tetrahedral grid that allows for regional patches of refined resolution (also uniform) but does not permit matching of model grid nodes to the input grid. For our experiments, by comparing the volumes of the input ice calculated on the input and GIA model grids, we found that the in-built Poisson interpolation scheme (Latychev et al., 2005) performed better in interpolating the finest resolution load grid onto the model grid compared to other tested schemes, suggesting that an understanding of the method in which the load in mapped onto the 555 model grid nodes is important. Additionally, we note that considerations such as the resolution of the input ice sheet model and treatment of the ice cover outside the region of interest also have an influence on the final GIA model predictions (see Supplementary Section S1) and should be explored further in future studies.

Conclusion
In this study, we present a comprehensive analysis of the influence of grid resolution on modelled GIA effects in response to 560 ice cover changes over the modern satellite era and through the 21 st century. We adopt a range of Earth models including models that capture lateral variations in Earth structure based on seismic tomography and GNSS analyses. These experiments showed that: (1) the grid resolution error introduced through adopting a 15 km grid relative to a 1.9 km model grid remains within centimetres to decimetres throughout our simulations; (2) the grid resolution error is the highest along load edges for purely elastic deformation cases, and along past and current grounding line positions for viscoelastic Earth models, and is 565 primarily associated with the representation of the surface load; (3) results with grid refinement beyond 3.75 km converged in our simulations, even when adopting a 1 km resolution input load, and this likely represents a conservative lower bound since the next coarser grid we considered was 7.5 km. The errors associated with the choice of grid resolution will decrease with time for longer simulations as the extent and magnitude of ice loss and associated GIA response increase. Furthermore, comparison of simulations adopting elastic and 3-D viscoelastic Earth models demonstrate that the contribution of viscous 570 deformation can be up to tens of metres over the 21 st century, or > 50% of the total deformation signal. Additionally, uncertainties in earth properties can contribute up to several metres of error. This indicates the importance of considering viscous uplift when modelling GIA over decadal to centennial timescales in the ASE. In comparison, the error due to grid resolution is negligible for grids of spacing of 3.75 km and less.

575
To supplement these findings with realistic loading, we conducted a sensitivity test with cylindrical loads with radii from 16 km to 0.5 km and grid resolutions from 7.5 to 0.5 km. These indicate a minimum 1:3 ratio between the required grid resolution and the load radius (or 1:6 with load diameter) to minimise grid resolution error. However, no significant spatially isolated loads occur in our adopted observation-and model-based ice loss scenarios, and a preliminary examination of other ice observation and modelled products suggest that significant ice loss with < 5 km wavelength is rare in the ASE. These results, 580 taken together, support the conclusion that km-scale resolution in GIA modelling is generally not necessary. However, as higher resolution sub-km ice observational and dynamic ice model grid products are released, this guidance may have to be revisited.

Code/Data Availability
We will make all model output from the sensitivity tests and more realistic simulations available on a public repository. The 585 3-D GIA model adopted here has been used in numerous previous studies, questions regarding the model or requests for additional output can be discussed with the corresponding author and K.L., the developer of the code. Additional data related to this paper may be requested from the authors.     Figure 5. Evolution of error in elastic sea level predictions due to grid resolution from 1950 to 2100 with the 1km input resolution ICE-RD ice model. (a) Box-whisker plots of the error and percent error (see methods) calculated from the difference in predicted sea level changes from the start of the simulation to the indicated time within 10 km of the grounding line at that time between a simulation with 1.9 km resolution and simulations adopting 15 km (light pink), 7.5 km (medium pink) and 3.75 km (dark pink) grid resolutions. The box represents (from left to right) the 25th percentile, median and 75th percentile of the distribution, whilst the 805 whiskers represent the "minimum" (25th percentile -1.5 x the interquartile range) and "maximum (75th percentile -1.5 x the interquartile range). Error (m) is the difference between sea-level predictions from the higher -lower resolution run. Percentage Error (%) is calculated as 100* (SL1.9km -SLlowres)/SL1.9km for each grid point.

810
Frames (a-c) shows the difference in sea level change predicted from simulations adopting 3-D viscoelastic earth model EM1_L and an elastic Earth model. Frames (d-f) shows the difference in sea level change predicted from simulations adopting two different 3-D viscoelastic earth models EM1_L and EM1_M. Note the difference in scale. Time frames and ice models are as indicated at the top of the columns. In each frame, the black and blue line indicates final and initial grounding lines respectively for each simulation and is annotated with the maximum and minimum data value within the frame. All simulations adopt a 1.9 km grid resolution.

835
we plot a classic box-whisker diagram overlain with a density curve. The box represents (from left to right) the 25 th percentile, median and 75 th percentile of the distribution, whilst the whiskers represent the "minimum" (25 th percentile -1.5 x the interquartile range) and "maximum" (75 th percentile -1.5 x the interquartile range). The diamonds outside the whiskers represent outliers.