Vertical bedrock shifts reveal summer water storage in Greenland ice sheet

GNSS data selection and preprocessing

The GNET consists of about 54 sites, most of which are located along the coast. Here we consider a subset comprising 22 stations that are mounted on bedrock and located close to outlet glaciers. GNSS stations on bedrock are sensitive to all of the processes that cause mass changes on/in/under the ice close to the station. A few stations (for example, UPVK, NUUK and so on) are excluded because they suffer from large mass variation signals from the ocean. In particular, this concerns stations in the northwestern and southeastern parts of Greenland.

We use coordinates and uncertainties from the GNSS daily solutions released by Technical University of Denmark, which are computed with the scheme proposed in ref. ⁴⁰ in the International Global Navigation Satellite System Service 2014 frame, after removing tidal deformations related to the solid Earth and oceans. The vertical displacements observed at GNET sites include not only the loading signals caused by the BWS changes but also several nuisance signals. The latter represent the crustal deformation caused by atmospheric pressure loading (ATM), land water storage (LWS) on ice-free land, non-tidal ocean loading (NTOL), SMB-related signals and thermal elastic expansion (TEM). To isolate variations in the BWS loading, we correct the observed vertical displacements for the nuisance signals using background models (Extended Data Table 1). The ATM, LWS and NTOL loading changes were computed at University of Luxembourg with the scheme proposed in ref. ⁴¹. The SMB contribution was computed by the time integration of signals from the RACMO2.3p2 model³⁵ for the entire ice sheet. The set of SMB signals consists of precipitation, runoff, sublimation and snow drift components. To compute the vertical displacements at the GNET sites caused by a given component, we used the Green’s function method^42,43. Bedrock deformations are also caused by thermal expansion owing to temperature variations. In this study, we account for this process by using the output from ref. ⁴⁴, which is computed from the ERA5 model with a finite element method. Subtraction of all these nuisance signals results in a time series of residual vertical displacements, which is the input for the further analysis.

As an example, Extended Data Fig. 3 shows the subtracted signals and residual vertical displacements for station KAGA located near Jakobshavn Isbræ (Greenlandic: Sermeq Kujalleq), one of the largest GrIS outlet glaciers. SMB and ATM signals are the largest among those subtracted⁴³.

Uncertainties of the mean annual cycle of residual vertical displacements

Next, we estimate uncertainties of the mean annual cycle of residual vertical displacements, which is shown in Fig. 2. We identify various nuisance signals and categorize them into those that can be completely neglected (for example, ice discharge, lakes) and those that need to be included in the uncertainty analysis. Note that we consider separately shallow groundwater in tundra areas and deep groundwater beneath the ice sheet, as detailed in the ‘Contribution of groundwater storage in Greenland’ section. To that end, we compare the ‘primary’ time series computed as explained above with a family of alternative ones. To produce the latter time series, we replace either the input GNSS data themselves or one of the background models of nuisance signals with one or two alternative ones. The alternative set of input GNSS data was provided by the Nevada Geodetic Laboratory at the University of Nevada, Reno⁴⁵. The alternative background models^{44,45,46,47,48} are listed in Extended Data Table 1, along with the primary ones. The resulting ATM and NTOL loading signals were downloaded from the International Mass Loading Service (IMLS⁴⁶). Alternative SMB loading signals were computed in house at DTU Space on the basis of the alternative SMB model (MAR) in the same way as in the case of the primary one. The alternative time series of bedrock deformations owing to thermal expansion were also computed in ref. ⁴⁴ from ERA5 but with a harmonic method.

Noise in the input GNSS data and in each of the nuisance signals is quantified by computing the standard deviation between the ‘primary’ time series and the corresponding alternative one^34,49 (see the caption of Extended Data Fig. 4). The total uncertainty of the mean annual cycle of residual vertical displacements is computed as the root sum square of the standard deviations of noise from all possible sources considered in this study, that is, GNSS data and models of ATM, NTOL, LWS, SMB, TEM and GWS signals.

Notably, the GWS signal in tundra areas is not a part of the SMB models (SMB is only defined over glacial ice). Therefore, this signal is not subtracted from GNSS data when residual vertical displacements are computed. To quantify the impact of that signal, we use hydrological models. We consider the difference between mean annual vertical displacements per calendar month with and without groundwater signal subtracted using PCR-GLOBWB and WGHM models. This basically results in the vertical displacements owing to modelled groundwater signal alone. We show that the root mean square signals computed from the PCR-GLOBWB and WGHM models are 0.25 mm and 0.02 mm, respectively, which is small compared with the BWS signal (see Extended Data Fig. 4).

A comparison of the total uncertainty computed for each station with the signal revealed in the mean annual cycle of residual vertical displacements demonstrates that the latter far exceeds the noise level (see Fig. 2) and, therefore, cannot be explained by inaccuracies in the input data or exploited models. Furthermore, we have chosen the KAGA station as the representative one to demonstrate that the use of alternative models or GNSS data leads to only minor changes in the observed annual cycle of residual vertical displacements (Extended Data Fig. 4).

Contribution of ice discharge

To obtain an upper bound of the contribution of seasonal ice discharge variations to the observed vertical displacements, we use Jakobshavn Isbræ because it shows the largest ice flow velocity seasonality among GrIS outlet glaciers⁵⁰. Note that, although the total change of mass owing to discharge is comparable with that resulting from SMB, its seasonal variability is much smaller⁵¹. We used in our analysis a time series of ice-discharge-related time series of mass anomalies obtained in ref. ⁵². That time series was obtained by combining two datasets: (1) monthly ice discharge at the Jakobshavn Isbræ flux gates (Extended Data Fig. 5a) computed from monthly ice velocities⁵³ and ice thickness values⁵⁴; (2) annual rates of mass-variation change from altimetry and SMB-based vertical displacements. Corresponding mass anomalies were obtained by weighting the latter annual patterns with the monthly ice-discharge estimates. Finally, the solid-rock vertical displacements were derived by means of the Green’s function method^42,43, for the location of KAGA station, which is only about 1 km from the 2015 Jakobshavn Isbræ calving front. We see that the mean contribution of Jakobshavn Isbræ ice discharge to vertical displacements does not exceed 1 mm (Extended Data Fig. 5), which is small compared with the magnitude of residual vertical displacements.

Contribution of lakes

There are many lakes in the Arctic region, none of which were considered in the LWS products used in this study. To further quantify the potential water-mass impact of lakes onto GNSS loading signal in Greenland, we consider three types of lake: lakes in the pan-Arctic region (north of 60° N) in general; supraglacial lakes (SGL) on the GrIS; and proglacial lakes in the coastal part of Greenland.

Lakes in pan-Arctic region

Specifically, the monthly water storage data of some large lakes were directly downloaded from a recently published dataset⁵⁵. Also, we used the intra-annual lake-level datasets from three portals: Hydroweb (https://hydroweb.theia-land.fr/), Global Reservoirs and Lakes Monitor (G-REALM; https://ipad.fas.usda.gov/cropexplorer/global_reservoir/) and Database for Hydrological Time Series of Inland Waters (DAHITI; https://dahiti.dgfi.tum.de/en/), respectively^56,57,58. For months without lake area data in the dataset provided in ref. ⁵⁵, we follow the empirical relationship between lake water level and lake surface area from a similar previous study⁵⁹ by interpolating the missing lake-area data with sampler-level data and converting the lake area and level data to water storage changes.

It is worth mentioning that the considered lakes (38 in total, indicated by blue circles in Extended Data Fig. 6a) are typically large, with 33 of them having a surface area greater than 500 km², accounting for 61.6% of the total surface area of Arctic lakes extracted in ref. ⁵⁵. For the remaining smaller inland lakes⁶⁰ (red circles in Extended Data Fig. 6a), the area data are very limited. Therefore, we use a scale-up strategy⁵⁹ to estimate their impact on mass loading.

Taking KAGA station as an example, we find that the vertical displacement induced by the load of Arctic lake water storage during the study period (2009–2015) ranges from −0.10 to 0.04 mm (Extended Data Fig. 7a). This magnitude is sufficiently small to allow for neglecting the potential influence of these lakes.

SGL on the GrIS

We also assessed the vertical displacement caused by water mass loading from SGL on the GrIS. First, we divided the GrIS into 10 × 10-km equal-area cells and identified all the cells with SGL using approximately 300,000 high-spatiotemporal-resolution optical images from Sentinel-2 and Landsat 8/9 over 2017–2022, as shown in Extended Data Fig. 6b. Then, the monthly SGL area changes during the melting season (that is, May–September) were derived. Note that it is problematic to obtain the monthly SGL area change for the entire GrIS before the launch of Sentinel-2B satellites in 2017 owing to a high cloud contamination. We assumed, however, that the SGL area changes over the 2009–2015 interval were similar to those over the 2017–2022 interval. When converting the SGL area changes to mass changes, we assumed that the maximum depth of GrIS supraglacial lakes is around 8.5 m (refs. ^61,62), which allowed us to estimate the upper limit of water mass changes (Extended Data Fig. 7b). The result shows that, even for the maximal depth of 8.5 m, the magnitude of loading signal caused by SGL mass changes is about 0.3 mm. Thus, supraglacial lakes provide only a minor contribution to the total BWS signal.

Proglacial lakes in Greenland

Similar to SGL, we also consider the impact of proglacial lakes in Greenland using the HydroLAKES database (Extended Data Fig. 6c). In total, 2,687 proglacial lakes are taken into account. Most of the proglacial lakes are smaller than 5 km², whereas the largest one could reach roughly 100 km². The loading signal caused by proglacial lakes is small. At the KAGA station, for instance, it is on the order of only 0.02 mm (Extended Data Fig. 7c).

Contribution of groundwater storage in Greenland

In the context of the impact from groundwater storage⁶³ in Greenland on the loading signals, we distinguish two types of groundwater:

1. 1.

   Shallow groundwater in tundra areas, which results from snow melting and rainfalls there. In principle, this is an ordinary component of the terrestrial water storage, which is described by various hydrological models, including those addressed in the manuscript (PCR-GLOBWB and WGHM); see Extended Data Fig. 7d,e. According to M. Bierkens, who is a developer of PCR-GLOBWB, it is very difficult to model the groundwater accurately in this region, but the model outcome is enough for a first-order estimate of its magnitude (personal communication, 2024). In this study, we analyse the uncertainty of groundwater estimates in the ‘Uncertainties of the mean annual cycle of residual vertical displacements’ section.

2. 2.

   Deep groundwater below the ice sheet (its presence was detected by a hydrological well down to the depth of hundreds of metres⁶³). It is a product of ice sheet melting (both at the surface and at the ice sheet base). To the best of our knowledge, little is known about variations in the deep groundwater mass. Here we consider the signal from deep groundwater as a part of the total BWS signal we detect in GNSS data. Unfortunately, elastic loading data do not allow the deep groundwater to be separated from the rest of the BWS.

Validation of the results using GRACE data

To validate the results based on elastic loading data, we compared them with water mass changes extracted from satellite gravimetry data. We used four GRACE-based mascon data products. Three of them are off-the-shelf products that were released by: (1) the Center for Space Research (CSR RL06 v02) of the University of Texas at Austin⁶⁴; (2) the Jet Propulsion Laboratory (JPL RL06 v02)⁶⁵; and (3) the Goddard Space Flight Center (GSFC RL06 v1.0)^66,67,68. The fourth is the mascon product computed in house^34,69,70. We corrected the time series for glacial isostatic adjustment using the model in ref. ⁷¹, subtracted the SMB signal and detrended the results. Both the GRACE-based monthly estimates of BWS and the monthly water-related elastic displacements were averaged over entire Greenland. On this basis, BWS mean seasonal cycles were obtained. The mean of the four GRACE-based mean seasonal cycles, as well as the mean seasonal cycles of elastic displacements, are shown in Fig. 3 for entire Greenland (the GRACE scale is inverted for clarity). The standard deviations δ_GRACE of the GRACE-based estimates were computed as

in which d_i (with i = 1, 2, 3, 4) represents the root mean square difference between the estimates based on the ith variant of the mascon data product and the mean ones. The factor 1/2 is present because of the fact that we address the error in the mean of the four time series, rather than errors in the individual ones.

The seasonal cycles based on GRACE and GNSS data are remarkably similar. In particular, both datasets show a mass increase from May to July/August, with a subsequent mass loss until February the following year. Minor differences between the GRACE-based and GNSS-based results can be explained by random errors and the different spatial resolutions of these two data types. We interpret the revealed similarity as a confirmation that both types of data show the signal of the same origin: an accumulation and release of water within the GrIS.

Analytic model of the BWS signal in GNSS elastic loading data

In this section, we explain the analytic function proposed to describe the BWS signal in the residual vertical displacement data. At most of the GNSS stations, a prominent signal in the residual vertical displacements is an upward trend, which reflects a slow mass loss caused by ice discharge. Superimposed to this slow mass loss, there is a seasonal BWS signal at many GNSS sites, which peaks in the middle of the summer (Fig. 2). We hypothesize that this signal is because of BWS. Our interpretation stems from the fact that buffered water, which is not refrozen in place, is not a part of SMB and, therefore, is not described by SMB models.

The total mass balance for grid cell j of an SMB model can be represented as³⁴:

in which M^(j)(t) is the total mass, D^(j)(t) is ice discharge, B^(j)(t) is SMB and S^(j)(t) is the BWS. In RACMO2.3p2, SMB is computed as a combination of four components²: SMB = P − SU − ER − R, in which P is total precipitative flux (sum of snowfall and rainfall), SU is sublimation, ER is erosion of snow by divergence of the drifting snow transport and R is runoff.

Figure 2 shows that, at many stations, the trend after the end of the melt season does not appear as a continuation of the trend observed before the melt season. This implies that the net change in the ice mass during the melt season is different from what could be expected if ice discharge were the only cause of that mass change (at least, under the assumption that the ice discharge is constant over the considered time period). Because the effect of temporal variations in ice discharge is probably minor (Fig. 3), we believe that there must be another explanation for this mismatch. We hypothesize that it can be explained by a difference between the mass loss owing to actual water runoff R and runoff R₀ from the SMB model. We assume that the true runoff, R, is related to the modelled runoff, R₀, as

in which factor ϵ_k accounts for errors in the modelled runoff and is estimated per year. This factor is assumed to be spatially invariant in the vicinity of a given GNSS station. Then, the true SMB can be represented as SMB = P − SU − ER − R = P − SU − ER − (1 + ϵ_k)R₀ = SMB₀ − ϵ_kR₀, in which SMB₀ = P − SU − ER − R₀ is the modelled SMB. For a cell j, the true SMB contribution B^(j)(t) to the total mass balance can be written as

so that the mass balance equation given by equation (2) can be rewritten as:

The integration of equation (5) over time yields the mass at a given time t relative to the mass at the initial epoch t₀. Let us assume for simplicity that t₀ coincides with the beginning of a calendar year. Then, the result of the integration is mass variation M^(j)(t), which is defined under the assumption that M^(j)(t₀) = 0:

in which K(t) is the number of the year containing the current time t, \({t}_{{k}_{0}}\) is the time at the beginning of the kth year and \({t}_{{k}_{{\rm{e}}}}\) is either the time at the end of the kth year (if k < K) or the current time t (if k = K). Notice that the integration of the term ϵ_kR₀^(j)(t) over time implies that the modelled runoff is integrated separately in each year; the results are scaled with factors ϵ_k and summed over all years (up to the year containing the current time t).

We assume that the contribution of the ice discharge does not change over time: −D^(j)(t) = A^(j). Then, equation (6) can be simplified as:

in which C^(j) = A^(j)t₀ − S^(j)(t₀). The mass variations given by equation (7) result in vertical elastic deformations of the solid Earth. Assuming a linear relationship between mass variations and elastic deformations, we can rewrite equation (7) in terms of vertical elastic deformations at the location of a GNSS station:

in which different signals in the time series of vertical elastic deformations (denoted with lowercase letters) are associated with the corresponding mass signals (denoted with capital letters in equation (7) and before). Technical details of the transformation of surface mass load into elastic vertical deformations can be found in the section below entitled ‘Computation of elastic vertical deformations and spatial sensitivity of GNSS loading data’. The left-hand side of equation (8) contains the residual GNSS measurements before the correction for the SMB signal. On the right-hand side, we see, among others, the elastic deformations associated with the SMB model, \(({\int }_{{t}_{0}}^{t}{b}_{0}(\tau ){\rm{d}}\tau )\), and the elastic deformations associated with the computed runoff \(({\int }_{{t}_{{k}_{0}}}^{{t}_{{k}_{{\rm{e}}}}}{r}_{0}(\tau ){\rm{d}}\tau )\). Both signals are computed on the basis of the RACMO2.3p2 model output. The unknown constant factors a, c and ϵ_k can be estimated using least squares from the observed time series m(t). The only term that requires a further discussion is the signal s(t) associated with the BWS.

Let us consider the total BWS S(t) in the drainage basins located around the current GNSS station (more specifically, in the drainage basins that substantially affect the elastic deformations at the current GNSS station). Temporal variations of that BWS mass are equal to the difference between the total runoff R(t) in those drainage basins (which describes the rate of liquid water production) and the rate Q(t) (which describes discharge of water from the drainage basin into the ocean):

Because transport of water from the location of production to the location of discharge can take weeks or even months, the BWS S(t) can be substantial. Let us assume that the discharge into the ocean is proportional to the BWS (that is, drainable water storage), which is a commonly used assumption in hydrology⁷²:

in which β is a certain constant proportionality coefficient. The substitution of this expression into equation (9) yields:

or

in which T_st = β^−1. Assuming that the runoff R(t) is given, we can readily find the solution of the differential equation given by equation (12) as

From this expression, it follows that the parameter T_st can be interpreted as the characteristic time of BWS. We refer to it as the ‘water storage time’.

The first term in equation (13) represents the impact of the initial BWS, S(t₀). Because t₀ is assumed to coincide with the beginning of a calendar year (here year 2009), we set S(t₀) = 0 (to minimize the impact of this assumption, we ignore year 2009 in the subsequent analysis as an ‘initialization’ year). Then, we obtain:

Equation (14) allows us to introduce an approximate relationship between runoff and BWS in terms of loading signals.

If the spatial pattern of runoff was similar to the spatial pattern of BWS, equation (14) could be rewritten in terms of the loading signal directly (as it was already done when equation (8) was introduced). In practice, of course, this is not the case. Notably, vertical deformation reduces as the distance between the GNSS station and the location of the surface load increases⁷³. Therefore, the buffered water signal measured at a GNSS station can show systematic deviations from the values predicted on the basis of the runoff in line with equation (14). To take this effect into account, we introduce an empirical scaling factor θ per GNSS station, so that an approximate relationship between runoff and BWS in terms of loading signals is:

Taking into account that the true runoff is defined as a scaled variant of the modelled one (see equation (3)), we can rewrite the expression above as:

After the substitution of this equation into equation (8) and the isolation of the known terms on the right-hand side, we finally obtain:

By considering this equation for all times t within the interval under consideration, we can form a system of nonlinear equations containing n + 4 unknown parameters per GNSS station: c, a, θ, T_st and ϵ_k (k = 1,…, n), in which n is the number of years in the considered time interval.

To estimate all the unknown parameters, an iterative least-squares adjustment could be directly applied. In the course of a preliminary study, we realized, however, that there is a trade-off between water storage time T_st and the mean value of the corrections ϵ_k. Each of the two can be used to explain signals in the input data, whereas an attempt to estimate them simultaneously frequently results in unphysical estimates (for example, a nearly zero water storage time). To solve that problem, we have introduced a constraint that forces the mean value of the corrections ϵ_k to be equal to zero. This constraint can be interpreted as an assumption that the runoff estimates provided by the SMB model are correct on average in the study period (even though they still may contain errors in individual years). Furthermore, we watch that ϵ_k ≥ −1. A violation of this inequality implies that the true runoff in a given year is negative. Then, in line with equation (14), the estimated BWS becomes negative as well. Of course, all of that is unphysical. In a few cases when this still happens, we refrain from estimating the true runoff. Instead, we force the corresponding estimates of ϵ_k to be exactly equal to −1. This corresponds to a zero runoff and a zero BWS.

Once all of the unknown parameters are estimated, the vertical displacements caused by variations in BWS can be readily computed with equation (16). For a comparison, we also present the vertical displacements computed under the assumption that the SMB-based runoff estimates are correct, so that ϵ_k = 0. Under this assumption, equation (16) simplifies to:

To distinguish the BWS-related vertical displacements computed with equations (16) and (18), we call them ‘calibrated’ and ‘uncalibrated’, respectively.

Notably, the last term on the left-hand side of the functional model given by equation (17) describes the accumulation and discharge of the BWS, which is a short-term process. This signal declines exponentially after the end of the melt season (that is, after the runoff-related signal r₀(t) turns to zero). This is fully consistent with the behaviour of BWS, which is primarily produced as a result of ice/firn/snow melting and ends up as discharge into the ocean. This term controls, in the first instance, the estimated water storage time T_st and empirical coefficient θ. By contrast, the third term on the left-hand side of equation (17) describes the long-term effect of inaccuracies in the runoff estimated as part of the SMB. The effect of these inaccuracies does not vanish in the course of time. The estimated corrections ϵ_k are mostly controlled by this term, whereas the impact of the fourth term on those estimates is minor. To demonstrate that, we have considered a modified functional model that lacks the BWS-related signal (that is, the aforementioned fourth term):

Of course, such a functional model is not applicable in the course of the melt season and immediately thereafter. Therefore, we have limited the input data time series to either 6 months per year (November–April) or even to 4 months per year (December–March). This allowed us to obtain two alternative estimates of corrections ϵ_k (as well as those based on the original functional model given by equation (17)). The mean of the three estimates, as well as the associated standard deviation, is reported in terms of scaling factor (1 + ϵ_k) per station per year in Extended Data Table 4. We can see that the standard deviation in most cases is less than 0.15. Only one GNET station—GMMA—may show standard deviations larger than 0.3. This means that the obtained estimates of corrections ϵ_k are sufficiently robust; intrinsic uncertainties associated with the spatial distribution of BWS during the melt season have only a minor effect.

Accuracy of the obtained estimates

Uncertainties for all reported estimates have been quantified. This concerns both water storage times (Fig. 6 and Extended Data Table 3) and the scaling factors to be applied to the SMB-model-based runoff estimates (Fig. 5 and Extended Data Fig. 2d). The input for the error propagation procedure was defined as the standard deviation of errors in the residual displacements. To that end, the post-fit residuals obtained after fitting the residual displacements with the analytic function given by equation (17) were considered as realizations of the aforementioned errors, which were assumed to be uncorrelated. In view of a nonlinearity of the inversion procedure and a skewness of the resulting error probability density functions, the error propagation was implemented by means of Monte Carlo simulations⁷⁴, in which 100 realizations of random errors were generated for each station. The resulting uncertainty intervals were quantified with the 1st quartile and the 3rd quartile as the lower and upper bounds, respectively.

Computation of elastic vertical deformations and spatial sensitivity of GNSS loading data

Because the solid earth is an elastic body, it experiences vertical deformations in response to a changing surface mass load. Let that load be defined in terms of equivalent water height (EWH) as h(φ, λ), in which φ and λ is geographical colatitude and longitude, respectively. We compute the resulting vertical deformations⁴² U(φ, λ) as the convolution of the surface mass load sources m(φ″, λ″) with the Green’s function G(ψ):

with

in which ψ is the spherical angular distance between the points (φ, λ) and (φ″, λ″); σ is the integration area; a_e and m_e are the mean radius and mass of the Earth, respectively; \({h}_{n}^{{\prime} }\) is the nth degree load Love number; and P_n(cosψ) are fully normalized Legendre polynomials. As follows from equation (20), vertical deformations are proportional to the magnitude of surface mass load.

To demonstrate a possible spatial variability of elastic vertical deformations, we have considered several surface mass loads, each of which is homogeneously distributed over a disc of a given radius. The thickness of each disc is defined in terms of EWH such that the deformation at its centre reaches 5 mm, which is similar to the magnitude observed in real data (Fig. 2). We can see (Extended Data Fig. 8a) that the spatial variability of the resulting deformations strongly depends on the spatial extent of the surface load. We believe that a disc of 200 km (or larger) radius, with a total mass of (at least) 35 Gt, gives the best approximation of the actual surface load distribution. This is because the total BWS mass per drainage system estimated from GRACE satellite gravimetry data is on the order of 20–40 Gt (ref. ³⁴), whereas the shape of the actual surface load distribution probably resembles half a disc rather than a disc (most of the GNSS stations are located near the coast, whereas the surface load over the ocean is nearly constant). This simple example also shows that the position of any realistic surface load relative to a given observation point must move laterally by at least a few tens of kilometres to change the elastic vertical deformation at the observation point substantially. This means, for instance, that a local redistribution of meltwater within the firn layer cannot affect deformations observed at a given GNSS station.

Also, we performed a sensitivity study to clarify how mass changes over the entire GrIS, including outlet glaciers hundreds of kilometres away, may affect the observed mass loading signal in practice⁷⁵. By taking the detrended SMB in July 2012 (Extended Data Fig. 8b) as the GrIS mass change signal, we analyse the sensitivity of GNSS loading observed at the KAGA station as an example. The processing strategy is to take the KAGA station as the centre, create ring-shaped zones outward with a step width of 50 km, calculate the vertical displacement caused by SMB inside each ring and then normalize with the SMB signal from entire Greenland. The obtained sensitivity curve (Extended Data Fig. 8c) reveals a substantial contribution from the ‘near field’: SMB changes within 200 km from KAGA contribute about 80% to the detrended SMB loading displacements. There is little sensitivity to SMB beyond 500 km from KAGA. This is evidence of a relatively high spatial resolution compared with GRACE data: the spatial resolution of the latter is on the order of 400 km (in terms of wavelengths, when the spherical harmonic expansion to degree 96 is considered) or even worse.