Contributors: Lin Gilbert (University of Leeds), Sebastian B. Simonsen (Technical University of Denmark), Jan Wuite (ENVEO IT GmbH)

Issued by: University of Leeds / Lin Gilbert

Date: 31/05/2020

Ref: C3S_312b_Lot4.D1.IS.2_v2.0_202001_Algorithm_Theoretical_Basis_Document_v1.0

Official reference number service contract: 2018/C3S_312b_Lot4_EODC/SC2

Note: This document provides the following three deliverables:

D1.IS.2-v2.0 Algorithm Theoretical Basis Document - Ice Velocity

D1.IS.4-v2.0 Algorithm Theoretical Basis Document – Gravimetric Mass Balance

D1.IS.6-v2.0 Algorithm Theoretical Basis Document – Surface Elevation Change (Antarctica) and Surface Elevation Change (Greenland)

Table of Contents

History of modifications

Issue

Date

Description of modification

Author

v0.1

13/01/2020

The present document was modified based on the document with deliverable ID: D1.IS.2_v1.0 (C3S_312b_Lot4.D1.IS.2_v1.0_Algorithm_Theoretical_Basis_v1.5)

CC

v1.0

31/05/2020

Updated dataset list and related documents to reference v2. Updated Section 3 with new cross-calibration method and addition of Sentinel-3A data. Section 4 (sections 4.1, 4.2 and 4.3.1.5) updated to include Sentinel-3A data. Executive summary modified to include Antarctic.

LG/JW/SS

List of datasets covered by this document

Deliverable ID

Product title

Product type (CDR, ICDR)

Version number

Delivery date

D3.IS.4

Ice velocity

CDR

2.0

31/01/2020

D3.IS.5

Gravimetric mass balance

CDR

2.0

31/01/2020

D3.IS.6.1

Surface elevation change, Antarctica

CDR

2.0

31/01/2020

D3.IS.6.2

Surface elevation change, Greenland

CDR

2.0

31/01/2020

Related documents

Reference ID

Document

D1.IS.2-v2.0

Development Milestone Plan and Status

D2.IS.2-v2.0

Product Quality Assessment Report

D3.IS.7-v2.0

Product User Guide

D1.IS.1-v2.0

System Quality Assurance Document

D2.IS.1-v2.0

Product Quality Assurance Document

Acronyms

Acronym

Definition

AIS

Antarctic Ice Sheet

ASAP

Austrian Space Applications Programmes

AT

Along-Track

ATM

Airborne Topographic Mapper

BISICLES

Berkeley – Ice Sheet Initiative for Climate Extremes

CATS

Circum-Antarctic Tidal Simulation

CCI

Climate Change Initiative

CDR

Climate Data Record

DEM

Digital Elevation Model

DLR

Deutsches zentrum fur Lüft- und Raumfahrt e.V. (German Aerospace Centre)

ERS

European Remote-sensing Satellite

ESA

European Space Agency

ESRI

Environmental Systems Research Institute

GDR

Geophysical Data Record

GIA

Glacial Isostatic Adjustment

GIMP

Greenland Ice sheet Mapping Project

GLAS

Geoscience Laser Altimeter System

GMB

Gravimetric Mass Balance

GrIS

Greenland Ice Sheet

ICDR

Interim Climate Data Record

IMBIE

Ice sheet Mass Balance Intercomparison Exercise

InSAR

Interferometric SAR

IV

Ice Velocity

IW

Interferometric Wide-swath mode

LRM

Low Resolution Mode

MODIS

Moderate Resolution Imaging Spectrometer

OT

Offset Tracking

POD

Precise Orbit Determination

POE

Precise Orbit Ephemerides

PROMICE

PROgram for Monitoring of the Greenland ICE sheet

RA

Radar Altimeter

REAPER

REprocessing of Altimeter Products for ERS

RT

Repeat Track

SAR

Synthetic Aperture Radar

SARIN

Synthetic Aperture Radar INterferometer

SEC

Surface Elevation Change

SIN

Synthetic aperture radar INterferometer

SIRAL

SAR/Interferometric Radar ALtimeter

SLC

Single Look Complex

S1

Sentinel-1

TOPS

Terrain Observation by Progressive Scans

XO

Cross-Over

Scope of the document

This document is the Algorithm Theoretical Basis Document for 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 describe their inputs, processing stages and outputs.

The gravimetric mass balance is brokered from the Greenland and Antarctic CCI projects. The CCI project's ATBD (Khvorostovsky et al. 2018; Nagler et al. 2018) is provided within Annex A. Only output specifics will be given in the main text, other topics will refer to Annex A.

1. Greenland Ice Sheet Velocity – D3.IS.4

1.1. Instruments

The primary satellite source for generation of the ice velocity (IV) fields is Copernicus Sentinel-1A and Sentinel-1B. The Sentinel-1 (S1) satellites carry a C-band synthetic aperture radar (SAR) instrument providing high resolution SAR images in different acquisition modes. The Interferometric Wide-Swath mode (IW) is the main acquisition mode over land areas, including Greenland. IW allows combining a large swath width (250 km) with a moderate geometric resolution (5 m x 20 m).

Sentinel-1A provides continuous coverage for the entire Greenland Ice Sheet margin since October 2014 with a 12-day repeat pass period. Sentinel-1B was launched in April 2016 and reduced the repeat pass period to 6 days. Besides the continuous coverage of the margin, the ice sheet is covered completely during dedicated annual (winter) campaigns, with 4 to 6 repeat observations per track. In 2019 additional tracks were added for continuous acquisition covering most of the interior of the ice sheet in ascending and descending direction, providing new opportunities for interferometric applications. The S1 acquisition plan for Greenland makes the constellation the primary source for year-round monitoring of IV. The main instrument parameters and IW mode are described in Table 1 and Table 2

Table 1: Sentinel-1 sensor specifications

Property

Value

Sensor

C-SAR (C-band Synthetic Aperture Radar)

Emitted frequency (GHz)

5.4

Pulse repetition frequency (Hz)

1 000 - 3 000

Pulse duration (microseconds)

5-100

Bandwidth (MHz)

0-100

Antenna length (m)

12.3

Antenna beamwidth (degrees)

0.23

Slant range resolution (m)

5

Azimuth resolution (m)

20

Table 2: Characteristics of Sentinel-1 IW mode

Parameter

Value

Swath width (ground range)

250 km

Nom. L1b product length

50 km

Full performance incidence angle range

20°- 46°

Data access incidence angle range

30° - 42°

Azimuth resolution

20 m

Ground range resolution

5 m

Polarizations

HH+HV, VV+VH, VV, HH

1.2. Input and auxiliary data

1.2.1. SAR data

Primary input data for IV retrieval are repeat pass Sentinel-1A and Sentinel-1B Level-1 Single Look Complex (SLC) images acquired in the Interferometric Wide (IW) swath mode. These products consist of focused SAR data geo-referenced using orbit and altitude data from the satellite and provided in zero-Doppler slant-range geometry. IW implements the Terrain Observation by Progressive Scans (TOPS) mode, consisting of 3 or 5 sub-swaths, each consisting of a series of bursts. Each burst has been processed as a separate SLC image. ESA provides Sentinel-1 data free of charge through the Copernicus Open Access Hub (SciHub) and various mirror sites.

1.2.2. Auxiliary data

Auxiliary data, used in the processing chain for IV production, include precise orbit files, a digital elevation model and ocean/land masks. Details are provided in the following:

1.2.2.1. Orbits

Orbit files (Precise Orbit Ephemerides - POE) are required for IV processing and enable accurate co-registration. Precise orbit files are available as separate products and are provided by the Precise Orbit Determination (POD) service for SENTINEL-1 approximately 3 weeks after acquisition. POE files cover approximately 28 hours and contain orbit state vectors at 10-second intervals.

1.2.2.2. DEM

A contemporary high-resolution digital elevation model (DEM) is required in the processing chain for various purposes, including co-registration, calculation of the vertical component of velocity (vz) and final geocoding of the product. The DEM used for the annual Greenland Ice Sheet velocity map is compiled from TanDEM-X 90m DEM tiles provided by DLR (Rizolli et al., 2017). This is a product variant of the original 12m (0.4 arcsec) DEM product with reduced pixel spacing (Wessel et al., 2018). The extent of the DEM is equal to the IV product. Further information can be found in the TanDEM-X Ground Segment DEM Products Specification Document (Wessel, 2018).

1.2.2.3. Land/Ocean Mask

This is an ESRI-shapefile of the ocean/land boundary, based on recent optical imagery and used for masking the retrievals in the ocean (usually coming from tracked sea ice). The shapefiles are internal datasets which are updated as required.

1.3. Algorithms

1.3.1. Coherent and Incoherent Offset Tracking using Sentinel-1

Offset tracking (OT) is the principal method for generating velocity fields on glacier surfaces from Sentinel-1 operationally. OT requires less operator interaction than interferometric techniques, i.e. can be more easily and reliably automated. Another advantage of the OT technique over InSAR is the provision of two components of the velocity vector (along track and line-of-sight, LOS) from data of a single track. Also, the Greenland Ice Sheet is subject to highly variable meteorological conditions, resulting in rapid temporal decorrelation of the radar signal phase impairing the comprehensive application of DInSAR. OT is less sensitive or insensitive to temporal decorrelation of the radar phase signal, albeit at a lower accuracy. In areas with distinct and stable surface features, as common on mountain glaciers and outlet glaciers of ice sheets, coherence of the repeat-pass SAR data is not required because the cross-correlation of image-templates is based on the amplitude signal. In the level interior sections of ice sheets, without distinct amplitude features, temporal stability of the target phase (coherence) is required for speckle tracking. Loss of stability of the target phase (temporal decorrelation) is less critical than for the InSAR method because a continuous path is not required for deriving velocities of an isolated patch in an image, as is the case with InSAR.

The ice velocity (IV) algorithm is developed as part of the ESA CCI and Austrian Space Applications (ASAP) Programmes. The IV system's core module performs coherent and incoherent offset tracking (OT), utilizing long stripes of Sentinel-1 data acquired in interferometric wide swath (IW) mode. OT refers to several related methods that include speckle tracking, coherence tracking and amplitude tracking or feature tracking. Feature tracking uses cross-correlation of image patches to find the displacement of surface features such as crevasses or rifts and edges, that move with the same speed as the ice and are identifiable on two co-registered amplitude images, and subsequently ice flow velocity. In coherence tracking, the offset which maximises the interferometric coherence within a certain window size is determined and used to derive the ice velocity.

For Sentinel-1 Terrain Observation by Progressive Scans (TOPS) mode acquisitions, OT algorithms need to support burst handling and stitching. Also, the azimuth-varying phase introduced by the TOPS beam steering must be accounted for, both for complex offset tracking and for incoherent offset-tracking, where a factor of two oversampling of the SLCs is required prior to detection.

OT accuracies depend on the level of coherence and on the correlation window size (Bamler and Eineder, 2005) as well as on the type of features that are being tracked. Ionospheric scintillations can have a large impact on offset-tracking by causing azimuth shifts that are introduced by the fluctuating electron density along the sensor path. Misregistration is not a problem, but the large azimuth shifts are interpreted as ice motion and observed offsets can exceed the azimuth pixel size. DEM errors have less impact on the algorithm performance compared to DInSAR, since there is no phase unwrapping. However, since offset-tracking can be applied with longer temporal baselines, and even though larger spatial baselines can be used for offset-tracking, in the typical scenario the sensitivity to DEM errors is much smaller.

For further details the reader is referred to the Algorithm Theoretical Basis Documents of the Greenland and Antarctic Ice Sheet CCI projects (Khvorostovsky et al. 2018; Nagler et al. 2018) and Nagler et al. (2015).

1.3.2. Generation of Regional Ice Velocity Maps from Multiple Tracks using OT

The generation of a regional ice velocity map requires the combination of offset tracking results, from several tracks, each providing slant range and azimuth displacements at local incidence and heading angle. Figure 1 illustrates the high-level processing line for IV production. The system includes 3 main modules, which are detailed in the System Quality Assurance Document, but briefly explained here.

Figure 1: High-level flow chart of the IV processing system. Green – input data, Blue – processing modules, Red - product and intermediate products, Yellow – product data base.

Within the IV module SAR data and orbit data are imported into the system and velocity maps are generated for pairs of repeat pass data of the same track. The long tracks are geocoded into the common map projection of the output grid, providing observations of slant range displacement and azimuth for each track separately. The local incidence and heading angle are calculated using the annotated image and orbit parameters and a DEM as input.

The Merge module combines all IV products from all tracks and image pairs over a specified time span (i.e. 1 year), applying a weighted average approach. The combination of displacement observations of multiple tracks applies a weighted least squares model to fit the multiple observed displacements according to:

$$y = Ax + \epsilon$$

where x is the horizontal velocity vector (in local Easting/Northing coordinates), y is a vector with the observed velocities (in slant range/azimuth), A is the matrix projecting the horizontal velocity to slant range/azimuth geometry, and is a noise vector. The matrices are defined as:

This equation is solved for each output pixel separately: n is the number of pairs with a valid slant range and azimuth displacement estimate at the current output pixel position, $(\Delta sr_n, \Delta a z_n)$ is the measured slant range/azimuth velocity from pair number n, $(\theta_n, \phi_n)$ is the elevation and azimuth angle (the latter measured counter-clockwise from East) of the look vector from the pixel to the sensor for pair n, $(\frac{\partial z}{\partial x}, \frac{\partial z}{\partial y})$ is the local height gradient at the pixel (from the DEM), and $(\nu_E, \nu_N)$ is the horizontal ice velocity to be calculated. The noise vector,$\epsilon$, is assumed normal distributed with zero mean and diagonal covariance matrix $\Sigma$:


where $(\sigma_{sr_i}, \sigma_{az_i})$ corresponds to the weighted pixel spacing in slant range and azimuth direction and is calculated according to: $$\sigma^2_{sr} = NCC \cos(\theta)/p_{sz}$$ $$\sigma^2_{sr} = NCC /p_{az}$$ Where $p_{sz}$ and $p_{az}$ correspond to the slant range and azimuth pixel spacing and NCC is the Normalised Cross-Correlation Coefficient: $$NCC(u,v)=\frac{\sum \left[f(x,y)-\overline{f}_{u,v} \right] [t(x-u,y-v)-\overline{t}]}{\sqrt{\sum \left[ f(x,y) - \overline{f}_{u,v} \right]^2 \sum [t(x-u,y-v)- \overline{t}]^2}}$$

where f() is the master, t() is the template image at the pixel position (u,v). The peak in the correlation matrix corresponds to the best matching shift of the slave window with respect to the master window.
The weighted least squares estimate of the horizontal velocity vector is then:

$$\hat{x} = (A^T \Sigma^{-1}A)A^T \Sigma^{-1} y$$

An advantage of the approach above is that if crossing tracks (ascending/descending) are used, the displacement measurement in the direction with higher resolution dominates in the velocity retrieval. On the other hand, it allows to reduce e.g. ionospheric azimuthal artefacts by adjusting the weights for azimuth shifts. In the case of only one track (single range/azimuth offset-tracking pair), the proposed method provides an exact inversion, using the assumption of surface-parallel flow by incorporating the local height gradient derived from the DEM. The standard deviation error estimate is propagated in this case.

Lastly the Validation module facilitates quality assessment of the IV products, by automating various standard validation tests, including internal consistency checks and intercomparisons with independent data sets (e.g. in-situ GPS, or other published velocity maps). The output is statistical information on the intercomparisons compiled in a quality assessment report.

1.3.3. Filtering and interpolation

For various reasons the OT algorithm sometimes fails to find matching features or erroneous matches are returned leading to gaps and/or errors in the velocity fields. This can for example be caused by a lack of traceable surface features or low coherence due to melt. For smoothing and removing outliers a 3x3 simple median filter is applied, after which a 9x9 first order plane fit filter is applied to fill in small data gaps. Further filtering/gap filling is left to the user if required.

1.3.4. Uncertainty characterisation

For incoherent offset-tracking (the general case), the uncertainty is characterised for each pixel by calculating a local offset-map standard deviation in a 5x5 neighbourhood. A plane fit to the offset map in the 5x5 neighbourhood is subtracted prior to calculating the standard deviation, so that an actual velocity gradient is not interpreted as a noise signal. The standard deviation estimate is corrected for any averaging carried out, as well as correlation between neighbouring samples (i.e. if the radar data are oversampled). The Greenland Ice Sheet velocity map is accompanied by its associated error standard deviation. The latter is also a map, in the same geometry as the associated measurement, providing a measure of uncertainty on a per-pixel basis.

1.4. Output data

The ice velocity product covers the Greenland Ice Sheet and is provided as a NetCDF file with 6 separate layers/maps representing different variables: the 3 velocity components: vx, vy, vz as well as vv (the magnitude of the horizontal components), and 2 gridded maps showing the valid pixel count and the uncertainty/standard deviation (Table 3). The ice velocity map is annually averaged and provided at 500m x 500m grid spacing in North Polar Stereographic projection (EPSG: 3413). A new NetCDF file of annually averaged IV is produced each year. The horizontal velocity is provided in true meters per day, towards easting (vx) and northing (vy) direction of the grid, and the vertical displacement (vz), is derived from a digital elevation model. More details are given in Section 1.3 of the C3S Product User Guide.

Table 3: IV main data variables

Variable name

Variable description

land_ice_surface_easting_velocity

Ice velocity East component [m/day]

land_ice_surface_northing_velocity

Ice velocity North component [m/day]

land_ice_surface_vertical_velocity

Ice velocity Vertical component [m/day]

land_ice_surface_velocity_magnitude

Ice velocity magnitude [m/day]

land_ice_surface_measurement_count

Valid pixel count [#]

land_ice_surface_velocity_stddev

Standard deviation [m/day]

2. Gravimetric mass balance – D3.IS.5

The gravimetric mass balance (GMB) is a brokered product from the ESA Climate Change Initiative (CCI). For details, the reader is referred to the to the GMB sections of the ESA Greenland_Ice_Sheet_cci and Antarctic_Ice_Sheet_cci products documentation, provided in Annex A.

2.1. Instruments

The primary satellite source for generation of the GMB is the Gravity Recovery and Climate Experiment (GRACE) (Tapley et al., 2004) satellite mission launched in 2002. The GRACE mission allows fluctuations in ice-sheet mass to be estimated through measurement of their changing gravitational attraction. Advantages of the GRACE method is that it provides regional averages without the need for interpolation, measures the effect of mass fluctuations directly, and permits monthly temporal sampling. However, a key challenge is to discriminate fluctuations in ice-sheet mass from changes in the underlying crust and mantle. The spatial resolution of GRACE observations derived from global spherical harmonic solutions of about 300 km in the Polar Regions is coarse in comparison to that of other geodetic techniques. Hence, a further complicating factor is that signals may leak into regional GRACE solutions as a consequence of remote geophysical processes.

In 2005, Velicogna and Wahr (2005) showed for the first time the possibility of using data from the GRACE mission to determine the mass balance of the Greenland ice sheet. Since then many mass balance estimates of both Greenland and Antarctica have been published, both on ice sheet scale (Chen et al., 2006; Ramillien et al., 2006; Forsberg and Reeh, 2007; Barletta et al., 2008; Velicogna, 2009; Shepherd et al. 2012) and drainage basin scale (Luthcke et al., 2006; Wouters et al., 2008; Schrama and Wouters, 2011; Sasgen et al., 2012; Barletta et al. 2013, Groh and Horwath 2016).

2.2. Input/auxiliary data and algorithms

Information on the data-use and specific algorithms can be found in Annex A, more specifically:

  • For the Antarctic Ice sheet: Annex A.1 Section 3 of the ESA Climate Change Initiative (CCI) Essential Climate Variable (ECV), Antarctic Ice Sheet (AIS) Algorithm Theoretical Basis Document (ATBD)
  • For the Greenland Ice sheet: Annex A.2 Section 6 of ESA Climate Change Initiative (CCI) Essential Climate Variable (ECV), Greenland Ice Sheet (GIS) Algorithm Theoretical Basis Document (ATBD)

2.3. Output data

The brokered data from the two ice sheet CCIs have been interpolated onto a common timeframe to provide a constantan product for both hemispheres. The GMB are given as a single NetCDF-file. Table 4 gives a brief description of the output fields, and a data sample.

Table 4: Summary of variables in the output fields, as given in the NetCDF-file for GMB, C3S_GMB_vers1_2018-12-18.nc

Variable

Short description

Unit

Long description

 Data sample

AntBasin

Basin definition

None

Zwally definition of major drainage basins of the Antarctic Ice Sheet

AntIS_x

Basin x GMB

Gt

Antarctic gravimetric mass balance for basin x

AntIS_x_er

Basin x GMB error

Gt

Error estimate of the Antarctic gravimetric mass balance for basin x

AntIS_total

Total GMB

Gt

Total Antarctic gravimetric mass balance

AntIS_total_er

Total GMB error

Gt

Error estimate of total Antarctic gravimetric mass balance

AntPeninsula_total

Total GMB

Gt

Total gravimetric mass balance for the Peninsula

AntPeninsula_total_er

Total GMB error

Gt

Error estimate of total gravimetric mass balance for the Peninsula

East_AntIS_total

Total GMB

Gt

Total gravimetric mass balance for the East Antarctica

East_AntIS_total_er

Total GMB error

Gt

Error estimate of total gravimetric mass balance for the East Antarctica

West__AntIS_total_er

Total GMB

Gt

Total gravimetric mass balance for the west Antarctica

West__AntIS_total_er

Total GMB error

Gt

Error estimate of total gravimetric mass balance for the West Antarctica

GrISBasin

Basin definition

None

Zwally definition of major drainage basins of the Greenland Ice Sheet

GrIS_x

Basin x GMB

GT

Greenland gravimetric mass balance for basin x

GrIS_x_er

Basin x GMB error

GT

Error estimate of the Greenland gravimetric mass balance for basin x

GrIS_total

Total GMB

GT

Total gravimetric mass balance for the Greenland ice sheet

GrIS_total_er

Total GMB error

GT

Error estimate of total gravimetric mass balance for the Greenland Ice sheet

time

Time

Hours since 01-01-1990


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

3.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 satellite and will extend into the future with CryoSat-2 and the Sentinel-3 series of satellites. Due to the ongoing baseline D reprocessing of Cryosat-2, the surface elevation change processing still utilises the current Baseline C. Baseline C provides data until April 2019, and after this date the surface elevation change record relies solely on elevation estimates from the Sentinel-3A satellite. 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. Until then, we flag the lower data coverage by giving the distance between the satellite observation and the gridded surface elevation change solution.

ESA recently switched its processing baseline from C to D. These two baselines differ in several respects, notably the slope mask used in low resolution mode. Because of the differences, baselines C and D cannot be combined. Baseline D has only been issued since April 2019, so its time-range is not yet long enough to cross-calibrate accurately with other missions. Consequently, only baseline C data is included in the SEC v2 CDR. When the baseline D full mission reprocessing is available, it can be incorporated into the next CDR. Sentinel-3A data is included in the v2 CDR, but a reprocessed release from ESA is expected, to resolve problems in transitioning from ocean to land ice surfaces

The Antarctic region surface elevation change dataset is derived from the crossover method (see Section 3.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 5 lists the missions/mission phases used as input to both the Antarctic and Greenland SEC datasets.

Table 5: Mission summary

Mission

Instrument

Period used

Southern observation limit

Repeat cycle

ERS1 phase C

RA

April 1992 to December 1993

81.5°

35 days

ERS1 phase G

RA

April 1995 to May 1996

81.5⁰

35 days

ERS2

RA

July 1995 to June 2003

81.5°

35 days

EnviSat

RA-2

October 2002 to October 2010

81.4°

35 days

CryoSat-2

SIRAL-2

November 2010 to April 2019

88.0°

369 days, with 30-day sub-cycle

Sentinel 3A

SRAL

December 2016 to present

81.4°

27 days

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

Table 6: ERS1 and ERS2 RA

Parameter

Value

Frequency

13.8 GHz (Ku band)

Pulse repetition frequency

1.02 kHz

Pulsewidth

20 µs chirp

Bandwidth in ice mode

82.5MHz

Range resolution

10cm

Beam width

1.3⁰

Footprint (pulse-limited)

16 to 20 km

Table 7: EnviSat RA-2

Parameter

Value

Frequency

13.575 GHz (Ku band)

Pulse repetition frequency

1.796 kHz

Pulsewidth

20 µs chirp

Bandwidth

320MHz

Range resolution

50cm

Beam width

1.3⁰

Footprint (pulse-limited)

2 to 10 km

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

Parameter

Value

Frequency

13.575 GHz (Ku band)

Pulse repetition frequency

1.97 kHz

Pulsewidth

50 µs chirp

Bandwidth

320MHz

Range resolution

45cm

Beam width

1.2⁰ across-track, 1.08⁰ along-track

Footprint (pulse-limited)

2km across-track, 2km along-track

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

Parameter

Value

Frequency

13.575 GHz (Ku band)

Pulse repetition frequency

17.8 kHz

Pulsewidth

50 µs chirp

Bandwidth

320MHz

Range resolution

45cm

Beam width

1.2⁰ across-track, 1.08⁰along-track

Footprint (pulse-limited)

2km across-track, 300m along-track

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

Parameter

Value

Frequency

13.575 GHz (Ku band)

Pulse repetition frequency

17.8 kHz

Pulsewidth

50 µs chirp

Bandwidth

350MHz

Range resolution

3cm

Beam width

~1.3⁰

Footprint (pulse-limited)

1.64km across-track, 300m along-track

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

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

3.2. Input and auxiliary data

3.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_v2.1, see Femenias (ed), 2018
  • CryoSat-2 SIRAL L2i Low Resolution Mode (LRM) and SAR-Interferometric (SIN) datasets, see Bouzinac, C, 2018
  • Sentinel-3A L2 SRAL, specifically the SR_2_LAN_NT products containing land ice data, see ACRI-ST IPF Team, 2018

It is planned to reprocess earlier data as new datasets become available – specifically Envisat's GDR_v3, CryoSat-2's baseline D and Sentinel-3A's land ice product. See the related document, the Development Milestone Plan for the intended reprocessing schedule.
Data is freely available on registration from these websites:

3.2.2. Auxiliary data

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

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

3.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. These are described in Zwally, Giovinetto et al, 2012, and freely available from:

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

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

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

3.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;

3.3. Algorithms

3.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 3.3.1 above, 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 3.3.1 of the related document the Product User Guide 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.

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.

3.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 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 2 shows the physical layout of a single crossover.


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

$$dH(t)=H_2-H_1+E $$

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.

$$\Delta h(x,t,t_{ref})= \frac{1}{2} \left[ \left(h_{At}-h_{Dt_{ref}} \right) + \left(h_{At_{ref}} - h_{Dt} \right) \right]_{t=t_1…t_N} $$

where h is height, Atref and At are the ascending reference and comparison tracks respectively, and Dtref and Dt 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-3A 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.

$$\Delta p(x,t,t_{ref})= \frac{1}{2} \left[ \left(p_{At}-p_{Dt_{ref}} \right) + \left(p_{At_{ref}} - p_{Dt} \right) \right]_{t=t_1…t_N} $$

where p is backscatter power, Atref and At are the ascending reference and comparison tracks respectively, and Dtref and Dt 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-3A datastream, a 2 year period has to be used.

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

$$dh=a_0+a_1dp $$

where dh is height difference to the reference cycle, dp is backscatter power difference to the reference cycle and a0 and a1 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

$$dh_{corr}=dh-(dp \ast \frac{dh}{dp})$$

where dhcorr 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

$$dh=a_0+(a_1×t)+a_2 \sin(a3 ×2\pi × t+a_4)$$

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.

3.3.3. Multi-mission cross-calibration

As expected in the Development Milestone Plan, the cross-calibration method was updated for CDR v2. A new algorithm, using multiple regression (for general details see Tabachnick and Fidell, 2019), was compared to an evolution of the v1 algorithm. The new method was judged better, as it enabled recovery of approximately 16% more SEC values from the v1 inputs than the original linear-fit method, with similar accuracy where both methods produced results. See Figure 3 below, which shows histograms of total uncertainty for both methods, and a plot of the percentage coverage of the Antarctic Ice Sheet, aggregated into one-year periods to accommodate CryoSat-2's drifting orbit. Note that an increase in data recovery is not equivalent to the same increase in area coverage, as gaps in the timeseries for single cells were also filled.


Figure 3: Histogram of total uncertainty on SEC values (left) and plot of percentage coverage of the AIS, aggregated per year (right), in datasets derived from the same inputs by the v1 and v2 methods. V1, linear fit, is shown in black, v2, multiple regression, in red.

See also the related Product Quality Assurance Document, Section 4.3, for comparison of the v1 and v2 methods with respect to the Operation IceBridge validation dataset.

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 4 below. This shows the 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 4: 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.

One of the major advantages of this method over the v1 algorithm is that all the datapoints in the timeseries can be used. Another is that even if one mission is unrepresented in the timeseries, it is still possible to cross-calibrate the rest. Mission datapoints are not inferred from the model. Biases are added to pre-existing datapoints only.

3.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 3.3.3 above). This is done to produce a basin level accuracy indicator (see Section 3.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 (although they may be in future revisions) but can be retrieved from the provided dh/dt timeseries.

3.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 3.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 3.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 (3.3.3 and 3.3.4) above.

3.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 3.3.1 of the related document, the Product User Guide.

4. 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 see Annex A.2. 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. 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).

4.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 3.1 describes the instruments used from each mission.
Table 5 in Section 3.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 6 to Table 10. The only difference is the choice of operational modes for CryoSat-2 which is limited to LRM, and SARIn over Greenland, in contrast to the Antarctic ice sheet.

4.2. Input and auxiliary data

4.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 3.2.1 for details.

4.2.2. Auxiliary data

4.2.2.1. Digital elevation model for Greenland

The Greenland surface elevation change applies the official level-2 data solutions provided by ESA, as seen above. The generation of such level-2 product, includes the use of a digital elevation model (DEM) in the geolocation of LRM data. We here refer to the mission documentation for the specific DEM used in the geolocation of the radar echo. Following Sørensen et al. (2018), no DEM is used for the combined cross-over and repeat-track solutions for ERS-1, ERS-2 and ENVISat. Both the along-track and plan-fitting algorithms need an a-priori DEM to be subtracted from the raw geolocated radar altimetric observation. Here, we use the Greenland Ice sheet Mapping Project (GIMP) DEM version 1 (Howat, Negrete, and Smith 2014).

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

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

4.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 3.2.3 for details.

4.3. Algorithms

The theoretical basis for the Greenland ice sheet SEC are following 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).

4.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 accuracy target is 0.1 m/year and the stability target is also 0.1 m/year. All product generation are in accordance to the equidistant grid in polar stereographic projection defined by EPSG: 3413 (see Product User Guide Section 4.3.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 Annex A.2. 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 to 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 algorithm specific equations for the four algorithms are illustrated in Figure 5 and summarised in the following subsections.


Figure 5: 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 is along the same exact track, and along-track when the tracks differ in position, then the "green" observations are projected onto the blue reference track. c) Plan-fitting, here the solution is modelled within regular-shaped grid cells. The figure is modified from Moholdt et al., (2010)

4.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 parameter, which have been shown to correlate to surface elevation change, e.g. backscatter (Bs) and seasonality in the altimetry data. Here, we according to Sørensen et al (2018) 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

$$\begin{split} H(x,y,t) & = H_0(\overline{x}, \overline{y}) + dH/dt(t - \overline{t}) \\ & + dB_s(Bs - \overline{Bs}) \\ & + sx(x - \overline{x}) + sy(y - \overline{y}) \\ & + \alpha \cos(\omega t) + \beta \sin( \omega t) \\ & + \epsilon(x,y,t), \end{split} $$ where x, y, t is the spatial and temporal components, H is the surface elevation as measured by the satellite, $H_0$ is the mean elevation of the evaluated grid cell, sx, sy are curvature terms along the along-track segment, the co-sine term including $\alpha$, $\beta$ and $\omega$ descripts the seasonality in the surface elevation changes. $\epsilon$ is the residual, which is minimised in the derivation of the surface elevation change along-track. We implement the along-track segments as overlapping segments to increase resolution.
4.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.

4.3.1.3. Cross-over

See Section 3.3.2 for a detailed description of the cross-over (XO) algorithm. Here, the time-series are produced for 0.2ox0.5o latitude and longitude grids (Khvorostovsky, 2012).

4.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,

$$\begin{split} H(x,y,t) & = H_0(\overline{x}, \overline{y}) + dH/dt(t - \overline{t}) \\ & + dLeW(LeW - \overline{LeW}) \\ & + sx(x - \overline{x}) + sy(y - \overline{y}) + cx(x^2- \overline{x}^2) \\ & + cy(y^2 - \overline{y}^2) + cc(x^2 - \overline{x}^2)(y^2- \overline{y}^2) \\ & + \alpha \cos(\omega t) + \beta \sin( \omega t) \\ & + b_{AD}(-1)^{AD} + b_m(-1)^m\\ & + \epsilon(x,y,t), \end{split}$$

where also the additional biases for CryoSat-2 (bAD ascending/descending and bm LRM/SARIN mode) were added alongside with replacing the Backscatter bias with biases in changing leading edge width of the recorded waveform (LeW).

4.3.1.5. Adapting the plane-fitting algorithm for Sentinel-3A

For the processing of Sentinel-3 the processing plan-fitting algorithm applied for Cryosat-2 was chosen. 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 plan-fitting algorithm:

$$\begin{split} H(x,y,t) & = H_0(\overline{x}, \overline{y}) + dH/dt(t - \overline{t}) \\ & + dLeW(LeW - \overline{LeW}) + dB_s(Bs - \overline{Bs}) \\ & + sx(x - \overline{x}) + sy(y - \overline{y}) + cx(x^2- \overline{x}^2) \\ & + cy(y^2 - \overline{y}^2) + cc(x^2 - \overline{x}^2)(y^2- \overline{y}^2) \\ & + \alpha \cos(\omega t) + \beta \sin( \omega t) \\ & + b_{AD}(-1)^{AD} \\ & + \epsilon(x,y,t), \end{split}$$

 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 plan-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 plan-fitting algorithm.

Merging of the created surface elevation grids and processing

The round robin preformed as a part of the ESA CCI Greenland ice sheet project (see Appendix A) showed collocation as the optimal interpolation method for combining heterogeneous data of different kinds.

For the second version of the surface elevation change processor the collocation-gridding procedure was replaced by an ordinary kriging method to preform 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.

Figure 6: Schematic drawing of the combined product generation. 

At the end of 2019, and with the inclusion of Sentinel-3A into the data procession, the internal production of surface-elevation-change grids consists of 98 Greenland maps. These grids need to be combined to create the seamless time-series of surface elevation change from 1992-present at monthly time steps. As illustrated in Figure 6 the processing is split in 3-phasis:

1. For the older missions, ERS-1, ERS.2, and ENVISat:

The internal product is a 5-year running mean, so surface elevation change grids are shifted 1-year apart, giving multiple maps of surface elevation change which cover a specific month in question. We apply a weighted averaging of the maps to produce the given map for the month in question. E.g. when the solution for December 1996 is derived, the maps of 1993-1997, 1994-1997, 1995-1998 and 1996-1999 are used. This ensures a stable progression in month-to-month solution of the surface elevation change for the older and nosier satellite data.

2. For CryoSat-2 alone:

The novel orbit of CryoSat-2 ensures that the altimeter provides enough coverage in Greenland for the PF-algorithm to be stable when 3-years of elevation observations are evaluated. Hence, the strategy for CryoSat-2 differs from the older satellites and monthly time-increments surface elevation change grids are derived by fitting 1.5-year of data. This gives a truly monthly solution only possible due to the novel altimeter design (Simonsen and Sørensen 2017).

3. For CryoSat-2 and Sentinel-3:

In similar fashion as for CryoSat-2, the processing for Sentinel-3A are done for monthly time-increments and the surface elevation change grids are derived by fitting ±1.5-year of data. The CryoSat-2 and Sentinel-3 grids are then merged by weighted-averaging the grids in accordance to the error estimates.

Based on the monthly surface elevation grids an accumulated surface elevation change is also provided in the C3S-product.

4.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 bias the derived solution toward more positive values. 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. Especially the errors induced from the time varying penetration of the surface snow and the geo-location of the radar echo is 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).

4.4. Output data

The output CDR is contained in a single netCDF4 file, which after the initial CDR will be 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 cantered on the time given for that grid. This 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, Section 4.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

Bamler, R., M. Eineder Accuracy of differential shift estimation by correlation and split-bandwidth
interferometry for wideband and delta-k SAR systems, IEEE Geosc. Rem Sens. Lett. 2, 151-155 (2005).

Barletta, V. R., Bordoni, A., and Sabadini, R.: Isolating the PGR signal in the GRACE data: impact on mass balance estimates in Antarctica and Greenland, Geophys. J. Int., 172, 18–30, doi:10.1111/j.1365- 246X.2007.03630.x, 2008.

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.

Bouzinac, C, 2018, CryoSat Product Handbook, https://earth.esa.int/documents/10174/125272/CryoSat_Product_Handbook

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

Chen, J. L., Wilson, C. R., Blankenship, D. D., and Tapley, B. D.: Antarctic mass rates from GRACE, Geophys. Res. Lett., 33, L11502, doi:10.1029/2006GL026369, 2006.

CryoSat product handbook, 2018. This is still 'under definition' in the ESA Earth Online document library, but see Bouzinac, 2018, above for an alternative source

Davis, C., and Ferguson, A., 2994. Elevation change of the Antarctic Ice Sheet 1995-2000 from ERS2 satellite radar altimetry. IEEE Trans. Geosci. Remote Sens. 42 (11), 2437-2445. DOI: 10.1109/TGRS.2004.836789

EnviSat product handbook, 2011, in the ESA Earth Online document library at https://eart.esa.int/web/guest/document-library

Femenias (editor), 2018, Envisat Altimetry Level 2 Product Handbook v2, ESA document CLS - ESLF - 18 -0003

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

Forsberg, R. and Reeh, N.: Mass change of the Greenland Ice Sheet from GRAC E, in: Proceedings, Gravity Field of the Earth – 1st meeting of the International Gravity Field Service, Harita Dergisi, Ankara, vol. 73, avaialble at: http://www.igfs.net, in review, 2007

Groh, A., and Horwath, M. (2016). The method of tailored sensitivity kernels for GRACE mass change estimates. Geophysical Research Abstracts, 18, EGU2016-12065.

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

Luthcke, S. B., Zwally, H. J., Abdalati, W., Rowlands, D. D., Ray, R. D., Nerem, R. S., Lemoine, F. G., McCarthy, J. J., and Chinn, D. S.: Recent Greenland ice mass loss by drainage system from satellite gravity observations, Science, 314, 1286–1289, 2006.

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, Rott, H., Hetzenecker, M., Wuite, J. and Potin, P., 2015, The Sentinel-1 Mission: New Opportunities for Ice Sheet Observations. Remote Sens., 7, 9371-9389.

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

Ramillien, G., Lombart, A., Cazenave, A., Ivins, E. R., Llubes, M., Remy, F., and Biancale, R.: Interannual variations of the mass balance of the Antarctica and Greenland ice sheets from GRACE, Global Planet. Change, 53, 198–208, doi:10.1016/j.gloplacha.2006.06.003, 2006.

Rastner, P., Bolch, T., Mölg, N., Machguth, H., Le Bris, R., Paul, F. (2012): The first complete inventory of the local glaciers and ice caps on Greenland. The Cryosphere, 6, 1483-1495. (doi:10.5194/tc-6-1483-2012)

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

Rizzoli, P., Martone, M., Gonzalez, C., Wecklich, C., Borla Tridon, D., Bräutigam, B., Bachmann, M., Schulze, D., Fritz, T., Huber, M., Wessel, B., Krieger, G., Zink, M., and Moreira, A. (2017): Generation and performance assessment of the global TanDEM-X digital elevation model.
ISPRS Journal of Photogrammetry and Remote Sensing, Vol 132, pp. 119-139.

Sasgen, I., Broeke, M. v. d., Bamber, J. L., Rignot, E., Sandberg Sørensen, L., Wouters, B., Martinec, Z., Velicogna, I., and Simonsen, S. B.: Timing and origin of recent regional ice-mass loss in Greenland, Earth Planet. Sci. Lett., 333–334, 293–303, doi:10.1016/j.epsl.2012.03.033, 2012.
Schrama, E. J. O. and Wouters, B.: Revisiting Greenland ice sheet mass loss observed by GRACE, J. Geophys. Res., 116, B02407, doi:10.1029/2009JB006847, 2011.

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.
Tapley, B. D., Bettadpur, S., Watkins, M., and Reigber, C. (2004). The gravity recovery and climate experiment: Mission overview and early results. Geophys. Res. Lett., 31, L09607.

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
Velicogna, I. and Wahr, J.: Greenland mass balance from GRACE, Geophys. Res. Lett, 32, L18505, doi:10.1029/2005GL023955, 2005.
Velicogna, I.: Increasing rates of ice mass loss from the Greenland and Antarctic ice sheets revealed by GRACE, Geophys. Res. Lett, 36, L19503, doi:10.1029/2009GL040222, 2009.

Wessel, B. (2018) TanDEM-X Ground Segment – DEM Products Specification Document", EOC, DLR, Oberpfaffenhofen, Germany, Public Document TD-GS-PS-0021, Issue 3.2. Available: https://geoservice.dlr.de/web/dataguide/tdm90/pdfs/TD-GS-PS-0021_DEM-Product-Specification.pdf

Wessel, B., Huber, M., Wohlfart, C., Marschalk, U., Kosmann, D., & Roth, A. (2018). Accuracy assessment of the global TanDEM-X Digital Elevation Model with GPS data. ISPRS Journal of Photogrammetry and Remote Sensing, 139, 171–182. doi:10.1016/j.isprsjprs.2018.02.01

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

Wouters, B., Chambers, D., and Schrama, E. J. O.: GRACE observes small-scale mass loss in Greenland, Geophys. Res. Lett., 35, L20501, doi:10.1029/2008GL034816, 2008.

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

Appendices– ESA CCI Ice Sheets ATBD

The GMB product is brokered from the ESA Climate Change Initiative two ice sheets projects. In the following the relevant ATBDs from the two projects are brokered.

Annex A.1 ESA Climate Change Initiative (CCI) Essential Climate Variable (ECV), Antarctic Ice Sheet (AIS) Algorithm Theoretical Basis Document (ATBD)


file not available from pdf in CDS documentation tab 

Annex A.2 ESA Climate Change Initiative (CCI) Essential Climate Variable (ECV), Greenland Ice Sheet (GIS) Algorithm Theoretical Basis Document (ATBD)


file not available from pdf in CDS documentation tab 

This document has been produced in the context of the Copernicus Climate Change Service (C3S).

The activities leading to these results have been contracted by the European Centre for Medium-Range Weather Forecasts, operator of C3S on behalf of the European Union (Delegation Agreement signed on 11/11/2014 and Contribution Agreement signed on 22/07/2021). All information in this document is provided "as is" and no guarantee or warranty is given that the information is fit for any particular purpose.

The users thereof use the information at their sole risk and liability. For the avoidance of all doubt , the European Commission and the European Centre for Medium - Range Weather Forecasts have no liability in respect of this document, which is merely representing the author's view.

Related articles