Greenland Ice Sheet Elevation Change: Direct Observation of Process and Attribution at Summit

Greenland Ice Sheet surface elevation is changing as mass loss accelerates. In understanding elevation change, the magnitudes of physical processes involved are important for interpretation of altimetry and assessing changes in these processes. The four key processes are surface mass balance (SMB), firn densification, ice dynamics, and isostatic adjustment. We quantified these processes at Summit, Greenland, where monthly Global Navigation Satellite System (GNSS) snowmobile traverses measured elevation change from 2008 to 2018. We find an elevation increase of 0.019 m a−1. The sum of the effects of the four processes reproduces the measured elevation time series, in linear trend and in intra‐annual variability. The short‐term variability in elevation is primarily explained by the variability in SMB. Since SMB has not changed significantly over the last century, and the other processes change over longer time scales, the elevation change likely has been ongoing for at least the last 100 years.


Introduction
The Greenland Ice Sheet (GIS) is losing mass at an accelerating rate (Shepherd et al., 2012;Smith et al., 2020). To measure large-scale ice sheet mass balance and derive the corresponding impact on sea level, many studies use surface elevation change dh dt to derive ice mass change (e.g., Johannessen et al., 2005;Thomas et al., 2000;Zwally et al., 1989), by determining volume change and inferring, measuring, or modeling the density of that volume change.
Surface elevation change in the dry-snow zone of an ice sheet is due primarily to four processes. The expression dh dt = w acc + w c + w low + w gia (1) relates change in surface elevation h over time t to the vertical velocity w associated with each process, where positive values correspond to elevation increases. w acc represents accumulation of snow through time, which in absence of any other process increases the surface elevation, with a corresponding ice mass increase. Accumulated snow compacts due to gravitational settling and increases its density. Consequently, firn densification (w fc ), in absence of other processes, decreases the surface elevation, without any corresponding decrease in ice mass. Below the firn, ice flows from the interior toward the margins due to the gravitational driving stress. In areas of flow divergence (w flow < 0), the ice surface elevation decreases, with a corresponding decrease in mass locally. In areas of flow convergence (w flow > 0), the elevation increases, with a corresponding increase in mass. Finally, the bedrock on which the ice sheet rests is not stationary and is subject to Glacial Isostatic Adjustment (GIA, represented by w gia ; Khan et al., 2016), which can either increase or decrease surface elevation, without any corresponding change in ice mass. Importantly, here we neglect the influence of melt, which adds significant complexity. As melt area increases in Greenland, this will become an increasingly important addition for future work.
Multiple studies have made direct and indirect measurements of dh dt and have attempted to relate dh dt to physical processes at the ice sheet surface and below. An early GPS survey at the Greenland Icecore Project (GRIP) site (at 74 • 34.5'N, 37 • 38.5'W, close to the highest point on the ice sheet) using a marker fixed at 80-m depth (Keller et al., 1995) estimated dh dt = 0 at the GRIP site. Hamilton and Whillans (2000) made direct measurements of w flow using GPS markers embedded deep in the firn and calculated the 3-year average dh dt = 2.5 cm a −1 at Summit, 30 km west of GRIP. On the scale of the full ice sheet, Thomas et al. (2000) used measured discharge and published accumulation to calculate dh dt = 1.6-2.1 cm a −1 for the Summit area. McConnell et al. (2000) argued that on a decadal scale, dh dt is dominated by accumulation, neglecting the contribution of w flow . Partitioning elevation change by physical process is important, as accumulation-driven elevation changes add thickness at a lower density, and thus a lower mass, than ice-dynamics-driven changes ). Kuipers Munneke et al. (2015 treated ice-dynamical elevation changes as a residual when comparing w fc + w acc to dh dt measured by laser altimetry and attributed the small increase in elevation in the interior to thickening of the firn layer. To date, studies partitioning elevation change among different physical processes (e.g., Kuipers Munneke et al., 2015;McConnell et al., 2000;Zwally et al., 2011) rely on measurements or model outputs of some of the processes, and treat the other(s) as residuals, particularly w flow . Here, we use measurements and model outputs of all four key physical processes to determine the contribution and associated uncertainty of each to elevation change in the central GIS. This is to our knowledge the first study to treat all four of these physical processes independently and in detail.

Methods
Near Summit Station, Greenland (N72.58 • , W38.45 • ), we measured or modeled all of the terms in Equation 1: Elevation change was measured using precision Global Navigation Satellite System (GNSS) transects, accumulation by measuring the burial rate of bamboo stakes along these transects, densification by numerical modeling (with initial and boundary forcing taken from field measurements), flux divergence calculated from measurements of relative and absolute motion of a surface stake network over the course of 4 years, and GIA from a published model (Khan et al., 2016) based on precise GNSS observations around Greenland's coastal bedrock.

Elevation Change
At Summit, science technicians have measured surface elevation and accumulation monthly since 2007 along a ∼11-km transect used for validation of ICESat Track 412 (Siegfried et al., 2011). Surface elevation along the transect is measured with a dual-frequency Trimble R7 receiver and Zephyr antenna attached to the top of a fiberglass sled (full details in Siegfried et al., 2011;Brunt et al., 2017). We include offsets for both the fixed height of the antenna above the bottom of the sled, and the varying depth to which the sled sinks into the snow surface. On each traverse the technicians measure this sink depth manually. We process the raw GNSS data into positions following the methods of Brunt et al. (2017) using the commercial software package Inertial Explorer (Version 8.60). To determine a single elevation for each monthly transect, we take advantage of the fact that along the transect the science technicians stop at each of 121 stakes to measure accumulation (section 2.2). At each of these stops the GNSS antenna is stationary for between 30 s and 10 min. By combining all data within a 25-m radius of stake locations to a single value, we eliminate any sampling or speed bias. The 25-m radius was chosen to accommodate the total migration of the transect due to ice flow over the time series duration in this study (2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018). We calculate an average elevation of the entire transect for any given month from the 121 averaged points. The time series of elevations thus obtained is shown in the top panel of Figure 1 and described further in the Supplemental Information.

Accumulation
At each of the 121 bamboo stakes along the transect, technicians measure the distance from a marker on the stake to the surface. After correcting for the vertical repositioning of the individual stakes, the change in this distance indicates burial by accumulation, in snow-equivalent units. Note that this is purely height change due to snow accumulation, and not a mass balance estimate, in the absence of a density measurement (section 2.3). To determine a time series averaged over the entire transect, we calculate the snow accumulation by fitting a smooth two-dimensional polynomial surface to the raw stake data at each time step. We used an iterative procedure that removed individual stake accumulation estimates at each time step that were found to be 3 standard deviations away from the computed surface. By averaging into a single surface and assuming that the individual stake measurements are uncorrelated, we can estimate the uncertainty in the resulting accumulation time series using the root mean square deviation between the stake measurements and the fitted surface at each time step. The time series of w acc measured in this way is shown in the middle panel of Figure 1. It is important to note that "accumulation" measured this way is actually "surface elevation change due to accumulation and densification of the top meter of snow," as this is a measure only of the height change relative to a stake. To convert to a water-or ice-equivalent accumulation, a snow density is required. Takahashi and Kameda (2007) and Eisen et al. (2008) showed that when measuring accumulation in this way, the appropriate density to use for conversions is the density at the bottom of the stake.

Densification
We simulated the time-dependent firn densification contribution to surface elevation change from July 2006 to May 2017 using the Community Firn Model (CFM Stevens et al., 2020) with the firn densification formulation of Kuipers Munneke et al. (2015). For the accumulation rate boundary condition, we used the records from the snow-stake measurements described in section 2.2. For the surface temperature boundary condition, we used the temperature record from Vandecrux et al. (2018), which is derived from automatic weather station (AWS) data at Summit and RCM outputs. We initialized the model using depth-density profiles from Lomonaco et al. (2011) (surface to 110-m depth) and Alley et al. (1997) (110 to 330 m). The initial depth-temperature profile was predicted by a separate, longer CFM run for Summit forced using outputs from the MAR3.9 (Fettweis et al., 2017) RCM.
Firn densification models depend critically on the prescribed surface snow density surf , for which there is no reliable parameterization . Snow-pit and shallow-core data from Summit (Montgomery et al., 2018) illustrate that a single universal value for surf does not exist and there is significant spatial variability. For our model simulation, we used a mean surface densitȳs ur = 330 kg m −3 , based on surface density data in the SUMup database (Montgomery et al., 2018).
The CFM predicts that the long-term average firn densification rate from 2006 to 2017 is −0.437 m a −1 . The monthly modeled firn densification rate as a function of time is shown in the bottom panel of Figure 1. The cumulative firn densification is shown in Figure 3. To assess the firn-model uncertainty, we repeated the model run described above but instead used constant surface densities of 300, 315, and 345 kg m −3 . We then repeated these four model runs using the firn densification equations from Herron and Langway (1980)

Flux Divergence
To measure flux divergence, we established a network of metal survey stakes in the vicinity of Summit (Figure 2). We visited the stakes annually and measured their positions with double-differenced precision GNSS using a Trimble NetR8 receiver with Zephyr Geodetic antenna over 4 years. Our stake network was designed to span a flowline and to extend from the east-west flow divide, ∼30 km east of Summit Station, to the west of our elevation transect. As shown in Figure 2, the easternmost stake moved slightly to the east while all other stakes moved west, at increasing rates with increasing distance from the divide, as expected. North and south of our centerline the stake positions also diverged slightly.
To determine surface elevation change due to flux divergence, we use conservation of mass: where u, v, and w are the components of velocity in the x, y, and z directions, respectively.
Because at Summit the ice is frozen to the bed (Cuffey et al., 1995), u and v equal zero at the bed, and we scale our surface measurements of du dx and dv d by approximating a velocity profile using Glen's flow law in laminar flow with an exponent of 3 (further details in the supporting information). We then integrate dw dz to determine w flow :

Glacial Isostatic Adjustment
Many models for GIA have been published (Wake et al., 2016). We choose Khan et al. (2016) as it is a recent, Greenland-specific, observation-driven model. Khan et al. (2016) used measured uplift rates from GNET, a network of 58 GNSS stations fixed on bedrock around the coast of Greenland, to improve mantle viscosity and lithospheric thickness estimates for the Green1 GIA model (Fleming & Lambeck, 2004). While the GIA term is small enough to be nearly negligible, we include it here for completeness. Khan et al. (2016) estimates w gia at Summit to be −0.0014 ± 0.0002 m a −1 .

Time Series of Process-Based Elevation Change
Some of our measurements of elevation change processes result in a single decadal-averaged value for vertical velocity due to the process at hand. In order to construct a process-based elevation time series, then, we must convert from vertical velocity to vertical displacement due to each process by integrating each value with respect to time. We can then sum the effects of each process to generate a process-based time series. The resulting time series are plotted in Figure 3.

Individual Process Rates and Uncertainty
The individual time series of elevation change processes are plotted in Figure 3 and summarized in Table 1. Consistent with expectations, accumulation increases elevation, densification, and ice flow decrease elevation, and GIA has an extremely small effect. When we sum the individual processes, they closely track the elevation signal measured by GPS (Figure 4). The two time series are significantly correlated (R 2 = 0.78, p ≪ 0.001). Note that in addition to the decadal elevation trend, intra-annual changes in elevation are closely tracked by the summed physical processes. Each process rate we consider has an associated uncertainty. We establish an independent uncertainty for each rate (Table 1; see supporting information for details). Note that the combined uncertainty value is larger than the summed dh dt processes. When viewing this combined uncertainty, it is important to note that it is uncorrelated in time. This means that while for any given monthly measurement epoch the measurement uncertainty seems large compared to the signal, the length of the time series allows us to be confident that the close correspondence between the two series seen in Figure 4 is meaningful.

Spatial Variability
The results for each point in the elevation and accumulation time series shown in Figure 1 are derived from the ∼11-km transect that measured elevation and accumulation across the ice sheet. We acknowledge that such averaging blinds our analysis to potential small spatial scale differences in accumulation or elevation change. To assess this, we experimented with multiple ways of subsetting and averaging the data spatially. When analyzed as 121 individual time series, the same trend emerges (standard deviation between the 121 series = 0.002 m a −1 ), though high-frequency noise is greater. When averaged in batches (north end of the transect versus south end of the transect, upslope vs. downslope, or grouped by direction of travel), the same trends and intra-annual features appear. Thus, we conclude that our results are representative of the ice sheet surface changes over our transect.
The elevation change rate calculated from the flux divergence measurements shown in Figure 2 is a single value derived from the velocities of individual stakes. Note that the velocities increase from east to west, and span a larger area than the ground traverse. Differentiating from velocity into strain yields a noisy divergence field owing to small errors in the velocity field. To account for spatial variability in the resulting divergence field, we averaged the divergence over the western half of the domain, which bounds the GPS elevation survey (also shown in Figure 2).

Elevation Change Trend
The elevation time series shown in the top panel of Figure 1 shows a significant (p ≪ 0.001) trend of 0.017 m a −1 . As noted above, ice flow carries our transect downslope at ∼1.8 m a −1 . Using the mean velocity and local slope (GIMP DEM; Howat et al., 2014) to correct for this effect yields a true elevation change trend of 0.019 m a −1 , a small but significant change. This elevation change trend indicates imbalance between the terms on the right-hand side of Equation 1. This imbalance could reflect a recent change in terms with short response times (e.g., accumulation increase due to warming atmosphere) or an ongoing change owing to Note. Note that while flux divergence and GIA were both determined as annual averages, accumulation and densification mean rates are based on the averages of the time series shown in the lower two panels of Figure 1, respectively. The "sum" uncertainty is actually the sum of w acc and w fc errors, combined with the w flow and w gia via root-sum-square (see supporting information for details). a Corrected for advection. terms with a long response time (e.g., w flow is related to the response time of an ice sheet to geometrical forcing factors and is thus likely to vary on a timescale of 10 3 years).
Initially, w fc at a given location was thought to be fixed (Sorge's law; Bader, 1954). However, given that the density profile shape is related to the time-dependent accumulation rate and surface temperature, annually averaged w fc could vary with a response time on the order of 10 years (Li & Zwally, 2015). Given the relatively small changes in both temperature and accumulation over the past several hundred years, including the time period of this study (Hanna et al., 2020;Meese et al., 1994;White et al., 1997), w fc is likely invariant on the timescale of 10 2 years. This leaves w acc as the most plausible driver of the observed dh dt . In section 3.4 below we show that short-term fluctuations in surface elevation are attributed directly to short-term changes in accumulation. Here, our focus is on the long-term trend, which we note is captured well by the red curve in Figure 4.
The mean annual accumulation from our stake measurements is identical to that reported (0.24 ± 0.05 m a −1 ice equivalent; note our value in Table 3.1 of w acc is expressed in snow equivalent) as both a Holocene average and most recent 100-year average by Meese et al. (1994). Additionally, there is no statistically significant trend in our accumulation measurements. While we cannot rule out multidecadal variability, our elevation change rate is close to the maximum expected from stochastic variability (0.02 m a −1 ; Cuffey, 2001). Thus, it is likely that the imbalance we observe between these processes (resulting in the decadal thickening we observe) has been in place for 10 2 -10 3 years. Over this time period Summit may have experienced an elevation increase of between 1.9 and 19 m. This is consistent with the findings of Marshall and Cuffey (2000), who modeled divide migration in the Summit region; importantly, this imbalance may predate the anthropogenic warming of the last 100 years.

Attribution of Elevation Change Variability
The two series match one another closely with no scaling factors (i.e., both series in Figure 4 show similar magnitudes in both intra-annual variability and decadal trend). This similarity indicates that the four major processes we have quantified are indeed the chief influences on elevation change. Figure 3 illustrates the magnitude of the effect of each individual process. The effect of GIA is small enough to be largely ignored. The effect of flux divergence, while much larger, is still significantly smaller than the effect of densification which is almost twice as large. The largest effect from a single physical process is accumulation. As accumulation also has the largest intra-annual fluctuations (the densification series in Figure 3 contains variability but it is a minor contributor to the fluctuations in the red curve in Figure 4; details in the supporting information, Figures S1 and S2), this indicates that the observed intra-annual elevation changes are caused by intra-annual accumulation changes. One implication of this is that altimetry measurements from space, if of sufficient vertical resolution, could be used to determine the short-term variability or even magnitude of accumulation in regions where w fc and w flow are known, could be estimated, or are changing slowly with time. Given recent advances in altimetry (Smith et al., 2020), this is an exciting possibility.

Application Over a Broader Area
This study is a detailed investigation of a single point on the ice sheet. It is natural to consider how the results can be applied to a broader region. Given the in situ nature of some of our measurements, direct upscaling to the remainder of the ice sheet is impractical. There remain, however, other means for measuring or modeling the variables in question. One significant contribution from this study is to confirm that we have indeed captured the key processes at Summit, and their relative importance. Put another way, we have in this study measurements or model results for all five terms in Equation 1, but future work over a broader area, given appropriate means, could measure or model four terms to solve for an unknown fifth. For example, in some regions where velocity change is slow (e.g., ice sheet interior), remotely sensed velocities may not produce acceptable results for w flow . In other regions, RCM output may not capture accumulation well enough to generate an acceptable w acc , particularly short-term fluctuations. Regions with significant melt add complexity and cause difficulty with models attempting to determine w fc . Fortunately, these challenges are not always geographically colocated, and recent advances in laser altimetry (e.g., Smith et al., 2020) produce high temporal and spatial resolution results for dh dt . Improvements in ice thickness estimates from more accurate radar and gravimetric measurements and novel mapping techniques (Morlighem et al., 2017) combined with improved surface velocity estimates from visual imagery feature-tracking (Gardner et al., 2018) and interferometric synthetic-aperture radar, particularly with the upcoming NISAR mission, will also help improve estimates of w flow . Thus, future investigations can use this framework to determine and constrain ice sheet processes that currently defy measurement and model characterization.

Conclusions
Elevation change high in the accumulation zone of the GIS is governed by four processes, which can independently have large effects but when summed produce very small changes. We have quantified each of these four processes at Summit, Greenland, and contemporaneously measured a detailed time series of elevation change using GNSS. We find that the sum of the four processes is significantly correlated (R 2 = 0.78, p ≪ 0.001) with observed elevation change. Accumulation is the single largest factor in elevation change, and also accounts for the largest part of the short-term variability in elevation. The surface elevation at Summit has increased at an average of 0.019 m a −1 over the period 2008-2018. Ice core records indicate that accumulation at Summit is unchanged through the Holocene. Coupled with the relatively long response times of firn densification and ice flow, this suggests the trend has likely been active on century to millennial time scales, predating anthropogenic warming. If so, central Greenland has yet to experience the direct impact of a warming climate which is so dramatic at lower elevations. Financial support for this work was provided by the NASA ICESat-2 Project Science Office and by the NSF, under award numbers PLR-1637003, ARC-1042358, and OPP-0352584. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the National Science Foundation. We thank Associate Editor Matthieu Morlighem, reviewer Vincent Verjans, and three anonymous reviewers for their handling of this manuscript. We are grateful to Jack Saba, Jonathan Chipman, and Ananya Mathur for data processing and analysis in support of this work. Most importantly, this project would not have been possible without the work of many Summit Station Science Technicians who collected the in situ GPS and accumulation data.