Last modified on Jun 27, 2024 17:17

Contributors:  Lin Gilbert (University Leeds), Sebastian Bjerregaard Simonsen (Technical University of Denmark), Jan Wuite (ENVEO)

Issued by: University of Leeds / Lin Gilbert

Issued Date: 04/05/2021

Ref:  C3S_312b_Lot4.D1.IS.6_v3.0_SEC_Algorithm_Theoretical_Basis_Document_i1.0

Official reference number service contract:  2018/C3S_312b_Lot4_EODC/SC2 

History of modifications

List of datasets covered by this document

Related documents

Acronyms

Scope of the document

This document is the Algorithm Theoretical Basis Document for Surface Elevation Change (SEC) as part of the Copernicus Ice Sheets and Ice Shelves service. It describes the satellite-mounted instruments, auxiliary datasets, auxiliary models and basic algorithms used to create the data products.

Executive summary

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

• Ice velocity is given for Greenland in product D3.IS.4
• Gravimetric mass balance is given for Greenland and Antarctica in product D3.IS.5
• Surface elevation change is given for:
  • Antarctica in product D3.IS.6.1
  • Greenland in product D3.IS.6.2

We document here their inputs, processing stages and outputs for the production of the two Polar region SEC products.

1. Surface elevation change, Antarctica – D3.IS.6.1

1.1. Instruments

In both hemispheres the SEC product uses the altimetry record from a series of overlapping ESA radar altimetry missions. The record started in 1992 with the ERS-1, ERS-2 and EnviSat satellites, 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 tacking by the radar altimeter as the satellite approaches the complex topography of the ice sheets from the oceans. When a dedicated land-ice processing is adapted for the Sentinel-3 mission, this will be the data product used for the surface elevation change processing.

In the v3 dataset two evolutions from v2 have taken place. The CryoSat-2 baseline used has been updated from C to D, since the full mission reprocessing was made available in the middle of 2020. These two baselines differ in several respects, notably the slope mask used in low resolution mode. Sentinel-3A data was included in the v2 CDR, and Sentinel-3B data is added in v3. Reprocessed releases for both Sentinels are expected from ESA, to resolve problems in transitioning from ocean to land ice surfaces, but unfortunately will not be released in time to be incorporated in the v3 product.

The Antarctic region surface elevation change dataset is derived from the crossover method (see Section 1.3 below and Zwally, Brenner et al, 2012). 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: Mission summary

Summary specifications for each instrument are provided in Table 2 to Table 6. Only operating modes that are used in the Antarctic and/or Greenland SEC products are listed.

Table 2: ERS1 and ERS2 RA

Table 3: EnviSat RA-2

Table 4: CryoSat-2 SIRAL, LRM (low resolution mode)

Table 5: CryoSat-2 SIRAL, SARIn (synthetic aperture radar interferometry) mode

Table 6: Sentinel-3A and B SRAL, in main mission will always run in SAR mode

For further information, these missions and their instruments are described by ESA's Earth Observation portal, with a top-level directory at:

• https://directory.eoportal.org/web/eoportal/satellite-missions

Radar altimetry is described in detail in a PDF tutorial document at:

• www.altimetry.info/filestorage/Radar_Altimetry_Tutorial.pdf

1.2. Input and auxiliary data

1.2.1. ESA Radar altimetry level 2 products

The initial main input data comes from four sources:

• ERS1 and 2 Reprocessing of Altimeter Products for ERS (Reaper) RA L2, see Brockley et al 2017.
• Envisat RA-2 L2 Geophysical Data Records GDR_v3, see Femenias (ed), 2018
• CryoSat-2 SIRAL L2i Low Resolution Mode (LRM) and SAR-Interferometric (SIN) datasets, see CryoSat-2 Product Handbook, 2019
• Sentinel-3A/B L2 SRAL, specifically the SR_2_LAN_NT products containing land ice data, see ACRI-ST IPF Team, 2018

Data is freely available on registration from these websites:

• ERS1 and ERS2: https://earth.esa.int/web/sppa/mission-performance/esa-special-series/reprocessing-of-altimeter-products-for-ers
• Envisat: https://earth.esa.int/web/guest/-/ra-2-geophysical-data-record-1470
• CryoSat-2: https://earth.esa.int/web/guest/-/how-to-access-cryosat-data-6842
• Sentinel-3A/B: https://scihub.copernicus.eu

Processing note – Sentinel-3 operational data processing was switched from baseline 3 to baseline 4 in spring 2020, without a full mission reprocessing. The change in baselines introduced a bias in the backscatter power, which is used in crossover processing. This bias must be removed when using data from a time range that includes both baselines. In practical terms, for Antarctic SEC, this is implemented by the Leeds data server prior to the service pipeline processing – the pipeline requests backscatter data from baseline 4 that has been rebiased to match baseline 3.

1.2.2. Auxiliary data

1.2.2.1. Surface type mask

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, Giovinetto et al, 2012. Further information and a freely downloadable copy of the mask are available from:

• https://icesat4.gsfc.nasa.gov/cryo_data/ant_grn_drainage_systems.php

1.2.2.2. Slope model

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:

• https://earth.esa.int/web/guest/missions/esa-operational-eo-missions/cryosat

1.2.2.3. Drainage basin mask

The Ice sheet Mass Balance Inter-comparison Exercise (IMBIE) 2 drainage basin outlines are used to split the Antarctic land mass into 27 basins. See Shepherd et al, 2012. These are described in Zwally, Giovinetto et al, 2012, and freely available from:

• http://icesat4.gsfc.nasa.gov/cryo_data/ant_grn_drainage_systems.php

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. The 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 freely available as an ancillary dataset from the IMBIE project, at:

• http://imbie.org/data-downloads/

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. In 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. AIS ice velocity map

During validation a bias value is calculated to compensate for the tendency of the ATM to preferentially sample fast-thinning ice. This bias is calculated from the 1km by 1km Berkeley – Ice Sheet Initiative for Climate Extremes (BISICLES) ice velocity map, see Rignot et al, 2011. It is freely available from:

• https://commons.lbl.gov/display/bisicles/BISICLES.

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 Airborne Topographic Mapper, a scanning laser altimeter flown on board aircraft by Operation IceBridge (Studinger 2014).

The validating dataset used for both hemisphere's 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, from;

• https://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 25km 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 a polar stereographic projection (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 Antarctic Ice Sheet Climate Change Initiative (AIS 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.

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 good geographical coverage. Each mission is treated individually. Figure 1 shows the physical layout of a single crossover.

Figure 1: 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 H₁ and H₂ 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, 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.

To calculate the heights, slope-corrected ranges are used, all geophysical corrections are made as applicable, and all ocean tides corrections, calculated with the CATS2008a model, are made.

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), and where the numbers of measurements from ascending and descending passes are greater than a given ratio (0.5).

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 dp timeseries of power differences relative to the reference cycle power at each crossover 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 recently established Sentinel-3 A and B datastreams, a 2 year period has to be 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 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. Values are removed if they amount to more than the annual change cycle magnitude plus the maximum change per annum for every year of time difference to the reference cycle.

Timeseries are removed if they contain 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, t is 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.

Finally, timeseries in any cell with a surface slope of more than 5° are removed, and elevation changes (dh) of more than 10m are removed as this has been found to be physically unlikely.

1.3.3. Multi-mission cross-calibration

The v3 cross-calibration method is the same as that for CDR v2. 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 2 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 2: 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.

Mission datapoints are not inferred from the model. Biases are added to pre-existing datapoints only.

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 optimum 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.

For test purposes, 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.

1.4. Output data

The output data is contained in a netCDF4 file, called a climate data record (CDR). After the initial CDR is made, monthly updated ICDRs (interim CDR) 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 polar stereographic grids at monthly intervals. Each grid cell contains the calculated surface elevation change rate for the five-year period cantered 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].

2. Surface elevation change, Greenland – D3.IS.6.2

The theoretical basis for the Greenland ice sheet SEC follows the R&D initiated in the ESA CCI project  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 to 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 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). With the upgrade to version 3 (CDRv3), all the FORTRAN scripts described in the literature have been ported to Python. This allows for a true seamless processing chain and simplifies the generation of CDRv3 products.

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 6. 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 25km 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.

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 3 and summarised in the following subsections.

Figure 3: 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. 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

With the upgrade to version 3, the Bs-correction (eq. x) utilized for ERS-1, ERS-2 and Envisat was replaced by a correction to account for changes in the leading-edge width of the recorded waveform (LeW) following our experience of implementing the plane-fitting algorithm for CryoSat-2 and Sentinel-3 (see Section 2.3.1.4).

2.3.1.2. 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.3. Cross-over

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

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, resulting in an elevation change model-fitting of,

where also the additional biases for CryoSat-2 (b_AD ascending/descending and b_m LRM/SARIN mode) were added, 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-3A

Sentinel-3 data were processed according to the plane-fitting algorithm applied to Cryosat-2. However, the experiments done in Simonsen et al. 2017 had to be redone in order to insure the optimal fitting algorithm for 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. For the CDRv2 and CDRv3 surface elevation change processor the collocation-gridding procedure was 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 is used to fit the surface of 1x1 km-segments including 3- or 5-years of satellite altimetric observations. The fit is performed every 3rd month. For each of these data segments the optimal parameters are used to drive the surface elevation (H(x, y,tm)) at the midpoint of the data segments (x, y) and at monthly (tm) temporal resolution for months inside the temporal data-window, as illustrated in the upper diagram of Figure 4.

2.3.1.7. Deriving the monthly time series.

As the position of the data segments (x, 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. 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 4: 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 1x1km 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-error 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 (lower elevation loss). It also hampers the solutions in valleys especially evident in the area of Jakobshavns Isbræ, 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.

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).

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.

The main variable, the surface elevation change rate, is stored in a stack of polar stereographic grids 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. 

References

ACRI-ST IPF Team, Product Data Format Specification – SRAL/MWR Level 2 Land products, ESA document reference S3IPF.PDS.003.2, accessed on 27/01/2020 from https://sentinel.esa.int/documents/247904/2753172/Sentinel-3-Product-Data-Format-Specification-Level-2-Land

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.

Brenner, A.C., J.P. DiMarzio, & H.J. Zwally, Precision and accuracy of satellite radar and laser altimeter data over the continental ice sheets, IEEE Trans. Geosci. Remote Sens. 45, 2, 321–331 (2007).

Brockley 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 Product Handbook, baseline D 1.1, C2-LI-ACS-ESL-5319, available from https://earth.esa.int/web/guest/document-library/browse-document-library/-/article/cryosat-baseline-d-ice-product-handbook, 2019

Femenias (editor), 2018, Envisat Altimetry Level 2 Product Handbook, ESA document CLS - ESLF - 18 -0003, available from https://earth.esa.int/web/guest/missions/esa-operational-eo-missions/envisat/news/-/article/envisat-altimetry-v3-0-full-mission-reprocessing

Flament, T. and F. Remy (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

Howat, I.M., Negrete, A., and Smith, B.E. (2014). The Greenland Ice Mapping Project (GIMP) land classification and surface elevation data sets. The Cryosphere, 8, 1509-1518. DOI: 10.5195/tc-8-1509-2014

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., 2012. Merging and analysis of elevation time series over Greenland Ice Sheet from satellite radar altimetry. IEEE Trans. Geosci. Remote Sens. 50 (1), 23-36. DOI: 10.1109/TGRS.2011.2160071

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

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.: Geometric changes and mass balance of the Austfonna ice cap, Svalbard, The Cryosphere, 4, 21-34, https://doi.org/10.5194/tc-4-21-2010, 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 http://www.esa-icesheets-antarctica-cci.org/

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 product handbook, 2014, in the ESA Earth Online document library at https://eart.esa.int/web/guest/document-library

Rignot, E., J., Mouginot and B. Scheuch, 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.: A Reconciled Estimate of Ice-Sheet Mass Balance, Science, 338, 1183–1189, 2012.

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, https://doi.org/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. https://doi.org/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., Using Multivariate Statistics, seventh edition, published by Pearson, 2019, 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. Jay, Mario B. Giovinetto, Matthew A. Beckley, and Jack L. Saba, 2012, Antarctic and Greenland Drainage Systems, GSFC Cryospheric Sciences Laboratory, at http://icesat4.gsfc.nasa.gov/cryo_data/ant_grn_drainage_systems.php

Zwally, H. Jay, Anita C. Brenner, Judy A. Major, Robert A. Bindschadler and James G. Marsh, 2012, Growth of Greenland Ice Sheet: Measurement, Science vol 246, issue 4937, pp, 1587-1589. doi:10.1126/science.246.4937.1587

3. Related articles

• Page:
  IV Version 1.2; SEC Version 2.0; GMB Version 2.0: Product User Guide Specification (PUGS)
• Page:
  IV Version 1.2; SEC Version 2.0; GMB Version 2.0: System Quality Assurance Document (SQAD)
• Page:
  IV Version 1.2; SEC Version 2.0; GMB Version 2.0: Product Quality Assurance Document (PQAD)
• Page:
  IV Version 1.2; SEC Version 2.0; GMB Version 2.0: Product Quality Assessment Report (PQAR)
• Page:
  IV Version 1.2; SEC Version 2.0; GMB Version 2.0: Algorithm Theoretical Basis Document (ATBD)