Last modified on Dec 01, 2023 17:38

Contributors: L. Gilbert (University of Leeds), S. B. Simonsen (Technical University of Denmark)

Issued by: University of Leeds / Lin Gilbert

Date: 26/04/2023

Ref: C3S2_312a_Lot4.WP2-FDDP-IS-v1_202212_SEC_ATBD-v4_i1.1

Official reference number service contract: 2021/C3S2_312a_Lot4_EODC/SC1

History of modification

List of datasets covered by this document

Related documents

Acronyms

General definitions

Backscatter: The portion of the outgoing radar signal that the target redirects directly back towards the radar antenna.

Baseline: A combination of processor versions, auxiliary data and other needed enablers that allows the generation of a coherent set of Earth observation products.

Climate Data Record (CDR): A time series of measurements of sufficient length, consistency and continuity to determine climate variability and change.

Retracking: Finding the range from the instrument to the point of closest approach on the ground by examining the shape of the radar echo.

Surface Elevation Change (SEC): The surface elevation of a point on an ice sheet is the height of the ice sheet surface above a reference geoid (a hypothetical solid figure whose surface corresponds to mean sea level and its imagined extension under land areas). Increase in surface elevation over time at a given location indicates a gain of ice or snow at that location, and conversely decrease indicates a loss. The surface elevation change product provides the rate of change given at monthly intervals at each location on a grid covering the ice sheet. The definition of the grid projection includes the geoid used. Given the rates of change, absolute change can be calculated for any time period.

Tracking: Retrieving the radar echo from a given radar pulse. 

Scope of the document

This document is the Algorithm Theoretical Basis Document for version 4 of the Surface Elevation Change (SEC) products made as part of the Copernicus Ice Sheets and Ice Shelves service. The products contain geographically-gridded timeseries of the rate of change of ice sheet and ice shelf surface elevation in Antarctica and Greenland, from 1992 to the present.

The document describes the satellite-mounted instruments, auxiliary datasets, auxiliary models and basic algorithms used to create the data products.

The products are hosted on the Copernicus Climate Data Store at https://cds.climate.copernicus.eu/cdsapp#!/dataset/satellite-ice-sheet-elevation-change?tab=overview.

Executive summary

The service addresses three essential climate variables (ECVs) by providing four separate products.

• Ice velocity is given for Greenland in product WP2-FDDP-IV-CDR 
• Gravimetric mass balance is given for Greenland and Antarctica in product WP2-FDDP-GMB-CDR 
• Surface elevation change is given for: 
  •  Antarctica in product WP2-FDDP-SEC-CDR-AntIS
  • Greenland in product WP2-FDDP-SEC-CDR-GrIS

We document here the inputs, processing stages and outputs of the two Polar region SEC version 4 (v4) products. 

Although the two products take the same input datasets and have the same output format, they necessarily use different map projections and processing methods, and so have been split to avoid possible confusion. Section 1 describes the Antarctic product and section 2 the Greenland product. Both sections follow the same format.

Subsection 1 describes the instruments used. Radar altimeters mounted on six different satellites provide the primary input data, and are discussed here, along with issues relating to the datasets they provide. As their datasets are continually being improved and new versions released, changes newly incorporated into the SEC v4 products are discussed. Links to further information and a tutorial on radar altimetry are given.

Subsection 2 describes the input data in more detail, and also describes the auxiliary data used during SEC processing. Tables of mission and instrument characteristics for each input dataset are given. The uses and sources of auxiliary data are listed. The types of validation data used and their sources are discussed.

Subsection 3 describes the algorithms used in the processing. It gives an overview of the methodology, highlighting changes between this version of the product and previous versions. Processing broadly follows a three-step procedure. The first step is derivation of surface elevation change over time, including compensation for instrument effects. The second is cross-calibration, which merges the results from each mission into one, multi-mission, consistent time-series. The third is derivation of a rate of change from each change timeseries. The algorithms for each step are described, followed by a discussion of the algorithms used to derive uncertainty measures at each step, and how these are combined into a total uncertainty value.

Subsection 4 describes the output of the processing chains, the product itself. The description includes the output file type, the difference between a Climate Data Record (CDR) and an interim CDR (ICDR), the file contents and format, and a table summarizing the variables contained.

1. Surface elevation change, Antarctica

1.1. Instruments

In both hemispheres, the SEC product uses the altimetry record from a series of overlapping European Space Agency (ESA) radar altimetry missions. The record started in 1992 with two European Remote-Sensing Satellites (ERS-1 and ERS-2) and Envisat, whose missions are now finished, and will extend into the future with CryoSat-2 and Sentinel-3 A and B, which are still in progress. Sentinel-3 does not yet apply the most optimal processing scheme for land ice. This results in a loss of tracking by the radar altimeter as the satellite approaches the complex topography of the ice sheets from the oceans. ESA has developed a specialised land-ice processor for the Sentinel-3 mission, with the product scheduled for release at the end of 2022. Unfortunately, that is too late for it to be incorporated in the v4 SEC product, but when it becomes available it is intended to be used.

In the v4 dataset one input data change from v3 has taken place. The CryoSat-2 processing baseline was changed from D to E in August 2021, and no reprocessing of the earlier data was made. Consequently, both baselines are included in the SEC product. These two baselines do not differ in major ways, see Geminale (2021). The most relevant changes are progressive improvements to the land ice retracking in both altimeter operating modes and resolution of an issue with the computation of backscatter power in low resolution mode (LRM) which caused a small drift in expected values. LRM is used only over the undynamic interior of the Antarctic Ice Sheet. A full mission reprocessing to baseline E is expected by the end of 2022, unfortunately, again, that is too late for it to be incorporated in the v4 product.

The Antarctic region surface elevation change dataset is derived from the crossover method (see Section 1.3 below and Zwally et al, 2012b). Only certain mission phases are suitable for use in crossover analysis, due to their orbital configurations. The mission's orbit inclination governs the area of the Antarctic region that can be observed, as it mandates the latitude of the most southerly observation possible.

Table 1 lists the missions/mission phases used as input to both the Antarctic and Greenland SEC datasets.

Table 1: Summary of mission-level parameters relating to datasets used as input to both the Antarctic and Greenland SEC datasets.

For further information on these missions and their instruments see the ESA's Earth Observation portal¹.

Radar altimetry is described in detail in a PDF tutorial document².

1.2. Input and auxiliary data

1.2.1. ESA Radar altimetry level 2 products

The following tables provide information on the primary input data sources for the SEC product, the satellite-mounted radar altimeters. Only operating modes that are used in the Antarctic and/or Greenland SEC products are listed.

Table 2: ERS1 altimeter parameters and data source information

Table 3: ERS2 altimeter parameters and data source information

Table 4: Envisat altimeter parameters and data source information

Table 5. CryoSat-2 altimeter parameters and data source information

Table 6: Sentinel-3A altimeter parameters and data source information

Table 7: Sentinel-3B altimeter parameters and data source information

1.2.2. Auxiliary data

1.2.2.1. Surface type mask

The surface type mask is used to split the region into areas of interest. A four-type surface classification is used, i.e. ice sheet, ice shelf, ice rise or island, and ocean. This classification scheme was prepared by NASA's ICESat mission team, by combining data from NASA's Terra satellite's Moderate Resolution Imaging Spectrometer (MODIS) and ICESat satellite's Geoscience Laser Altimeter System (GLAS). The mask is described by Zwally et al, 2012a. Further information and a freely downloadable copy of the mask are available from:

• https://earth.gsfc.nasa.gov/cryo/data/polar-altimetry/antarctic-and-greenland-drainage-systems ³

1.2.2.2. Slope model

Radar altimetry depends on the echo of a radar pulse from the ground surface below. The echo return begins from the highest point encountered within the pulse-limited footprint (see Table 2 to Table 7), which is not necessarily the point immediately below the satellite. In low slope regions, corrections are applied to locate the echo return point. In high slope regions the topography makes these corrections uncertain. The Slater et al, 2018 slope model is used to both filter out input data from areas of very high slope and to flag areas of high surface slope. It is freely available via the ESA CryoSat Operational Portal⁴.

1.2.2.3. Drainage basin mask

Drainage basins denote the ice sheet regions feeding each area of ice shelf or ocean outlet. The international Ice sheet Mass Balance Inter-comparison Exercise (IMBIE) 2 drainage basin outlines are used to split the Antarctic land mass into 27 basins. For details on IMBIE see Shepherd et al, 2012. The drainage basins are described in Zwally et al, 2012a, and freely available from:

• https://earth.gsfc.nasa.gov/cryo/data/polar-altimetry/antarctic-and-greenland-drainage-systems ⁵

1.2.2.4. Glacial isostatic adjustment

The Antarctic land mass is slowly uplifting as it reacts to the removal of some of its ice cover after the last glacial maximum approximately 10000 years ago, which is known as glacial isostatic adjustment (GIA). The Antarctic Ice Sheet (AIS) SEC product includes a correction term to account for this uplift, taken from the IJ05 R2 model, described in Ivins et al, 2013. The Antarctic solution is planned to be freely available as an ancillary dataset from the IMBIE project⁶.

1.2.2.5. Tide model

The service provides information from the Antarctic ice shelves, which are affected by tidal motion. Although the radar altimeter products used as input include corrections to the elevation measurement due to tides, the corrections are applied to locations based on land masks that are too low resolution to accurately outline the ice shelves. As a consequence, rectangular artefacts due to the land mask are visible in the input data.

To compensate for this, tidal corrections in the input data are removed, and a new set calculated using the Circum-Antarctic Tidal Simulation (CATS) 2008a polar tide model. The model is an update to that described in Padman et al, 2002. The CATS 2008a polar tide model is freely available from:

• https://www.esr.org/research/polar-tide-models/list-of-polar-tide-models/cats2008/ ⁷

1.2.2.6. Antarctic Ice Sheet ice velocity map

During validation a bias value is calculated to compensate for the tendency of the Operation IceBridge Airborne Topographic Mapper (ATM), which provides the validation data, to preferentially sample fast-thinning ice. This bias is calculated from the 1 km by 1 km Berkeley – Ice Sheet Initiative for Climate Extremes (BISICLES) ice velocity map⁸, see Rignot et al, 2011.

1.2.3. Validation data

Independent estimates of the rate of surface elevation change at discrete locations and over specific time periods are provided by the ATM, a scanning laser altimeter flown on board aircraft by Operation IceBridge (Studinger 2014).

The validating dataset used for both hemispheres SEC products is a level 4 product, IceBridge ATM L4 Surface Elevation Rate of Change V001. This can be obtained free of charge on registration, fromhttps://icebridge.gsfc.nasa.gov/

1.3. Algorithms

1.3.1. Introduction

The product requirements are a stack of gridded surface elevation change rates from the Antarctic region, at 25 km resolution, from the start of the ERS1 mission to the present. The grids should be given at monthly intervals, and flag grids for steep terrain and terrain type should also be provided. The change rate units should be m/year. The accuracy target is 0.1 m/year and the stability target is also 0.1 m/year.

To make the SEC product, elevation data from the radar altimetry satellite missions listed in Section 1.2.1, over a long period of repeated observations, is necessary. The Antarctic region, including ice sheets, ice shelves, ice rises and islands, is divided into regular grid cells on the standard European Petroleum Survey Group (EPSG) polar stereographic projection for Antarctica, EPSG:3031⁹, (see Section 1.3.1 of the related document the Product User Guide [D2] for full details), and a timeseries of surface elevation change is derived in each cell for each mission. The mission timeseries are cross-calibrated in each grid cell to produce a single long-period timeseries, and this is used to derive the rate of surface elevation change in that grid cell for the product. The processing chain is described in the related System Quality Assurance document [D3].

A review of the scientific background of surface elevation measurement is given in the AIS Climate Change Initiative (CCI) ATBD, Section 2.1, see Nagler et al. 2018.

There are several methods for recovering surface elevation change from a timeseries of surface elevations, including crossover and along-track analysis. In this project, the crossover method has been selected as its results are invariant with respect to the period it uses. Along-track analysis fits a model to data over a long period as a whole. When new data is added, the model changes, and this affects the results for the whole period. Crossover timeseries analysis only compares data in pairs from short, specific periods, and thus earlier product values will not be changed by updates to the data timeseries.

In the v4 algorithms only one processing improvement will be included, and that is a review of the filtering limits on elevation change values (described in 1.3.2). This is due to accelerating ice loss around the coastline of Thwaites and Getz glaciers leading into their ice shelves (see Figure 1.1), which was found during the 2022 IMBIE exercise to exceed limits set in previous exercises.

Figure 1.1: Map of Antarctic ice shelves. Credit:  By Paleo nim - Own work, CC BY-SA 4.0, https://commons.wikimedia.org/w/index.php?curid=76810164

1.3.2. Crossover processing

Missions run in cycles, where the same orbital tracks are repeatedly followed. For ERS1, ERS2 and Envisat 35-day long cycles are used in the AIS SEC product. CryoSat-2 has a more complex orbit that 'drifts' with time and repeats only every 369 days but can be split into nearly-repeating pseudo-cycles that are 30 days long. Sentinel-3A/B has 27 day cycles. A timeseries of elevation changes for each mission is built up by comparing each cycle to a reference cycle from the same mission, which is selected for its extensive geographical coverage and high data content. Each mission is treated individually. Figure 1.2 shows the physical layout of a single crossover.

Figure 1.2. Example crossover situation (taken from AIS CCI ATBD, see Nagler et al, 2018)

The radar altimeter measures the surface elevation at closely spaced discrete points along the satellite’s ground track. In a small region around where an ascending and a descending track cross, an average height can be derived for each track, and this height will only be applicable to the very short overflight period. The difference between these heights is the change in surface elevation plus a random measurement error, which includes errors in the altimeter’s height measurements and its orbital position.

where H1 and H2 are the heights at the crossover point and E is the measurement error.

There are several crossover methods available, but we use the dual crossover method (see Wingham et al, 1998) where one cycle of repeating tracks is used as reference, and height changes are derived at all crossovers between that cycle and each of the other cycles. This produces a timeseries of height differences relative to the reference cycle height at each crossover.

where h is height, t is the set of all cycles, t_ref the reference cycle, A_tref and A_t the ascending reference and comparison tracks respectively, and D_tref and D_t the descending equivalents at the same crossover, x.

To calculate the heights, fully corrected elevations from the input datasets are used. Fully corrected elevations have slope correction, which compensates for the source of a radar echo not being directly below the satellite, and a set of geophysical corrections, which are supplied by the input products (see the source dataset documentation listed in Table 2 to Table 7).  The geophysical corrections compensate for tidal effects on water and land at the Earth’s surface, and atmospheric effects on the propagation of the radar pulse¹⁰. The ocean tide corrections are then replaced with a consistent set calculated with the CATS2008a model, as described in section 1.2.2.5.

There may be several crossover locations in each regional grid cell. To regularise the timeseries in each cell, the average dh per cell per cycle-period is calculated. Filtering restricts the crossover locations used to those with a minimum number of nearby measurements (2 per pass). Crossovers are only used where the numbers of measurements from ascending and descending passes are broadly similar (the number in one pass no more than twice the number in the other), to avoid one pass dominating the results.

The standard deviation of the inputs to the averaged dh is noted – this is the uncertainty from the input measurements, which will be used as a component of the total uncertainty of the product.

Two main error sources in the crossover method are the distance the radar signal penetrates the surface, and its volume scattering, which can vary with time. To mitigate this the mission datasets have applied retracking to the received signal, for details see the product handbooks for REAPER, Envisat, CryoSat-2 and Sentinel-3 A and B in the references below. To mitigate further, in this processing each grid cell’s dh timeseries is individually corrected for backscatter power fluctuations, which are influenced by surface height (see Khvorostovsky, 2018). Similarly to the dh timeseries, a timeseries of power differences relative to the reference cycle power at each crossover (dp) can be made, using a similar equation.

where p is backscatter power, A_tref and A_t are the ascending reference and comparison tracks respectively, and D_tref and D_t are the descending equivalents at the same crossover, x.

A fixed five-year period is chosen during each mission, and the correlation between dh and dp is modelled with a linear fit. Five years has been chosen to allow sufficient seasonal cycles to mitigate anomalous periods. Once the period is chosen, and the correlation calculated, it will be applied to the whole mission, so new input data will not affect the backscatter power correction process. In the case of the more-recently established Sentinel-3 data streams, Sentinel-3A is able to use a five year period for the first time in the v4 product, while Sentinel-3B has not been operating long enough for a five year period, and will use three years instead. Neither of the selected ERS1 operational phases lasted five years, so their full extents are used.

To calculate the correlation, in each grid cell we make a linear fit to the model

where dh is height difference to the reference cycle, dp is backscatter power difference to the reference cycle and a₀ and a₁ are constants. We check the Pearson’s R correlation coefficient between dh and dp, and if it is less than 0.5 (indicating that at least a quarter of the variation in dh is accounted for by changes in dp, as used in the AIS CCI project) we assume there is no correlation, and so do not perform any correction.

If correlation exists for a grid cell, we perform the correction on each dh in its timeseries using this equation

where dh_corr is the corrected dh and dh/dp is the gradient of the linear model.

More filtering is used to remove elevation change measurements that are excessively far from their main distribution. Any timeseries from a particular grid cell is removed if it contains less than 10 datapoints. Further filtering is based on a modelled fit to a seasonally cycling signal on an overall linear trend represented by

where dh is change in elevation relative to the elevation of the reference cycle, t is difference in time relative to the reference time and a0 to a4 are constants.

If a timeseries has a datapoint more than 3 standard deviations from the model fit, then that datapoint is rejected. Also, following practise from an internal study of crossover timeseries in the AIS CCI project, an upper limit on the model parameters was set, where the linear trend is 12m/yr and the seasonally cycling signal has a maximum of 1.2m. Datapoints are rejected if their dh represents a change from the reference cycle that exceeds these upper limits.

Finally, timeseries in any cell with a surface slope of more than 5° are removed. In v3 absolute dh values of more than 10 m were filtered out, as this was deemed to be physically unlikely, but experience from IMBIE has shown that these very large elevation changes are possible at the extreme, glacierised edges of the ice sheet. Therefore, for v4 a different filtering step is used, based on the multi-mission timeseries, as described in the next section.

1.3.3. Multi-mission cross-calibration

The v4 cross-calibration method is the same as that for Climate Data Records (CDRs) v2 and v3. It uses a multiple regression algorithm. For general details see Tabachnick and Fidell, 2019.

Before crossover processing begins each mission's dh data is arranged in a stack of grids with a time array corresponding to the stack order. Cross-calibration values are calculated for each grid cell separately. Basin timeseries, which are not part of the Antarctic SEC product but are used to track the basin-level accuracy target, are calculated for each mission using ice-velocity guided interpolation on each grid, and then cross-calibrated using exactly the same algorithm as for the single cells.

The algorithm uses multiple linear regression to fit a model to a multi-mission timeseries. It is assumed that the timeseries follows a cubic polynomial form over time. The independent variables are:

• time
• time squared
• time cubed
• a flag array for each mission except the first.

The regression produces a coefficient for each independent variable. The coefficients for each mission are the cross-calibration bias values. The first mission is unbiased and the rest biased with respect to it. When these biases are applied to the data, it clusters around the cubic polynomial model, as seen in Figure 1.3 below. This shows six mission timeseries in order, from an example grid cell inland of the Thwaites Glacier. On the left they are not cross-calibrated. On the right cross-calibration has been applied so that the earliest mission's position is unchanged. The model polynomial is also shown, as a solid line.

Figure 1.3: Multi-mission timeseries before and after cross-calibration

The regression algorithm in use is the REGRESS function in the IDL v8 software package. This function returns the standard deviation of each term, given input uncertainty estimates. In this case, the input uncertainties are the standard deviations of each datapoint in the timeseries. Thus, estimates of cross-calibration uncertainty can be obtained for each mission, assuming no error on the first.

The model is not used to add extra datapoints into any data gaps. The biases calculated as necessary for cross-calibration are added to pre-existing datapoints only.

For the v4, a final filtering step, calculated in each cell from its multi-mission timeseries, is applied. Similar to the single-mission filtering, any elevation change (dh) datapoint with more than 3 standard deviations from the regression model fit is removed. This changed filtering is more responsive to regional conditions, and consequently constrains the central ice sheet, where changes are smaller, more tightly than before, while allowing for more rapid ice movement at the ice sheet fringes.

1.3.4. Surface elevation change rate processing

Once a cross-calibrated multi-mission dh timeseries is made for each grid cell, the surface elevation change rate is calculated in a 5 year window that progresses in monthly steps. The datapoints within the cell/time period have their cross-calibrations applied if necessary - this is not always the case, as in some periods data from only one mission is available. If cross-calibration is necessary but could not be calculated, or if there are less than 10 datapoints in the timeseries, or if the data within the time window covers less than a 3 year period, then the change rate is not calculated. Otherwise, a linear least squares model is fitted to the timeseries and its gradient taken as the change rate.

The surface elevation change rate calculations are repeated at basin level, with caveats as discussed for cross-calibration (see Section 1.3.3 above). This is done to produce a basin level accuracy indicator (see Section 1.3.5 below), as users may wish to combine data from specific regions. Users should be aware that the ice velocity-guided method used may not be the ideal interpolation method over their specific region of interest. Proper scientific analysis should be tailored to the geographical regions and time periods involved. The dh timeseries for the grid cells are not explicitly provided in the current dataset but can be retrieved from the provided dh/dt timeseries.

1.3.5. Uncertainty processing

There are three contributors to the uncertainty of the surface elevation change rate

• input data
• cross-calibration
• modelling

The input data contribution depends on the distribution of elevation measurements within each grid cell/cycle. The standard deviation of these measurements does not formally account for all uncertainty sources, but will include residual errors from radar penetration and volume scattering that are not removed by retracking and backscatter power correction (see Section 1.3.2) and factors such as radar speckle, satellite location uncertainty and atmospheric attenuation uncertainty which decorrelate within the cycle period, see Wingham et al, 1998. When calculating rates of surface elevation change, the input component is formed from the individual errors on each datapoint used from the surface elevation change timeseries. It is taken as the root mean square of the input standard deviation, divided by the length of the time window.

The cross-calibration contribution accounts for errors in the biases calculated between each pair of missions, see Section 1.3.3. In any 5 year surface elevation change rate period, data from one or more radar altimeters may be used.

If only one mission is used, then no cross-calibration is necessary, and the cross-calibration uncertainty contribution is zero.

If two missions are used then the cross-calibration uncertainty between the two missions is converted to an uncertainty on the rate of change, by dividing by the time period over which the rate is calculated, which in this dataset is always 5 years.

If more than two missions are used then the root mean square of the cross-calibration uncertainties between each consecutive pair of missions is calculated and divided by the time period over which the rate is calculated, similarly to the two-mission case discussed earlier.

The modelling contribution is the standard deviation of the model fit. This is also the measure of stability.

The three contributions are summed in quadrature to give a total uncertainty, and this total uncertainty is provided in the product.

To provide one of the product's key performance indicators, basin accuracy, the uncertainty calculations are repeated at basin level, with caveats as discussed in the cross-calibration and surface elevation change rate processing sections (1.3.3 and 1.3.4) above. Results are reported in section 2.1 of related document D1, the Product Quality Assessment Report.

1.4. Output data

The output data CDR is contained in a netCDF4 file. After the initial CDR is made, monthly updated interim CDRs (ICDRs) will be issued, each containing the whole dataset from the previous release as well as the new additions, i.e. the product outputs are accumulative.

The main variable, the surface elevation change rate, is stored in a stack of EPSG:3031 polar stereographic grids (standard for Antarctica) at monthly intervals. Each grid cell contains the calculated surface elevation change rate for the five-year period centred on the time given for that grid, which is in a separate time array. A corresponding status array flags whether valid data exists in a grid cell or not. The grid projection coordinates of each cell and its longitude and latitude equivalents are given. Two masks are included, one of surface type (i.e. ice sheet, ice shelf, ice rise or island, and ocean), and one of slope levels (i.e., 0° to 2°, 2° to 5°, above 5°). More detail is given in Section 1.3.1 of the related document, the Product User Guide [D2]. The product fields are summarised in Table 8 below.

Table 8: Antarctic SEC output product variables summary

2. Surface elevation change, Greenland

The theoretical basis for the Greenland ice sheet SEC follows the R&D initiated in the ESA CCI project¹¹ and is detailed in the published ATBD¹². All applied algorithms are published in the relevant literature (Sørensen et al., 2015; Levinsen et al., 2015; Simonsen and Sørensen, 2017; Sørensen et al., 2018). The round-robin exercise completed during the ESA CCI project concluded the Greenland ice sheet surface elevation change be best represented by a combination of repeat-track, along-track, cross-over or plane-fitting algorithms (Levinsen et al., 2015). As seen in Section 1.3, the Antarctic ice sheet surface elevation change is only derived by the cross-over algorithm. The Greenland implementation of all four algorithms follows tightly Sørensen et al., (2018), in which the applied framework is presented. In the following, we present a summary of Sørensen et al., (2018). As of CDRv3, all the FORTRAN scripts described in the literature have been ported to Python, which allows for a true seamless processing chain and simplifies the generation of products. Similarly, to the Antarctic SEC the development in v4 is the consolidation of the CryoSat-2 processing baseline E.

2.1. Instruments

The Greenland surface elevation change product generation uses the radar altimetry record from the same series of overlapping ESA radar altimetry missions as for the Antarctic product. Section 1.1 describes the instruments used from each mission.

Table 1 in Section 1.1 lists the missions/mission phases used as input for the surface elevation change product generation in both hemispheres, and instrument specifications are given in Table 2 to Table 7. The only difference is the choice of operational modes for CryoSat-2 which for the Greenland ice sheet is limited to LRM, and SARIn, in contrast to the Antarctic ice sheet.

2.2. Input and auxiliary data

2.2.1. ESA Radar altimetry level 2 products

The input datasets for the Greenland surface elevation change product are the same as those for the Antarctic product. See Section 1.2.1 for details.

2.2.2. Auxiliary data

2.2.2.1. Digital elevation model for Greenland

The Greenland surface elevation change uses the official level-2 data products provided by ESA for all missions, as seen above (Section 2.2.1). The generation of such level-2 product includes the ESA processing facilities use of a digital elevation model (DEM) in the geolocation of LRM data. As the DEMs used might be satellite-mission dependent, we here refer to the individual mission documentation for the specific DEM used in the geolocation of the radar echo.

2.2.2.2. Ice extent

The processing is done for all Greenlandic grid-cells defined within the Greenland ice sheet or ice bodies with a strong connection to the ice sheet, by the ESA CCI glaciers project (Raster et al. 2012).

2.2.2.3. Glacial isostatic adjustment

No glacial isostatic adjustment is applied to the dataset, due to the large discrepancy in the model GIA signal in Greenland, and the limited bias in the resulting SEC (Barletta et al., 2013).

2.2.3. Validation data

The input datasets for the Greenland surface elevation change product are the same as those for the Antarctic product. See Section 1.2.3 for details.

2.3. Algorithms

The theoretical basis for the Greenland ice sheet SEC follows the R&D initiated in the ESA CCI project and can be found in the relevant literature (Sørensen et al., 2015; Levinsen et al., 2015; Simonsen and Sørensen, 2017; Sørensen et al., 2018). Here, we give an overview of the applied algorithms with special emphasis on the merging of data from different satellite missions.

2.3.1. Introduction

The product requirements are a stack of gridded surface elevation change rates from the Greenland ice sheet, at 25 km resolution, from the start of the ERS1 mission to present. The grids are given at monthly intervals, and flag grids for steep terrain and terrain type are also provided. The change rate unit is m/year. The GCOS accuracy target is 0.1 m/year, and the stability target is also 0.1 m/year. All products are provided on the equidistant grid in polar stereographic projection defined by EPSG: 3413 (see Product User Guide [D2] Section 2.1 for full details), and a timeseries of surface elevation change is derived in each cell for each mission. The processing chain is described in the related System Quality Assurance document [D3].

A review of the scientific background of surface elevation measurement for the Greenland ice sheet can be found in the ATBD¹³ provided by the ESA CCI project. Here, a combination of repeat-track, along-track, cross-over or plane-fitting algorithms were found to be the most optimal for the Greenland ice sheet. The mission-specific application of repeat-track, along-track, cross-over or plane-fitting algorithms is described in Sørensen et al. (2018). All algorithms for the Greenland ice sheet, are implemented in accordance with Sørensen et al (2018). In general terms, the input data are used to minimise the residual fit to the surface elevation change (dH/dt) model. The differences between the four algorithms are illustrated in Figure 2.1 and summarised in the following subsections.

Figure 2.1: Schematics of the algorithm differences, here the elevation change is estimated in the cyan hexagons. a) Cross-over algorithm. b) Repeat track, when all observations are along the same exact track (blue), and along-track, when the tracks differ in position. In this case, the “green” observations are projected onto the blue reference track. c) Plane-fitting, here the solution is modelled within regular-shaped grid cells. The figure is modified from Moholdt et al., (2010)

2.3.1.1. Cross-over

See Section 1.3.2 for a detailed description of the cross-over algorithm.

2.3.1.2. Repeat-track

The repeat-track (RT) algorithm is used to minimise the residual fit to the surface elevation change (dH/dt) model along segments of the repeated ground track. Included in the model are different biases for different parameters, which have been shown to correlate to surface elevation change, e.g. backscatter (Bs) and seasonality in the altimetry data. According to Sørensen et al (2018), here we follow the formulation of Legresy et al., 2005; Sørensen et al., 2011, 2015; Flament and Rémy, 2012, who model the time evolving elevation of the ice sheet given by

2.3.1.3. Along-track

The along-track algorithm follows the RT solution, but the formulation defines a reference track instead of segments on a repeated track. Data is then grouped and referenced to the reference track by a spatial search and all data within 2 km of the reference track is used. This formulation of the RT allows for the incorporation of multiple satellite missions and satellite mission with changes in the orbit configuration.

2.3.1.4. Plane-fitting

The plane-fitting (PF) method proposed by Simonsen and Sørensen (2017), has its foundation in the RT-algorithm presented above and in Sørensen et al. (2015). However, instead of solving for elevation change along-track, the drifting orbit of CryoSat-2 is utilised to solve in the full plane spanned by the applied grid. To account for the drifting orbit, Simonsen and Sørensen, (2017) proposed to add parameters, to solve for the 2d-topography within regular grid-cells (adding the parameters cx, cy, cc), to the equation presented above for the RT-algorithm. This results in an elevation change model-fitting of,

the naming convention follows eq. 2.1, with the addition of biases for CryoSat-2 ascending/descending orbits (b_AD) and LRM/SARIN mode (b_m) were added (AD and m is boolians), alongside with replacing the backscatter bias with biases due to changing leading edge width of the recorded waveform (LeW).

2.3.1.5. Adapting the plane-fitting algorithm for Sentinel-3

Sentinel-3 data were processed according to the plane-fitting algorithm applied to Cryosat-2. However, the experiments done in Simonsen and Sørensen, (2017) had to be redone in order to ensure the optimal fitting algorithm for the data from the new SAR sensor onboard Sentinel-3. Starting from the basic plane-fitting algorithm:

including a bias term for backscatter (Bs) correction, as the initial model guess, we performed 32 model perturbation experiments. Each of the 32-surface elevation change solutions were validated against Operation IceBridge to find the optimal choice of model-parameters in the plane-fitting algorithm. This model exercise showed that the most optimal surface elevation change solution was found by adding a LeW bias. This was also the case for the CryoSat-2 plane-fitting algorithm. 

2.3.1.6. Merging of satellite missions.

The round robin performed as a part of the ESA CCI Greenland ice sheet project showed collocation as the optimal interpolation method for combining heterogeneous data of different kinds. Since CDRv2 surface elevation change processor the collocation-gridding procedure has been replaced by an ordinary kriging method to provide model estimates for all Greenlandic ice-covered grid-cells. As ordinary kriging is capable of extrapolation over unrealistic distances, a distance flag has also been added to the product. This is given to enable the end-user to perform filtering of undesired data if needed.

With the CDRv3 surface elevation change processor upgrade, the cross-calibration of satellite missions went from being a part of the merging algorithm of surface elevation change grids to be included in the all-python processing scheme. This is illustrated by the updated RT-algorithm, given by

where dh_AB is the elevation bias between satellite mission A and B, given by the satellite mode parameter m_AB (0 for satellite A and 1 for satellite B). The LeW-bias is applied for each of the satellites A and B. For the Greenland ice sheet, this inter-mission implementation of the RT- and PF-algorithms derived at 1x1 km grid resolution, from 3- or 5-years of satellite altimetric observations. A new estimate of SEC (running average) is derived every 3^rd month, based on the optimal combination of cross-over, along-track, and plane-fitting method.

2.3.1.7. Deriving the monthly time series.

As the position of the data segments ( \( \overline{x}, \overline{y} \) ) varies from data-window to data-window, we need to combine multiple elevation change estimates to derive the Greenland-wide (gridded) estimates of surface elevation change, as illustrated in Figure 2.2. This is done for every month m_m by choosing the two elevation change estimates with the closest data-window mid-point to m_m. Ordinary kriging is then applied to the resulting point-cloud of dh/dt estimates to obtain the final 25x25-km gridded C3S-product.

Figure 2.2: Schematic drawing of the combined product generation. The upper panel shows the solutions available at a selected data segment of 1x1 km and 3- or 5-years of satellite altimetric observations. Whereas the lower panel represents the merging of the raw 1x1 km dh estimates to the final 25x25km C3S product. 

2.3.2. Uncertainty processing

There are two main contributors to the uncertainty in the surface elevation change rate:

• Measurement errors in the input data and spatial distribution.
• Fitting errors in the surface elevation change modelling.

The measurement error introduced by the input data depends on the spatial distribution of elevation measurements within each grid cell/cycle. Additional measurement-errors include radar penetration, volume scattering, radar speckle, satellite location uncertainty and atmospheric attenuation uncertainty which decorrelate within the cycle period and geolocation of the echo (Wingham et al. 1998). For the Greenland ice sheet, the geolocation of the echo is the largest error-source, as the nature of the radar-altimeter restrict the measurements to the "highest" point within the footprint of the satellite (Sørensen et al. 2018b). Hence, the surface elevation change estimate is a determination of the time-evolution of the highest points within the radar footprint. This potentially biases the derived solution toward more positive values (less elevation loss). It also hampers the solutions in valleys especially evident in the area of Jakobshavns Isbræ (see Figure 2.3), where the fixed grid solution (PF) struggles to provide accurate estimates of surface elevation in the main trough. In general, the radar altimeter performs better in the simple topography of the central, flat areas of the Greenland ice sheet compared to the coastal areas characterised by a more complex topography. Therefore, the error is generally larger in areas with steeper surface slopes.

Figure 2.3: Surface slope of the Greenland ice sheet (Simonsen and Sørensen 2017)

The fitting error has a well determined contribution to the total error as this error is directly derived when the residuals in the surface elevation change modelling are minimised. As the internal elevation change grids are generated the resulting errors are independent and can be summed as independent errors by the root-mean-square error. This summed error estimate is given in the present product release.

The number of unknowns in the measurement errors limits the end-to-end uncertainty characterisation. The errors induced from the time varying penetration of the surface snow and the geo-location of the radar echo are especially difficult to quantify. Here, we rely on independent dataset to estimate this contribution to the error. Applying a similar approach as for the C3S validation effort, Simonsen and Sørensen (2017) estimated a bias of 9 cm per year for the entire Greenland ice sheet with lower bias in the interior (6 cm per year), a direct implication of changing scattering horizon observed by the radar altimeter, not seen by an airborne altimeter measuring the true ice sheet surface.

2.4. Output data

The output CDR is contained in a single netCDF4 file, which is updated monthly as ICDRs (intermediate CDR). Each of the ICDRs contains the whole dataset from the previous release as well as the new additions, i.e., the product outputs are accumulative. See Figure 2.4 for an example of long-term elevation change rates derived by averaging the elevation change product.

Figure 2.4: The C3S surface elevation product averaged from 1992-2009, at the native 25x25 km grid in EPSG 3413 projection.

The main variable, the surface elevation change rate, is stored in a stack of polar stereographic grids (25x25 km grid, EPSG:3413) at monthly intervals. Each grid cell contains the calculated surface elevation change rate for the five-year period centered on the time given for that grid. This time information is stored in a separate time array. A corresponding status array flags whether valid data exists in a grid cell or not. The grid projection coordinates of each cell and its longitude and latitude equivalents are given. Two masks are included, one of surface type (i.e. ice sheet, ice shelf, ice rise or island, and ocean) and one of slope levels (i.e., 0° to 2°, 2° to 5°, above 5°). More details are given in the Product User Guide [D2], Section 2.3.1, and summarised in Table 9 below.

Table 9: Greenland SEC primary output variables summary

References

ACRI-ST IPF Team (2020). Product Data Format Specification – SRAL/MWR Level 2 Land products. ESA document reference S3IPF.PDS.003.2. Available from https://sentinel.esa.int/documents/247904/2753172/Sentinel-3-Product-Data-Format-Specification-Level-2-Land [URL resource last accessed 26th April 2023]

Barletta, V. R., Sørensen, L. S., and Forsberg, R. (2013). Scatter of mass changes estimates at basin scale for Greenland and Antarctica. The Cryosphere, 7(5), 1411–1432.

Brockley, D. et al. (2017). REAPER: Reprocessing 12 Years of ERS-1 and ERS-2 Altimeters and Microwave Radiometer Data. IEEE TGRS, June 2017. DOI: 10.1109/TGRS.2017.2709343

CryoSat-2 Team (2019). CryoSat-2 Product Handbook, baseline D 1.1. ESA document reference C2-LI-ACS-ESL-5319. Available from

https://earth.esa.int/eogateway/documents/20142/37627/CryoSat-Baseline-D-Product-Handbook.pdf/c76df710-2a5c-c8b8-00c1-13c8db0e9f51 [URL resource last accessed 26th April 2023]

Femenias, P. (editor) (2018). Envisat Altimetry Level 2 Product Handbook. ESA document CLS - ESLF - 18 -0003. Available from https://earth.esa.int/handbooks/ra2-mwr/ [URL resource last accessed 26th April 2023]

Flament, T. and Remy, F. (2012). Antarctica volume change from 10 years of Envisat altimetry.

Conference paper, Geoscience and Remote Sensing Symposium (IGARSS), 2012 IEEE International. DOI: 10.1109/IGARSS.2012.6351149 [URL resource last accessed 26th April 2023]

Geminale, T. (2021). CRYOSAT Ground Segment Instrument Processing Facility Baseline E Evolutions, C2-RP-ACS-ESL-5330, Issue 2. Available from https://earth.esa.int/eogateway/documents/20142/37627/Cryosat-Baseline-E-Evolutions.pdf/bdd53679-1de5-ebd4-3d18-b72096f7b7c8 [URL resource last accessed 26th April 2023]

Ivins E. R., James T. S., Wahr J., Schrama O., Ernst J., Landerer F. W., and Simon K. M. (2013). Antarctic contribution to sea level rise observed by GRACE with improved GIA correction. J. Geophys. Res. Solid Earth, 118(6), 3126–3141.

Khvorostovsky, K. et al. (2018). Algorithm Theoretical Baseline Document (ATBD) for the Greenland Ice Sheet CCI Project of ESA’s Climate Change Initiative, version 3.2. ESA document ID: ST-DTU-ESA-GISCCI+-PMP001

Legresy, B. et al. (2005). Envisat radar altimeter measurements over continental surfaces and ice caps using the ICE-2 retracking algorithm. Remote Sensing of Environment, v95, p150-163.

Levinsen, J.F. et al. (2015). ESA ice sheet CCI: derivation of the optimal method for surface elevation change detection of the Greenland ice sheet – round robin results. International Journal of Remote Sensing, v36 (2), p551-573. DOI: 10.1080/01431161.2014.999385

Moholdt, G., Hagen, J. O., Eiken, T. and Schuler, T. V. (2010) Geometric changes and mass balance of the Austfonna ice cap, Svalbard. The Cryosphere, 4, 21-34. DOI:10.5194/tc-4-21-2010.

Nagler, T. et al. (2018). Algorithm Theoretical Baseline Document (ATBD) for the Antarctic Ice Sheet CCI Project of ESA’s Climate Change Initiative, version 3.0. Available from https://climate.esa.int/en/projects/ice-sheets-antarctic/key-documents/ [URL resource last accessed 26th April 2023]

Padman, L. et al. (2002). A new tide model for the Antarctic ice shelves and seas. Annals of Glaciology v34 pp247-254. DOI: 10.3189/172756402781817752

REAPER Team. (2014). REAPER product handbook. Available from https://earth.esa.int/eogateway/documents/20142/37627/reaper-product-handbook-for-ers-altimetry-reprocessed-products.pdf [URL resource last accessed 26th April 2023]

Rignot, E., Mouginot, J. and Scheuch, B. (2011). Ice Flow of the AIS. Science 333 (6048):1427-1430

Shepherd, A., Ivins, E. R., Geruo, A., Barletta, V. R., Bentley, M. J., Bettadpur, S., Briggs, K. H., Bromwich, D. H., Forsberg, R., Galin, N., Horwath, M., Jacobs, S., Joughin, I., King, M. A., Lenaerts, J. T. M., Li, J. L., Ligtenberg, S. R. M., Luckman, A., Luthcke, S. B., McMillan, M., Meister, R., Milne, G., Mouginot, J., Muir, A., Nicolas, J. P., Paden, J., Payne, A. J., Pritchard, H., Rignot, E., Rott, H., Sorensen, L. S., Scambos, T. A., Scheuchl, B., Schrama, E. J. O., Smith, B., Sundal, A. V., van Angelen, J. H., van de Berg, W. J., van den Broeke, M. R., Vaughan, D. G., Velicogna, I., Wahr, J., Whitehouse, P. L., Wingham, D. J., Yi, D. H., Young, D., and Zwally, H. J. (2012). A Reconciled Estimate of Ice-Sheet Mass Balance, Science, 338, 1183–1189,

Simonsen, S.B., and Sørensen, L.S., (2017). Implications of changing scattering properties on Greenland Ice Sheet volume change from CryoSat-2 altimetry. Remote Sensing of Environment, v190, p 207-216.

Slater, T., Shepherd, A., McMillan, M., Muir, A., Gilbert, L., Hogg, A. E., Konrad, H., and Parrinello, T. (2018). A new digital elevation model of Antarctica derived from CryoSat-2 altimetry. The Cryosphere, 12, 1551-1562. DOI:10.5194/tc-12-1551-2018.

Sørensen, L.S., et al. (2011). Mass balance of the Greenland ice sheet (2003-2008) from ICESat data – the impact of interpolation, sampling and firn density. The Cryosphere, v5, p173-186. DOI: 10.5194/tc-5-173-2011

Sørensen, L.S., et al. (2015). Envisat-derived elevation changes of the Greenland Ice Sheet, and a comparison with ICESat results in the accumulation area. Remote Sensing of Environment, v160, p56-62. DOI: 10.1016/j.rse.2014.12.022

Sørensen, L.S., et al. (2018). 25 years of elevation changes of the Greenland Ice Sheet from ERS, Envisat and CryoSat-2 radar altimetry. Earth and Planetary Science Letters, vol 495, p 234-241. DOI: 10.1016/j.epsl.2018.05.015

Sørensen, L.S, Simonsen, S. B., Langley, K., Gray, L., Helm, V., Nilsson, J., Stensengm, L., Skourup, H., Forsberg, R. and Davidson, M. W. J. (2018b). Validation of CryoSat-2 SARIn Data over Austfonna Ice Cap Using Airborne Laser Scanner Measurements. Remote Sensing, 10(9), 1354. DOI:10.3390/rs10091354

Studinger, M. (2014). IceBridge ATM L4 Surface Elevation Rate of Change, Version 299 1, Antarctica subset. N. S. a. I. D. C. D. A. A. Center. Boulder, Colorado, USA. DOI: 10.5067/BCW6CI3TXOCY

Tabachnick, B.G. and Fidell, L.S. (2019). Using Multivariate Statistics, seventh edition. Published by Pearson, ISBN-13: 978-0-13-479054-1

Wingham, D.J., Ridout, A.J., Scharroo, R., Arthern, R.J., and Shum, C.K. (1998). Antarctic Elevation Change from 1992 to 1996. Science vol 282, issue 5388, pp 456-458, DOI:10.1126/science.282.5388.456

Zwally, H. J., Giovinetto, M.B., Beckley, M.A. and Saba, J.L. (2012a), Antarctic and Greenland Drainage Systems, GSFC Cryospheric Sciences Laboratory. Available at https://earth.gsfc.nasa.gov/cryo/data/polar-altimetry/antarctic-and-greenland-drainage-systems [URL resource last accessed 26th April 2023]

Zwally, H. J., Brenner, A.C., Major, J.A., Bindschadler, R.A. and Marsh, J.G. (2012b). Growth of Greenland Ice Sheet. Measurement, Science vol 246, issue 4937, pp, 1587-1589. DOI:10.1126/science.246.4937.1587

Related articles

• Page:
  ERB TSI: Product User Guide and Specification (PUGS)
• Page:
  ERB TSI: Product Quality Assurance Document (PQAD)
• Page:
  ERB TSI: Product Quality Assessment Report (PQAR)
• Page:
  ERB TSI: System Quality Assurance Document (SQAD)
• Page:
  ERB TSI: Algorithm Theoretical Basis Document (ATBD)