Glacial Isostatic Adjustment (GIA) in Greenland: a Review

This review provides updated estimates of the glacial isostatic adjustment (GIA) component of present-day uplift at a suite of Global Navigation Satellite System (GNSS) sites in Greenland using the most recently published global ice sheet deglaciation histories. For some areas of Greenland (e.g. the north-west and north-east), the use of GNSS to estimate elastic uplift rates resulting from contemporary mass balance changes is more affected by the choice of GIA correction applied compared to other regions (e.g. central-west). The contribution of GIA to GRACE estimates of mass imbalance is becoming increasingly insignificant for large areas of Greenland as it enters a period of extreme warmth, and in total represents <5 % contribution (−6 to +10 Gt/year) to the observed Greenland-wide mass trends over the last decade. However, differences between deglacial histories and uncertainties in their assumed viscoelastic Earth structure combine to result in significantly different region-by-region estimates of GIA.


Introduction
The surface of the solid Earth is continually adjusting and responding to external (e.g. atmospheric loading, tidal loading) and internal (mantle flow) forces exerted upon it. Whilst many of the short-term elastic readjustments are tangible (e.g. tectonic plate friction and resultant earthquakes), the Earth is still trying to reach isostatic equilibrium in response to deglaciation of the Pleistocene Ice Sheets that occupied a significant proportion of the northern hemisphere and the advance and retreat cycles of the Antarctic and Greenland ice sheets since the Last Glacial Maximum (LGM). The process of ongoing viscoelastic relaxation in response to this redistribution of (specifically) ice (glacio-isostasy) and water (hydro-isostasy) on the Earth's surface is termed 'glacial isostatic adjustment' (henceforth 'GIA'). Regions located both inside and outside former ice sheet centres are still responding to the deglaciation of many of the northern hemisphere ice complexes (Laurentide, Cordilleran, Innuitian, Eurasian and the British-Irish) that concluded several thousand years ago.
The process of GIA produces four observable geodetic perturbations on the region it acts upon: (1) vertical and (2) horizontal motion (deformation) of the Earth's surface; (3) downward deflection of the geoid in previously glaciated regions with respect to the geoid of a uniform, topographically invariant Earth; and (4) a linear contribution to the secular rate of geoid height change. As GIA is a result of retreat/regrowth cycles of grounded ice, vertical (and horizontal) rates of deformation are quantities that may be inverted to provide constraints on ice thickness changes over millennial timescales, for an assumed viscoelastic configuration and rheology of the Earth. Horizontal plate motion has, in the past, been utilised to provide extra constraints on deglaciation histories, e.g. [1,2]. However, owing to uncertainties in the horizontal deformation signal associated with choice of GIA model, larger or equivalent in magnitude horizontal deformation signals and associated contaminants such as tectonic and/or rotational strain [3] or deformation resulting from proximal tectonic plates, vertical uplift rates are generally prioritised as GIA constraints. This review paper focuses specifically on the role of GIA in Greenland and aims to bring the reader up to date with the current state of the research regarding quantification of the GIA signal and modelling attempts to constrain GIA and to remove its effect from geodetic datasets such as the Gravity Recovery and Climate Experiment (GRACE) and Global Navigation Satellite System (GNSS) data. Specifically, this paper will 1. Present estimates of present-day uplift patterns by correcting Greenland network (GNET) network GNSS uplift observations with predicted GIA estimates 2. Conclude that the GIA correction due to past changes is small compared to the present elastic signal, and although being a relatively small component of the GRACE mass balance data, it still remains an important quantity to constrain further in specific areas of Greenland 3. Highlight areas of disagreement in the field and suggest areas for future study  (Fig. 1a). The GNET (Greenland GPS network) of GNSS stations is composed of 59 mostly ice-proximal sites ( Fig. 1a) and is an international project led by Ohio State University, the National Space Institute at the Danish Technical University and the University of Luxembourg, and receives technical support by the NSF-supported facilities UNAVCO and Polar Field Services. Station coordinates for GNET are available from https://www.unavco.org. Regional GNSS networks used to monitor crustal deformation tend to get denser over time, and this has certainly been true in Greenland, though GNET is limited to the peripheral portions of Greenland where bedrock is exposed. A smaller subset of GNSS stations was set up in south-west Greenland in 1995 and remained until 2002 as part of a 'campaign-type' study by [5]. Coordinates of ten stations in southwest Greenland were monitored at five intervals between 1995 and 2002 at the same time each year to minimise seasonal variations on uplift rate from, e.g. ice and air mass changes, later shown by [4] to contribute significantly to surface deformation depending on the time of year. We review briefly the findings of [5] but keep in mind the epoch-specific nature of the measurements made in this study mean that the uplift rates measured between 1995 and 2002 are not directly comparable to the results presented in Fig. 1 and, as such, should be treated separately and are presented in Fig. 2. The contribution of GIA to the uplift rates of [4] and [5] is reviewed in BModelling of the Glacial Isostatic Adjustment Component of Vertical Deformation in Greenland^section. Figure 1 presents uplift data from a network of 59 stations around the periphery of the Greenland Ice Sheet. The basic daily analysis of station coordinates was performed using MIT's GAMIT and GLOBK software [6,7]. The GNET stations were processed along with hundreds of global stations, as described in [8]. These velocities are referred to the geometrical reference frame associated with VREF and HREF, whereas the GIA predictions presented here utilise a centre of mass (CM) reference frame. The differences in vertical GPS velocities expressed in this geometrical frame, and those expressed in a CM frame, are of order of 1 mm/year in the vertical [9].
The black circles in Fig. 1b-e show that all GNSS sites have been uplifting since their individual deployment dates. Some, but not all, of largest mean uplift rates are found at GNSS sites HEL2 (Fig. 1e), KAGA (Fig. 1c) and KUAQ (Fig. 1e) situated close to the grounding lines of major outlet glaciers (Helheim (HH), Jakobshavn Isbrae (JI) and Kangerdlugssuaq (KL), respectively, see Fig. 1a, inset). These fast-flowing glaciers, along with the Petermann Glacier (PG; Fig. 1a, inset), are sourced by drainage basins that comprise approximately one fifth of the area of the Greenland Ice Sheet [10].
Over the time span of the HEL2 data  [11]. For Helheim and Kangerdluqssuaq,the dynamic component of mass loss increased sharply during 2005, facilitating a shift to a thinning regime for both glaciers in response to a period of warm ocean and air temperatures [12]. Uplift rates quoted for HEL2 display significant accelerations Differentials of ∼5 mm/year between ice margin and coastal sites (e.g. DGJD/VFDG and SCOR/SCOB, Fig. 1d) highlight the localised nature of the elastic response to present-day ice unloading, also demonstrated in [13][14][15][16]. The north-east of Greenland displays the lowest average uplift rate (7 mm/year) and least spatial variability in uplift rates, given similar time spans and number of observations (Fig. 1d). SCBY station, located approximately 80 km upstream of the current calving front of Petermann glacier (PG; Fig. 1a, inset) on its western flank, has experienced a moderate deceleration in uplift of 0.62 ± 0.31 mm/ year 2 between 2008 (12 mm/year) and 2013 (8.5 mm/year). Prior to this period of deceleration, [16] predict an increase in uplift rate in this overall region of from 0 to 5 mm/year in 2004-2007 to 5-10 mm/year by 2006-2009, the latter spanning a period in which the nearby Humboldt glacier (HG, Fig. 1a, inset) lost 1-5 km of floating ice from its northern flank [17]. For this region, [11] show that the dynamic component of mass balance contributed positively to the regional ice mass budget to the tune of approximately 10 [11].
Between stations DKSG and UPVK, where the area of icefree land is minimal (Fig. 1a), uplift rates vary by a factor of 2 ( Fig. 1b). All stations along this transect have station data for equivalent time periods, so the variability in this area results from proximity of the individual GNSS station to the large number of marine-terminating, fast-flowing outlet glaciers in the area. Several outlet glaciers in this region have demonstrated increases in velocity of 20-60 % from 2005 to 2010 with a regional mean increase in velocity of 18 % [18,19], whilst the Upernavik Isström (near UPVK) underwent a rapid retreat and thinning phase during 2005-2009 dominated by dynamic ice loss [20]. Moreover, a regional analysis carried out by [21] demonstrated that of the 46.6 Gt of ice loss in the north-west drainage basin during 2005-2010, 55 % originated from dynamic ice loss concentrated at the margins.
In the west of Greenland, the current GNET network has fewer stations compared to other regions (Fig. 1c). One of the longest-standing GNSS stations in Greenland, KELY, records a mean uplift for the time period 1995.6 to 2013 of 2.  (1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002), the campaign dataset from west Greenland presented in [5] (Fig. 2) displays surface deformation in the range of −3.4 to 1.6 mm/year, with ice-proximal sites showing the greatest degree of subsidence, interpreted to be largely a result of late Holocene advance. This is in contrast to the trends observed by [4], which displayed no subsidence at similar sites (KELY, KAPI, NUUK) in western Greenland, reflecting a change in mass balance regime in the intermittent period between measurements. Indeed, this is demonstrated in basin-scale analysis of mass balance by [22] showing a 50 % increase in discharge in the Nuuk fjord drainage basin from 7.4 to 11.1Gt/year from 1996 to 2007. On a smaller scale, [23] detected decreasing SMB in the Nuuk catchment area from 2002 to 2010 accompanied by an increase in discharge.
In BModelling of the Glacial Isostatic Adjustment Component of Vertical Deformation in Greenland^section, we produce estimates of the GIA component of the uplift rates described in this section using two recent deglacial histories: Huy3 [24] and ICE-6G_C (VM5a) [25].

Tide Gauges
Tide gauges record temporal changes in height of the water column, resulting from deflection of the solid Earth surface and changes in the height of the sea surface (which Fig. 2 Vertical velocity (a, b) map of GNSS sites contained within the west Greenland campaign network of [5] for 1995-2002 (black circles). In this case, the error bars represent the root-mean-square error. Colour conventions for symbols as in Fig. 1 approximates to the geoid). As such, sea-level data from tide gauges have been analysed in a number of studies in order to estimate the component of sea-level change resulting from solid Earth deformation prior to the installation of expansive GNSS networks and also as constraints on Earth viscosity profiles [26][27][28]. In Greenland, the Permanent Service for Mean Sea Level (PSMSL: www.psmsl.org) has archived data from nine tide gauges, situated mainly in the south of Greenland. No tide gauge datasets in Greenland are categorised by PSMSL as Revised Local Reference ('RLR') datasets (i.e. they have a well-documented local tide gauge benchmark and are quality controlled by PSMSL) and should therefore be used with caution if analysing for secular trends. The longest dataset in Greenland is found in Nuuk (near NUUK; Fig. 1a), but the tide gauge is now decommissioned. The University of Hawaii Sea Level Centre (UHSLC) provides research-quality tide gauge data from three sites for 2008-2014 (Thule, Scoresby Sund) and 2006-2014 (Qaqortoq). A trend analysis of the Nuuk tide gauge was performed by [29], who showed that sea level rose at Nuuk at a rate of 1.93 mm/year between 1958 and 2002, and sites from the wider network were used in combination with other local pressure gauge measurements by [30] to decompose the tidal signal in western Greenland into its harmonic constituents. To the authors' knowledge, aside from [29], data from the current tide gauge network have not been used for any other geodetic analyses, likely due to uncertainties surrounding benchmarks and the lack of intra-network temporal continuity of records.

Modelling of the Glacial Isostatic Adjustment Component of Vertical Deformation in Greenland
The most common method to estimate GIA for any region is to have a constrained a priori estimate of the global deglaciation history ideally extending as far back to the Last Interglacial (LIG) e.g. [14]. Typically, this global ice evolution dataset is then prescribed as primary input into a sea-level/ solid Earth response model [31,32] to provide a forward model estimate of present-day uplift rates for a given viscoelastic configuration of a compressible one-dimensional, spherically symmetric gravitationally self-consistent rotating Earth. However, newer studies have revisited the assumptions of linear [33] and composite [34,35] rheology and demonstrated that rheological assumptions, along with simplistic 1D depth-varying seismic Earth structure and velocity profile [36], provide measureable differences between predictions of surface deformation, and we refer the reader to the aforementioned publications for more detail. To the best of the authors' knowledge, GIA studies related specifically to the Greenland Ice Sheet have not applied non-Maxwellian rheologies to predict solid Earth response. In summary, present-day rates of GIA produced by a priori models are primarily a result of (1) the deglacial history and (2) the assumed Earth configuration, i.e. the elastic (lithospheric thickness, L) and viscous (upper, ν um , and lower mantle viscosity, ν lm ) profile and rheological characteristics of the Earth's subsurface. The latter is usually determined by errors in the deglacial history or often used as a tuning parameter to fit constraining relative sea-level (RSL) data.
In this section, we briefly review the GIA signal associated with deglacial histories by produced by Fleming and Lambeck in 2004 ('GREEN1': [37]) and Simpson et al. in 2009 ('Huy2' [38]) where possible by placing them into context with the model output shown in this paper but focus detailed review of the most recent iterations of the ICE-6G_C (VM5a) [25] and Huy3 models [24]. The Huy3 deglacial history is constrained by several hundred observations, e.g. ice extent (terminal moraines, trimlines) and RSL data, that define the deglacial history of Greenland from LGM to present. Huy3 is therefore considered to be one of the more robust representations of Greenland deglaciation since the LGM. In the Huy3 model, the non-Greenland Ice Sheet component of the global deglacial history is derived from the ICE-5G model [39] with alternate North American ice complex reconstructions from a Bayesian-calibrated model [40] to quantify the role of non-Greenland ice on near-field Greenland predictions. The optimal Earth configuration of Huy3 provides predictions of GIA using a Maxwell viscoelastic Earth composed of three layers-a lithosphere with a thickness of 120 km, ν um = 0.5 × 10 21 Pa s extending to 670-km depth with ν lm = 2 × 10 21 Pa s. The only difference between this preferred Earth model for Huy3 and its predecessor, Huy2, is that lower mantle viscosity of the latter is a factor of 2 smaller (1 × 10 21 Pa s). The most recent iteration of W. R. Peltier's series of global deglaciation models, ICE-6G_C (VM5a) [25], does not have an updated Greenland component compared to ICE-5G [39] and is derived from [41] and named 'GrB'. The accompanying earth model, VM5a, is a model that more crudely represents the complex VM2 model into five subsurface layers. The uppermost viscoelastic layer is composed of a 60-km-thick elastic lithosphere underlain by a 40-km sublithospheric layer of viscosity 1 × 10 22 Pa s. The remainder of mantle is subdivided into three layers: an upper mantle (ν um ) with viscosity 7 × 10 20 Pa s extending to 660 km and a two-layer lower mantle (660-1250 km; 2 × 10 21 Pa s; 1250-2900 km; 5 × 10 21 Pa s) extending to the core-mantle boundary [42,43].
The Huy3 deglacial history produces a spatial pattern of GIA (Fig. 3a) for Greenland that is broadly not too dissimilar to its predecessor, Huy2, and to an extent GREEN1 [37]. On a Greenland-wide scale, the GIA signal in Greenland is mostly formed by the two distinct contributions from both the evolution of the Greenland Ice Sheet and the deglaciation of the Laurentide Ice Sheet-the latter more dominant in the west of Greenland, along with a smaller contribution from the postdeglaciation ocean loading and associated rotational feedback. A large area of intense subsidence is located in the south-west of Greenland (Fig. 3a), with maximum subsidence rates of 6 mm/year predicted at KAPI (Fig. 1c; the difference between the solid black and solid red circles). This area of intense subsidence is the result of neoglacial advance of the western margin of the ice sheet [44], after reaching a minimum configuration at the Holocene Thermal Maximum (HTM). Differences in the timing, magnitude and maximum areal extent of the readvance exist between the regional deglacial models reviewed in this text. For example, in the GREEN1 model, a 40-km readvance of the ice margin from its post Holocene minimum configuration to its present-day margin in western Greenland is required to reproduce data from isolation basins in western Greenland that predict sea-level rise from ca. 2.5 kyr ago to present [44][45][46][47]. The nominal Earth model used in this instance was composed of a lithosphere of 80-km thickness and upper and lower mantle viscosities of 4 × 10 21 and 10 × 10 21 Pa s, respectively. The readvance of the western margin in the Huy2 [38] model occurs after the model has reached a minimum configuration between 5 and 4 kyr BP. A readvance of 80 ± 20 km occurs between 4 kyr BP and present compared to its predecessor, HUY1 [48], whose minimum extent is reached at 3 ka BP with a subsequent advance of 50 km. The most up-to-date model in the series, Huy3, has a neoglacial regrowth in the region of 20-60 km from its minimum extent, which was reached between 5 and 3.5 kyr BP. The GrB model of [41] produces the largest retreat behind the present-day margin (60-120 km) and reaches a minimum, the earliest of all the models considered at 8-9 kyr BP, with the majority of the readvance completed by 4 kyr BP.
The station at the centre of the subsidence 'bullseye' predicted using Huy3, KAPI, has a GIA contribution from Huy3 of −6 mm/year and from ICE-6G_C (VM5a) of −0.09 mm/ year to present-day uplift rates (Fig. 1c). Applying the correction of Huy3 (ICE-6G_C (VM5a)) revises the elastic contribution to average uplift rates from present-day mass variations from the ice sheet as of 2013 to 15 mm/year (9) at KAPI, 5.8 mm/year (2.2) at KELY and 8.4 mm/year (3.7) at NUUK. Predicted rates of GIA-related uplift at the southern margins of the ice sheet also differ greatly with GIA-related uplift regime predicted by ICE6G_C compared to weak subsidence in Huy3. Corrected rates of present-day uplift for the Fig. 3 Comparison of present rates of GIA-related vertical deformation in Greenland calculated using the deglacial chronology originally presented in [24] (Huy3). (a, b) The ICE-6G_C (VM5a) global deglaciation model and Earth model VM5a [25,70]. White lines on the plots denote the zero contour and, therefore, the boundary between areas of predicted GIA-related subsidence and uplift. Black lines provide contours of uplift at the 1 mm/year interval. GNSS sites contained in the GNET network are marked as grey triangles and the campaign network of [5] denoted as grey circles. Comparison of present-day uplift rates from this figure with rate of change of geoid height in Fig. 4 illustrates the relative contributions of each process (surface deformation and geoid height change) to overall GIA-induced trends in present-day sea level in Greenland south of Greenland differ by ∼5 mm/year at SENU and QAQ1 (Fig. 1c) and NNVN (Fig. 1e). The differences in timing, magnitude and extent of readvance between the two models considered here cause significant differences in predicted presentday uplift rates in western Greenland. This is mainly as a result of the timing and the magnitude of the ice advance rather than differences in the preferred Earth model. A review of threshold lake records by [49] confirms that Huy3 provides the most realistic reconstruction of late Holocene margin dynamics in the west and south-west of Greenland. In the Jakobshavn Isbrae region, Huy3 predicts that the ice margin remained behind present between 7 and 4 kyr ago, which is largely in agreement with threshold lake data in this region. GrB predicts a retreat behind present 2000 years earlier than Huy3 and those recorded in lake sediments. This is not a surprise given the large quantity of geological observations employed to constrain the deglacial history of Huy3 in this region compared to earlier models (GREEN1 and GrB) whose observational constraints are not well spatially and temporally distributed, nor were they available in such a large quantity. Both models are unable to resolve a short-lived retreat between 1.5 and 1 kyr ago in the west. In the south-west, the region between KAPI and QEQE-D (Fig. 1a), threshold lakes indicate a margin retreat outside of their present catchments around 9 kyr ago. GrB predicts retreat past the present-day margin at 8 kyr, whereas the margin in Huy3 retreats behind present at 6 ka. Discrepancies between modelled and observed retreat also exist for both GrB and Huy3 in the southern, eastern and south-east regions of Greenland, where the retreat is mistimed by ±2-4 kyr in some cases.
In the north-west of Greenland (e.g. from AASI to THU2), both models predict that the solid Earth is experiencing a net subsidence (−2 to 0 mm/year), mainly due to the sinking of the peripheral bulge region of the Laurentide Ice Sheet (Fig. 3a). THU2, DKSG, ASKY, KULL and UPVK all lie on, or close to, the predicted location of the −2 mm/year contour (Fig. 3a) and, as a result, estimate present-day elastic uplift rates in the region of 5-20 mm/year

The Contribution of Glacial Isostatic Adjustment to Gravity Recovery and Climate Experiment Glacial Isostatic Adjustment Models as Forward Models
The GRACE mission was set up in 2002 in order to monitor secular trends in the Earth's gravity field, targeting specifically mass changes associated with movement of water on the Earth's surface. The mission has proved highly successful at providing observations of regional-and global-scale variations in the Earth's hydrosphere. GRACE has also significantly contributed to knowledge of the mass changes within the cryosphere and estimated that the combined Antarctic and Greenland mass loss rate between 2003 and 2013 is 460 ± 59 Gt/year [50]. Over the GRACE monitoring period, Greenland Ice Sheet mass loss increased from 75 ± 26 Gt/year during 2002-2004 [51] to 179 ± 25 to 232 ± 13 Gt/year by the end of 2008 [52,53,54]. Since 2008, this ice loss has further accelerated to 280 Gt/year as of 2013 [50,55], thus contributing 0.7 ± 0.15 mm/year to present-day rates of sea-level change (3.2 ± 0.4 mm/year [56]; correct for 1993-2009). Figure 4 shows the contribution of the Huy3 and ICE-6G_C (VM5a) GIA correction to the present-day geoid change observed by GRACE. One important feature to note in these plots is the order-of-magnitude difference compared to uplift in Fig. 3, with the change in geoid height not necessarily 'mirroring' the patterns of GIA. Both Huy3 and ICE-6G_C (VM5a) (Fig. 4a, b, respectively) show geoid lowering in the southern part of Greenland, but with differing magnitudes with ICE-6G_C returning a stronger rotational feedback resulting from polar wander [57]. The Huy3 deglacial history shows a more noticeable geoid height change in the west of Greenland, a result of mantle flow into the region occupied by the Laurentide Ice Sheet and a weaker rotational feedback component in this region. This is a consequence of Huy3 employing an alternate glaciologically self-consistent North American ice component. Both models predict geoid uplift across the north of Greenland in the range of 0 and 0.2 mm/ year As GRACE data is a measure of the total sum of contributions to the secular change of the Earth's gravity field, a number of corrections [58] must be made before the data can be interpreted as resulting from changes in water storage on the Earth's surface, resulting in the ∼15-59 Gt/year uncertainty attached to the estimates of mass imbalance quoted thus far. A comprehensive summary of error sources and uncertainties in the GRACE gravity fields is provided by [59], which includes corrections for (1) GIA, (2) atmospheric loading, (3) harmonic truncation and signal smoothing and (4) mass leakage from outside the geographical area of interest. For this review, we focus on (1). For a detailed discussion on (2)-(4), we refer the reader to [59] and [58].
The magnitude of the GIA correction applied to GRACE data, along with the uncertainty bounds supplied (e.g. −2 ± 21 Gt/year for Greenland in [59]), is dependent on several factors, namely the following: (1) choice of deglacial history; (2) subsurface viscosity profile [59]; and (3) the mechanical characteristic physical model that one applies to simulate the rheology of the subsurface, e.g. linear vs. composite rheology, and the compressibility of the Earth [60].
For example, [58] provided a comparative error analysis of three Greenland GIA models with differing compressibility assumptions (compressible/incompressible mantle) and with/without rotational feedback included and found that for Greenland, the corrections for the models considered (ANU with GREEN1 Greenland component [37] and ICE-5G VM2 with an incompressible/compressible Earth [61]) ranged from −5.3 to +9.5 Gt/year for the ice sheet as a whole. On a regionby-region basis, specifically, the central-west portion of Greenland (between QAAR and PAA1 stations; Fig. 1), the disagreement between three different GIA models considered in [58] was pronounced and amounted to a range of 9 Gt/year between the maximum and minimum corrections (−6.8 Gt/year for ICE-5G VM2 with a compressible Earth; to +2.4 for ANU computed using a range of viscoelastic parameters (50 km ≤ L ≤ 100 km; 2 × 10 20 Pa s ≤ ν um ≤5 × 10 20 Pa s; 0.5 × 10 22 Pa s ≤ ν lm ≤ 2 × 10 22 Pa s)). However, as a proportion of the observed mass loss trends in the west of Greenland (70 Gt/year for 2003(70 Gt/year for -2012, the GIA correction is small (<10 %). According to [58], the largest contribution to observed mass trends from GIA is in the north-east of Greenland (40 % of an observed mass loss of ∼10 Gt/year). In [62], this is also a region where there is a significant disagreement between GIA corrections provided by [38,63] and those by [64] with the GIA correction from [64] producing a residual signal of ice sheet thickening (5cm/yr) during 2003-2011, compared to thinning of 2-5 cm/yr for this time period when applying [38] or [63].
In summary, the GIA correction applied to GRACE data collected for Greenland varies between −6 and 9.5 Gt/year [51,58,59,62] with an average uncertainty of ∼26 Gt/year (mainly a result of uncertainty in Earth structure) in the current published literature. Compared to other errors originating from (2)-(4) and methodological differences between data processing centres, the GIA correction itself, combined with uncertainties attached, comprises <10 % of the observed mass loss signal in Greenland.

Recovering Glacial Isostatic Adjustment from Gravity Recovery and Climate Experiment and Other Geodetic Methods
As a result of the increased remote monitoring of the polar ice sheets, coupled with increasingly sophisticated modelling approaches, it is possible to constrain present-day changes in Greenland Ice Sheet mass balance directly from observations. In light of these advances, some groups have adopted inversion-type approaches to solve for the GIA component of mass change in Greenland. For example, [64] simultaneously inverted solid Earth deformation rates from a global network of GNSS stations and mass changes from ocean bottom pressure data-assimilated model output to calculate the GIA and surface mass components of the GRACE mass signal, producing a GIA component in Greenland of −69 ± 19 Gt/ year as calculated by [62], an order of magnitude larger than corrections derived from mainly forward models. [16] used geodetic observations of mass balance provided by the ICESat as input to a regional elastic rebound model (with no ocean loading component) of [15] and compared the resultant uplift rates to GNSS data from a small subset of GNSS stations around Greenland [65]. The rates of GIA found for KELY, QAQ1 and SCOR derived from [16] were found to be most compatible with the GIA prediction supplied using the Huy2 model of [38].
The methods used by [16,64,66], contrary to the forward modelling approaches of [24,38], eliminate the uncertainty of Earth model configuration and deglaciation history but are Fig. 4 Comparison of present rates of GIA-related geoid height variation in Greenland calculated using the deglacial chronology and viscoelastic Earth structure of [24] (Huy3), a, and b the ICE-6G_C (VM5a) global deglaciation history and Earth model VM5a [25,70]. White lines on the plots denote the zero contour between areas of geoid uplift (subsurface mass input) and subsidence (subsurface mass output ). Geoid height trends are contoured every 0.1 mm/year dependent upon the accuracy and resolution of geodetic datasets such as ice height measurements gained from altimeters (e.g. ICESat), assumptions regarding densification of the snowpack and the error this propagates within later altimetry measurements. Nevertheless, increasing numbers of studies demonstrate that recovery of GIA trends from remote sensing of the cryosphere is achievable.

Conclusions
The importance of an accurate GIA correction when interpreting geodetic data has been demonstrated in this review paper. However, for the Greenland Ice Sheet, although there is some agreement regarding the magnitude of this correction, and the broad spatial pattern of Greenland-wide GIArelated uplift and subsidence, the contribution of subsurface mass movement as a result of solid Earth adjustment to present-day mass change trends measured by GRACE is becoming increasingly insignificant as the Greenland Ice Sheet enters a phase of increasing global temperatures. In this respect, the priority to further constrain Greenland's GIA contribution to present-day geoid variations is not high, but differences in GIA-derived uplift rates between models do exist at the regional scale as demonstrated by Figs. 1, 2 and 3, which pose an issue when deriving changes in ice height using remote sensing techniques. Estimations of the GIA component in Greenland (regional and Greenland-wide) would benefit from advancements made in the following areas: 1. In light of increasing computational power, the use of higher-resolution (<5 km) glaciological process models could be used to better capture retreat-readvance cycles such as that in the south-west of Greenland. Currently, areal extent of late Holocene retreat is estimated by deglacial models with a spatial resolution that is one fourth of the modelled magnitude of margin readvance. 2. Application of a 3D Earth structure to calculate GIA to investigate the degree to which a regionally varying Earth structure affects the present-day uplift rates. The approach of using a heterogeneous rheology to model GIA has recently been applied to the Antarctic Ice Sheet [63,67]. 3. In northern Greenland, observed uplift rates are lower at present compared to other regions in Greenland. Therefore, it is more important to produce accurate GIA corrections for this region. The treatment of the Innuitian Ice Sheet in GIA models, especially those that 'bolt' together a global ice distribution with a thermomechanical ice sheet constrained regional deglacial history (e.g. Huy3), is often left with issues with respect to how and if the Innuitian Ice Sheet coalesced with the Greenland Ice Sheet in the past. This has consequences for predicted uplift rates in the far north of Greenland. For example, in Huy3, the Innuitian Ice Sheet component is sourced from the North American ice complex in [40], with no explicit glaciological interaction with Huy3 considered. In [68], a deglaciation history for the Ellesmere/Queen Elizabeth Island region (named SJD15) produced an improved fit to a suite of RSL and GNSS data compared to the ICE-5G deglaciation history in Greenland and therefore should be implanted in regional GIA models in the future.