MASS BALANCE CHANGES AND ICE DYNAMICS OF GREENLAND AND ANTARCTIC ICE SHEETS FROM LASER ALTIMETRY

During the past few decades the Greenland and Antarctic ice sheets have lost ice at accelerating rates, caused by increasing surface temperature. The melting of the two big ice sheets has a big impact on global sea level rise. If the ice sheets would melt down entirely, the sea level would rise more than 60 m. Even a much smaller rise would cause dramatic damage along coastal regions. In this paper we report about a major upgrade of surface elevation changes derived from laser altimetry data, acquired by NASA’s Ice, Cloud and land Elevation Satellite mission (ICESat) and airborne laser campaigns, such as Airborne Topographic Mapper (ATM) and Land, Vegetation and Ice Sensor (LVIS). For detecting changes in ice sheet elevations we have developed the Surface Elevation Reconstruction And Change detection (SERAC) method. It computes elevation changes of small surface patches by keeping the surface shape constant and considering the absolute values as surface elevations. We report about important upgrades of earlier results, for example the inclusion of local ice caps and the temporal extension from 1993 to 2014 for the Greenland Ice Sheet and for a comprehensive reconstruction of ice thickness and mass changes for the Antarctic Ice Sheets.


INTRODUCTION
The Greenland and Antarctic ice sheets contain more than 99 % of the freshwater ice on Earth, and thus, they are important contributors to global sea level rise.If they were melted completely, sea level would rise by 66 meters.Rapid warming has caused a dramatic decrease of ice in the Artic region during last few decades.Sea ice extent rapidly declined; glaciers sped up, thinned and retreated, and permafrost warmed up and thinned.Arctic warming is further amplified as surface temperature increases due to the melt of the protective snow and ice cover.The Greenland Ice Sheet (GrIS) has lost an average 250 Gt/yr ice annually, equivalent to 0.7 mm/yr sea-level rise, since 2003.There remains much uncertainty in estimating the current mass loss of the Antarctic Ice Sheet (AIS).While vast regions of the West Antarctic Ice Sheet exhibit increasing ice loss in its marine-based regions, increasing snowfall or long-term dynamic thickening of the East Antarctic Ice Sheet might balance this ice loss.
Here we present a comprehensive update of Greenland and Antarctic Ice Sheet surface elevation changes obtained from the laser altimetry record using our Surface Elevation And Change detection (SERAC) software suite (Schenk and Csatho, 2012).We have developed SERAC to derive information from laser altimetry data, particularly time series of elevation changes and their partitioning into changes caused by surface processes and ice dynamics.This partition allows direct investigation of ice dynamic processes that is much needed for improving the predictive power of ice sheet models.
This latest reconstruction of GrIS elevation changes consists of 130,000 ice-sheet elevation change time series, spanning 1993-2014, derived from satellite laser altimetry data acquired by NASAs Ice, Cloud and land Elevation Satellite mission, airborne laser observations obtained by Airborne Topographic Mapper (ATM) and the Land, Vegetation and Ice Sensor (LVIS).This work provides a major upgrade of the results presented in our latest study * Corresponding author (Csatho et al., 2014), by extending the temporal coverage and by including local ice caps and glaciers.We also present the first synoptic reconstruction of ice sheet elevation changes of the AIS, derived from ICESat observations spanning the period of 2003-2009.For both ice sheets, elevation changes are partitioned into components due to vertical crustal deformations, firn-compaction and ice dynamics, and are converted into mass changes.

METHODOLOGY
We have developed the comprehensive Surface Elevation Reconstruction And Change detection (SERAC) method to determine surface changes by a simultaneous reconstruction of the surface topography (Schenk and Csatho, 2012).Although general in the design, SERAC was specifically developed for detecting changes in ice sheet elevation from ICESat crossover areas.It addresses the problem of computing a surface and surface elevation changes from discrete, irregularly distributed samples of the changing surface.In every sampling period, the distribution and density of the acquired laser points is different.The method is based on fitting an analytical function to the laser points of a surface patch, for example a crossover area, size 1 km × 1 km, or within repeat tracks, for estimating the ice sheet surface topography.The assumption that a surface patch can be well approximated by analytical functions, e.g.low order polynomials, is supported by various roughness studies of ice sheets (Van Der Veen et al., 1998, Van Der Veen et al., 2009).Considering physical properties of solid surfaces and the rather small size of surface patches suggest that the surface is only subject to elevation changes but no significant deformations.That is, the shape of surface patches remains constant and only the vertical position changes.It is important to realize that SERAC determines elevation changes as the difference between surfaces, unlike other methods that take the difference between identical points of two surfaces.
Laser points of all time epochs of a surface patch contribute to the shape parameters while the laser points of each time period determine the absolute elevation of the surface patch of that period.Since there are many more laser points than unknowns, surface elevation and shape parameters are recovered by a least squares adjustment whose target function minimizes the square sum of residuals between the fitted surface and the data points.The large redundancy makes the surface recovery and elevation change detection very accurate and robust.Moreover, the confidence of the results is quantified by rigorous error propagation.Moreover, while most other studies do not compute more than an average elevation change rate, SERAC reconstructions incorporate elevation observations from all different time periods covering a crossover area to generate a time series of elevation changes that is referred to the centroid of the crossover region.SERAC also provides error estimates for each reconstructed surface elevation within the time series.Thus, the SERAC method is able to capture important interannual variations that are not recognized by methods that assume a constant elevation change rate.

Antarctica
We reconstructed surface shape and elevation changes from the complete ICESat observational record at all ICESat crossover regions to generate approximately 110,000 elevation change histories for the time period of 2003-2009 (Fig. 1).The crossovers are densest around 86 o S, where ground track spacing approaches 0 km in separation, and less dense in coastal regions where ground track spacing can be as wide as 30 km.While crossovers provide a less dense coverage of the ice sheet than repeat-track measurements, they allow for a robust reconstruction of surface elevation changes, whereas the correlation between cross-track slope and elevation changes between repeat observations limits the accuracy of the along-track solutions.
The resulting elevation change histories are filtered in a multistep process that identifies and removes erroneous change histories based upon spatial boundaries as well as the elevation change data contained in each history.First, all crossover locations outside of grounded ice, with boundaries defined by the MODIS Mosaic of Antarctica (e.g.ice-shelves, and ice-sheets boundaries and grounding lines; freely distributed through the National Snow and Ice Data Center ( (Haran, 2005), https://nsidc.org/data/moa/)), are removed from further processing.To account for observations within the grounding line, but over exposed bedrock, a landcover mask, also derived from the MODIS Mosaic of Antarctica, was used to remove all elevation change histories within 500 m of exposed rock outcrops.
Remaining SERAC time series are examined for goodness of fit between the original laser points and the third order polynomial surface fitted by SERAC as well as for the distribution of elevation observations contributing to the reconstruction, as measured by the spatial fitting error and the condition number, respectively (Schenk and Csatho, 2012).Those time series exhibiting large elevation uncertainties are removed from further processing.In order to assert long enough temporal coverage to adequately examine inter-annual changes over Antarctica, any of the remaining elevation change histories with less than three years of total coverage between the first and last observation were removed from further processing.Additionally, of those remaining histories, any that did not have at least five unique time epochs worth of coverage were also removed from further processing.To account for variations in elevation changes at the calving fronts of outlet glaciers due to ice retreat during the observational time period, any crossover regions whose centroid was within 10 m from the EGM 2008 geoid surface was excluded from further processing.
During the final filtering step, elevation change histories exhibiting localized non-linear dynamic thickness change trends are identified, and removed from further processing.This was done for two reasons: 1.These events reflect very localized changes in ice sheet thickness, in contrast to regional surficial processes or dynamic changes which are more important factors in determining regional ice behavior.2. Short term, large amplitude changes are likely caused by hydrologic events such as the direct redistribution of subglacial water, through infilling and drainage, or by localized changes in subglacial basal shear stress, as in the case of sticky spots (Sergienko and Hulbe, 2011).
After performing the filtering steps described above and removing the curves indicating short-term rapid changes, the remaining 92,000 elevation change histories were corrected for vertical crustal motion using postglacial rebound rates from the w12a Antarctica GIA model (Whitehouse et al., 2012) and partitioned into thickness changes related to surface processes and ice dynamics (Figure 2).Elevation changes due surface processes were taken from a semi-empirical firn densification model (FDM) of Antarctica ( (Ligtenberg et al., 2011)), driven by the Regional Atmospheric Climate Model 2 (RACMO2, (van Meijgaard et al., 2008)).This model accounts for elevation changes due to surface mass balance anomalies and the processes of firn densification into ice over time.Dynamic ice thickness changes are computed as the difference between the total ice thickness change and the ice thickness change caused by surficial processes.It should be noted that other studies have reported the presence of additional instrument-related elevation bias between ICESat campaigns (Hofton et al., 2013, Siegfried et al., 2011, Zwally et al., 2015).However, as no current agreement exists regarding an inter-campaign bias correction (IBC) the data in this study have not be adjusted for it during processing.
The total and dynamic thickness changes were fit with a third order polynomial approximation in the time domain and thickness changes were sampled from these polynomial approximations at the beginning and the end of each Antarctic balance years (March 1 and February 28/29), during the ICESat operational period (March 1, 2003-February 28, 2009) to calculate annual thickness change rates at each crossover location.In order to identify potential outliers from the polynomial approximation process, mean annual thickness change rates were calculated for each year and those crossover with at least one annual value outside three standard deviations from the annual mean was flagged and manually inspected.
It is important to note that the annual changes from 2003-2004 are more limited in spatial resolution than other balance years.The annual change rates were gridded at 25 km postings for each balance year from 2004-2009 (Fig. 3), using Ordinary Kriging in Geosoft Oasis Montaj 8.2.The annual change rate grids were then resampled to 5 km postings, and filtered using a land-cover classification to remove regions of exposed bedrock and regions outside of the grounding line.

Greenland
We used our SERAC system to determine ice sheet elevation changes with up to four repeat ICESat missions per year and with additional ATM and LVIS flights to obtain a comprehensive record for the ICESat period 2003-2009.SERAC is able to handle simultaneously different laser altimetry sensors, such as ICESat, and NASA's airborne sensors ATM and LVIS.This feature enables the increase of the temporal domain of time series and at the same time densifies their spatial positions.The spatial densification is particularly useful in coastal areas where the ICESat tracks are widely separated.Moreover, he magnitude of thickening and thinning patterns of the outlet glaciers require denser sampling than ICESat can provide.Equally important is the temporal extension to 1993-2014, going back in time to the first ATM observations (1993) and forward to current times by including the most recent laser altimetry observations provided by NASA's Operation IceBridge mission.
The final reconstruction consists of 130,000 ice-sheet elevation change time series (Fig. 4), spanning 1993-2014.We followed the methodology described in (Csatho et al., 2014) to compute elevation and mass change grids partitioned into components related to surficial and dynamic changes by using the firn densification model from (Kuipers Munneke et al., 2015).We also examined the recent total and dynamic thickness changes of a the 130 outlet glaciers that we categorized according to their dynamic thickness change patterns in (Csatho et al., 2014).

Antarctica
Fig. 3 shows the annual total, surface process related (FDM) and dynamic thickness change rates between 2004 March 1 and 2009 February 28 and the mean annual dynamic thinning rate is presented in Fig. 6.
These results illustrate the complex nature of ice behavior across Antarctica.In the Pine Island Glacier region (Fig. 6: 1), where warming ocean temperatures quickly erode the ice lying below sea level, surface elevation exhibit widespread regional thinning of over 2 meters per year during the full time period of the study.Further along the coast from Pine Island Glacier, the Sulzberger Ice Shelf demonstrates steady dynamic thinning, in contrast to More interesting is the widespread number of observations of large ice sheet changes extending into the interior of the West Antarctic Ice Sheet where previous studies found no such changes.Additionally, examination of surface elevation change histories identifies regions where large surface changes are likely happening in response to changes in subglacial conditions, including changes in subglacial hydrology due to the episodic infilling and drainage of subglacial lakes (Fig 6: 5), transport of water through subglacial pathways, and locations where the orientation of thickness change anomalies along geologic boundaries (e.g.faults) suggests regions where subglacial geology has a profound influence on ice sheet behavior.

Greenland
The new reconstruction improves the spatial resolution of previous annual thickness change rate grids (Fig. 7).It is especially obvious along the SE coast of Greenland, where the local thinning and thickening patters are confined to individual drainage basins, rather than occur over larger regions as shown in earlier work (Csatho et al., 2014).The annual ice sheet volume change rates determined from the new reconstruction are very similar to our previously published values (Csatho et al., 2014) for the period of the ICESat mission (Table 1).The fact that the increased spatial density of time series does not change the overall volume change rate suggests that the sampling provided by the ICESat ground tracks was adequate for ice sheet mass balance monitoring at the scale of the entire GrIS.Volume change (km 3 /yr balance year Csatho et al., 2014  The extended time series of Greenland outlet glaciers indicates that outlet glacier behavior remained relatively unchanged during the 2012-2014 time period.Thus, glaciers that belong to the category according to the their ice thickness change patterns in 1993-2012 (Fig. 8), still exhibit similar behaviors (Fig. 9).For example Zacharyae Isstrøm in NE Greenland has continued to accelerate in NE Greenland, while Kjer Glacier in W Greenland continues to thin with a uniform rate.However, after a short period of decreasing thinning, Jakobshavn Isbrae has resumed rapid thinning after 2012.It is worth to note that the re-stabilization of Helheim, Køge Bugt and A.P. Barnstorm glaciers in SE Greenland appears to be complete as these glaciers show very little elevation change in the 2012-2014 period.

CONCLUSION
We describe a new methodology which provides improved capabilities of detecting and monitoring surface elevation changes over ice sheets and glaciers where elevation data are collected.We demonstrate the modular flexibility of the SERAC detection framework by applying different firn-densification and glacial isostatic adjustment models between the Greenland and Antarctic Ice Sheets to investigate the uncertainty of ice-sheet mass loss estimates derived from altimetry.
We show the first detailed reconstruction of annual dynamic ice thickness changes in Antarctica, generated on an observation-byobservation basis with quantified error estimates.We highlight regional complexity of dynamic ice thickness behavior in Antarctica for the Pine Island region as well as the Siple Coast, related to the stagnation of the Kamb Ice Stream and continued thinning of Mercer Ice Stream likely caused by long-term dynamic imbalance.We also illustrate the utility of SERAC to describe very localized changes in ice dynamics, particularly those likely related to changes in subglacial hydrology or subglacial geology.
In Greenland, the large spatial and temporal variations of dynamic mass loss and widespread intermittent thinning depicted by our annual elevation change rate maps indicate the complexity of ice sheet response to climate forcing, strongly enforcing the need for continued monitoring at high spatial resolution and for im-proving numerical ice sheet models.Moreover, the sensitivity of the altimetry-based mass balance solutions to the applied climate models and corrections highlights the need for a synergistic application of different methods (gravimetry, GPS, altimetry, inputoutput method).While all of the results described in this work involved crossover-type solutions, efforts are currently underway to densify observations, in both Greenland and Antarctica, using along-track solutions and by fusing different sensor data with laser altimetry.Additionally, the modular SERAC framework should prove valuable in future ensemble-of-models approaches to examine the sensitivity of surface elevation changes to different models of glacial isostatic adjustment, firn-densification, and sensor-based biases.

Figure 1 .
Figure 1.Location of times series locations in Antarctica generated by SERAC at ICESat ground track crossovers.

Figure 3 .
Figure 3. Annual total, surficial process (SMB and firncompaction) and ice dynamic related thickness change rates, Antarctica.

Figure 4 .
Figure 4. Left panel: ATM (2013, red lines; 2014, blue lines) and LVIS (2013, green lines) trajectories over the Greenland Ice Sheet (yellow) and local ice glaciers and peripheral ice caps (blue).Ice sheet, ice caps and glacier boundaries are from Rastner et al., 2012.Ice sheet is shown in yellow and local glaciers and ice caps are blue.Right panel: location of SERAC times series (black dots) generated from the combined 1993-2014 ICESat, ATM and LVIS data sets.

Figure 5 .
Figure 5. SERAC elevation time series illustrating the fusion of ICESat (blue dots) and ATM (red dots) altimetry data from the Jakobshavn drainage basin, West Greenland.The good agreement between the elevations derived from ICESat and ATM observations confirms the excellent accuracy of the laser systems.

Figure 6 .
Figure 6.Upper panel: mean annual dynamic thickness change rate (March 1 2004 -February 28 2009).Numbers mark regions with distinct thickness change patterns described in the text.Lower panel: error of mean annual dynamic thickness change rate map (left) and ice velocity from InSAR measurements (right).

Figure 7 .
Figure 7. Annual total, FDM and ice dynamics related ice thickness changes in Greenland.
However, there is a relatively large difference (30 %) between the previously determined and revised 2003-2004 volume change rates.The first complete 30-day data acquisition phase of ICE-Sat started in September 2003.In our previously reconstruction, we have not extended the record back to 2003 September 1 for those time series that started later in 2003 September with an L2a observation.Therefore, the distribution of the annual rates for the 2003-2004 balance years was very sparse.Our new work has remedied this situation by extrapolating all time series with one month before the first and one month after the last observation.

Table 1 .
Volume change rates of the entire GrIS during the ICE-Sat mission period fromCsatho et al., 2014 and this study.