Estimating Arctic sea ice thickness and volume using CryoSat-2 radar altimeter data

Arctic sea ice is a major element of the Earth’s climate system. It acts to regulate regional heat and freshwater budgets and subsequent atmospheric and oceanic circulation across the Arctic and at lower latitudes. Satellites have observed a decline in Arctic sea ice extent for all months since 1979. However, to fully understand how changes in the Arctic sea ice cover impact on our global weather and climate, long-term and accurate observations of its thickness distribution are also required. Such observations were made possible with the launch of the European Space Agency’s (ESA’s) CryoSat-2 satellite in April 2010, which provides unparalleled coverage of the Arctic Ocean up to 88 (cid:1) N. Here we provide an end-to-end, comprehensive description of the data processing steps employed to estimate Northern Hemisphere sea ice thickness and subsequent volume using CryoSat-2 radar altimeter data and complementary observations. This is a sea ice processor that has been under constant development at the Centre for Polar Observation and Modelling (CPOM) since the early 1990s. We show that there is no signiﬁcant bias in our satellite sea ice thickness retrievals when compared with independent measurements. We also provide a detailed analysis of the uncertainties associated with our sea ice thickness and volume estimates by considering the independent sources of error in the retrieval. Each month, the main contributors to the uncertainty are snow depth and snow density, which suggests that a crucial next step in Arctic sea ice research is to develop improved estimates of snow loading. In this paper we apply our theory and methods solely to CryoSat-2 data in the Northern Hemisphere. However, they may act as a guide to developing a sea ice processing system for satellite radar altimeter data over the Southern Hemisphere, and from other Polar orbiting missions. 2017 COSPAR. Published by Elsevier


Introduction
Satellite passive microwave observations of Arctic sea ice have recorded a decline in the summer extent of $40% since 1979 (Cavalieri et al., 1996, updated yearly;Fetterer et al., 2002, updated daily). The decline is coincident with abrupt global and Arctic warming over the last 30 years (Hartmann et al., 2013). It is crucial to observe and understand changes in the Arctic sea ice cover, as it is a major element of the Earth's climate system. Sea ice influences the freshwater (Aagaard and Carmack, 1989;Serreze et al., 2006) and surface heat (Sedlar et al., 2011) budgets of the Arctic, and subsequently the global climate. The melting of sea ice could disrupt the oceanic global thermohaline circulation (Vellinga and Wood, 2002) and atmospheric circulation patterns (Singarayer et al., 2006;Schweiger et al., 2008;Francis and Vavrus, 2012), with knock-on effects for regional weather patterns in Europe, America and much of the northern hemisphere, and potentially the southern hemisphere (Vellinga and Wood, 2002). Arctic sea ice cover, long-term and accurate observations of the ice pack as a whole are required. However, it has previously been difficult to quantify trends in sea ice volume because detailed thickness observations have been spatially sparse and temporally sporadic (McLaren, 1989(McLaren, , 1992Wadhams, 1990).
In April 2010 the European Space Agency's (ESA) CryoSat-2 radar altimeter satellite (Wingham et al., 2006) was launched. CryoSat-2 provides unparalleled coverage of the Arctic Ocean with a region of coverage (ROC) extending to 88°N. In 2013, a study led by scientists at the Centre for Polar Observation and Modelling (CPOM) produced the first estimates of Arctic sea ice thickness and volume from CryoSat-2 (Laxon et al., 2013). The estimates were produced within a fixed central Arctic region that covers an area of $7.2 Â 10 6 km 2 . The region was first defined by scientists at the National Aeronautics and Space Administration (NASA) for use with the NASA ICESat satellite (Kwok et al., 2009), and will herein be referred to as the ICESat domain. The CPOM processor has since been updated to cover all northern hemisphere sea ice, defined as sea ice at latitudes above and including 40°N . Since 2013 a number of institutions have published results from their own CryoSat-2 sea ice processing systems, including the National Aeronautics and Space Administration (NASA) Goddard Space Flight Centre , the Alfred Wegener Institute (AWI) (Ricker et al., 2014), NASA Jet Propulsion Laboratory (JPL) (Kwok and Cunningham, 2015), the Finnish Meteorological Institute (FMI) (Rinne and Similä, 2016), and the Ulsan National Institute of Science and Technology (UNIST) (Lee et al., 2016). These all differ slightly in their areal coverage and processing technique. Like CPOM, both NASA teams and AWI provide a publicly available sea ice thickness product. The NASA data are available within the ICESat domain whereas the AWI data cover the area of the Arctic Ocean where values of snow depth and density from a climatology (Warren et al., 1999) are considered realistic. Another key difference is the retrackers applied to CryoSat-2 data for each processor. At CPOM we estimate the elevation of ocean and ice surfaces by applying a Gaussian-exponential retracker to ocean waveforms and a threshold retracker to sea ice waveforms (Section 4.4). NASA Goddard have developed a waveform fitting model, NASA JPL select the first unambiguous peak for all waveforms, and AWI apply a threshold retracker to all waveforms. Other variations include the method used to account for the mean sea surface when estimating instantaneous sea surface height, and to classify sea ice type.
Here we provide an end to end, comprehensive description of the data processing steps that we currently employ at CPOM to estimate northern hemisphere sea ice thickness and volume using CryoSat2 radar altimeter data. This is a sea ice processing chain that has been under constant development at CPOM since the early 1990s (Laxon, 1994). Past studies have documented aspects of its evolution (Giles et al., 2008;Laxon et al., 2013Laxon et al., , 2003Peacock and Laxon, 2004;Tilling et al., 2015) and provided a detailed analysis of sources of error and uncertainty in the retrieval of sea ice freeboard from satellite radar altimetry (Tilling et al., 2016;Giles et al., 2007). In this paper we also develop an uncertainty budget for northern hemisphere sea ice thickness and volume, evaluate of our thickness product by comparison with in situ and airborne Arctic sea ice measurements, and present an assessment of the changes in sea ice thickness and volume from CryoSat-2.

The CryoSat-2 satellite
The primary aim of ESA's CryoSat-2 mission is to accurately determine the inter-annual fluctuations and longerterm trends in Earth's continental and marine ice fields (Drinkwater et al., 2004;Wingham et al., 2006). Its primary payload is a Synthetic Aperture Interferometric Radar Altimeter (SIRAL). SIRAL is an altimeter/interferometer system operating in the . The SIRAL antenna system comprises two nadir looking antenna mounted 1 m apart in the across-track direction (Wingham et al., 2006). It operates in three modes ( Fig. 1) -low resolution mode (LRM), synthetic aperture radar (SAR) mode, and SAR interferometric (SARIn) mode -depending on the type of surface that is being observed (ESA/MSSL, 2013).
In LRM a single antenna is used to transmit and receive the radar signal and SIRAL acts as a conventional nadirpointing, pulse-limited altimeter. This means that the radar footprint size is dependent on the length of the compressed pulse. The typical CryoSat-2 orbital velocity is 7.4 km s À1 and the interval between pulses in LRM is approximately 500 ls, which corresponds to a pulse repetition frequency (PRF) of 2 kHz (Wingham et al., 2006). This ensures that returning echoes are uncorrelated. The LRM footprint is the SIRAL pulse-limited footprint (PLF), which is approximately 1.7 km based on an altitude of 730 km (Scagliola, 2013). CryoSat-2 operates in LRM over areas of the continental ice sheets, and the majority of the Earth's ice-free oceans and land (Fig. 1).
When operating in SAR mode, the SIRAL instrument on-board CryoSat-2 uses a single antenna to transmit and receive pulses, but emits a burst of 64 phase-coherent pulses as opposed to a single pulse. By exploiting the slight frequency shifts caused by the Doppler effect, in the forward-and aft-looking parts of the burst, the data processor can separate the burst into narrow beams arranged across-track. The along-track sampling resolution of each beam is approximately 250 m. As with other pulse-limited radar systems, the surface area illuminated by SIRAL continues to grow as an annulus following the time that the PLF is reached. Therefore the PLF is smaller than the full antenna illumination pattern, or antenna-limited footprint. The across-track footprint in SAR mode is simply the antenna-limited footprint, which can reach 15 km depending on satellite altitude (Scagliola, 2013). The beams from each burst are made to overlap with the beams from other bursts in exact coincidence by adjusting the look angle of the central beam (the angle at which the beam ''looks" at the surface), and in turn all other beams. Over multiple bursts a given beam location is sensed multiple times. All waveform returns from beams at the same location are gathered together in a 'stack', and their delay time and power have a dependence on look angle. Therefore, waveform delay times are adjusted to account for the different ranges to the surface in a process known as slant range correction, and powers are weighted depending on look angle. The power of each return waveform in the stack also has a dependence on incidence angle (the angle between the incident beam and the normal to the intercepted surface), which results partly from variations in surface backscattering. The standard deviation of backscattering variation for all waveforms that make-up the stack is known as the stack standard deviation (SSD), which is an important parameter for the distinction of different surface types. The waveforms from each stack are averaged into one composite waveform -a 20 Hz waveform -in a process know as 'multi-looking' (ESA/MSSL, 2013). The 20 Hz waveforms are provided with the CryoSat-2 Baseline-C Level 1b (L1b) data product (ESA/MSSL, 2013;ESA/ACS, 2011), along with SIRAL instrument information relating to each waveform. The burst-mode SAR processing adopted by SIRAL is described in detail in Raney (1998) and Wingham et al. (2006). SAR mode is implemented over sea ice areas as they are relatively flat, and over some ocean basins and coastal zones (Fig. 1). The application of SAR technology over sea ice provides a larger number of independent measurements with finer resolution compared to a conventional pulselimited system, to improve the retrieval accuracy.
In SARIn mode, the along-track processing remains the same as SAR mode. However, the second antenna is activated to receive additional across-track information. Following pulse emission, both antennae receive the return signal almost simultaneously. If the signal originates from anywhere on the ice surface other than satellite nadir, then there will be a difference in the path lengths of the signal received at each antenna. It is therefore possible to derive the angle between the baseline (joining the antennas) and the echo direction (origin) by interferometry. SARIn mode is used to estimate the across-track surface slope to account for slope-induced errors when estimating range. SARIn mode is mainly implemented over ice sheet margins, small ice caps and mountain glaciers, which often have large slope variations. It is also used over some geostrophic ocean currents, major hydrological river basins, and coastal zones. Prior to October 2014, SARIn mode was also used over a small area of sea ice north of Ellesmere Island, Nunavut (ESA/MSSL, 2013) ( Fig. 1) and has been used to investigate the influence of off-nadir ranging to leads in sea ice elevation estimates (Armitage and Davidson, 2014).

CryoSat-2 data
To estimate Arctic sea ice thickness and volume, we use CryoSat-2 Baseline-C L1b SAR and SARIn mode data (ESA/ACS, 2011;Scagliola and Fornari, 2015), available from ESA via an ftp client (ftp://science-pds.cryosat.esa. int). The L1b data contain measurement information for each 20 Hz waveform along the orbit ground track of the Fig. 1. CryoSat-2 operation modes across the Arctic. CryoSat-2 operates in low resolution mode (LRM; blue tracks) over the continental ice sheets and the majority of the Earth's ice-free oceans and land, synthetic aperture radar (SAR) mode (yellow tracks) over sea ice and some ocean basins and coastal zones, and SAR interferometric (SARIn) mode (pink tracks) over ice sheet margins, small ice caps and mountain glaciers, and some geostrophic ocean currents, major hydrological river basins and coastal zones. (a) The mode mask for orbit tracks in April 2011, when the SARIn mode was in operation over a small area of sea ice north of Ellesmere Island. (b) The mode mask for orbit tracks in April 2015, once the SARIn mode was no longer in operation over sea ice. This is the mode mask that has been implemented since October 2014.
satellite. Each ground track constitutes around 3000 20 Hz waveforms, on average. For SAR and SARIn modes, the L1b data include the average multi-looked echo power, satellite altitude, window delay, measurement time, geolocation, geophysical corrections, a measurement confidence data (MCD) flag (to indicate a number of potential problems that may arise with each waveform), surface type model flag (open ocean, enclosed sea, continental ice, or land), SSD (Section 2), and numerous other SIRAL instrument measurements (ESA/MSSL, 2013). In the L1b data ESA have applied the precise satellite orbit and instrument corrections to the window delay and satellite altitude computations, respectively, and they are also provided in the product. The geophysical corrections provided in the product have not been applied to the window delay or satellite altitude. In SARIn mode the multi-looked echo is complex, and therefore the L1b product also contains phase differences and coherence terms for each bin of the echo, computed by comparing the echoes received by both antennas. In SAR mode the CryoSat-2 range window explored is about 60 m, which corresponds to 256 range bins in the waveform data. In SARIn mode the range window is increased to about 240 m in order to capture the slope variation in ice sheet margins, which corresponds to 1024 range bins. ESA also produce a Level 2 (L2) data product, which provides waveform information derived from the L1b data. L2 data includes estimates of elevation for all surface types, as well as other surface parameters such as the radar backscattering coefficient.

Sea ice concentration
We use daily sea ice concentration data that are generated at the NASA Goddard Space Flight Centre (GSFC) and are available through the National Snow and Ice Data Centre (NSIDC) (Cavalieri et al., 1996, updated yearly). The data are generated from brightness temperature data derived from satellite passive microwave sensors using the NASA Team algorithm. Gridded averages of the percentage of ocean area covered by sea ice are provided on a Polar stereographic projection with a grid size of 25 km square.

Sea ice type
Sea ice type data are provided by the Norwegian Meteorological Service (NMS) Ocean and Sea Ice Satellite Application Facility (OSI SAF) (Andersen et al., 2012). The NMS OSI SAF data defines sea ice type as first year ice (FYI), multi-year ice (MYI), ambiguous, ice free, or unclassified. The sea ice type is determined by combining satellite measurement of brightness temperature and backscatter. For brightness temperature, the gradient ratio (the normalised difference in brightness temperature between the 37 and 19 GHz vertically polarised channels) can be used to distinguish between ice types. At vertical polarization the brightness temperature of FYI is similar to that of MYI at radiation frequencies of $5 GHz. The brightness temperature of MYI then drops as frequency increases, due to increased internal scattering (Eppler et al., 1992;Svendsen et al., 1983). MYI tends to be rougher than FYI, meaning that backscatter is larger over MYI. MYI also has an additional backscatter signature compared to FYI as a result of volume scattering, particularly during winter (Onstott, 1992).

Sever expedition sea ice data
Sever expedition data (NSIDC, 2004) provide in situ measurements of sea ice freeboard, sea ice thickness and snow depth. These data consist of mean values of each parameter on sea ice runways from 689 aircraft landings between 1982 and 1988. The measurements were made in spring over level, predominantly FYI in the Eurasian Russian Arctic (Alexandrov et al., 2010).

Snow loading climatology
We obtain snow loading data from a climatology (Warren et al., 1999). The climatology was compiled from in situ measurements of snow depth and snow density collected over MYI in the central Arctic between 1954 and 1991, with a two-dimensional quadratic function fitted to all measurements to represent the spatial variability of depth ( Fig. 2a-d) and density ( Fig. 2e-h). However, these quadratic functions are not constrained by in situ measurements outside of the central Arctic ( Fig. 2a and e) and there are known differences between the climatology and the current snow depth on younger Arctic sea ice (Kurtz and Farrell, 2011;Webster et al., 2014). Therefore, for each month we apply the mean climatology values of snow depth and density from within the central Arctic region to all sea ice measurements. We then half snow depth over FYI (Kurtz and Farrell, 2011), to account for the reduced snow accumulation compared with MYI.

Sea ice thickness and volume processing method
The Arctic sea ice thickness and volume processor employed at CPOM (Fig. 3) reads the input L1b files, then analyses and outputs one file at a time. The theory and methods described may act as a guide to developing a sea ice processing system for any Polar orbiting satellite radar altimeter.

Pre-processing of CryoSat-2 data
To begin, we set the latitude range within which northern hemisphere sea ice is found as 40°N-90°N. Any datapoints outside of this range are removed from the processing. Next, the surface type flag is used to remove all nonocean waveforms, and the MCD flag is used to remove all waveforms that are fatally degraded. Waveforms may be fatally degraded due to a computation error with the window delay or automatic gain control (AGC), or if the value for either parameter falls outside a specified range.
The window delay refers to the two-way delay time between pulse emission and the reference point at the centre of the range window, which for CryoSat-2 is centrally located; at bin 128 for LRM and SAR mode and bin 512 for SARIn mode. CryoSat-2 implements AGC to adjust the receiver sensitivity for the best reception of signals with varying return powers. AGC uses information about the previous return signal level to adjust the receiver gain in anticipation of the next. The aim is to keep the signal level as constant as possible (ESA/MSSL, 2013). Waveforms can also be fatally degraded due to inaccurate or missing information regarding the time that the radar pulse was reflected from the surface.
To allow for identical processing of both SAR and SARIn mode data acquired over Arctic sea ice, we crop all waveforms to 128 bins. As the location of the return waveforms within the range window is variable, especially in SARIn mode, we position each waveform at a consistent location within this cropped range window. We select bin numbers from b max À50 to b max +77, where b max is the range bin number corresponding to the maximum return power, and bins are counted from zero. The additional phase and coherence information available in the SARIn mode product is not required for sea ice processing. The individual SAR and SARIn mode files are then merged into single files for each Arctic crossing of the satellite, using the timestamps of the first and last waveforms in each individual file. Each timestamp is expressed as three integers; the day, seconds of the day, and microseconds. The orbital period of CryoSat-2 is extremely stable (Wingham et al., 2006), allowing the time of the ascending equator crossing before each Arctic pass and the descending equator crossing after each Arctic pass to be accurately computed.

Discrimination between sea ice and ocean waveforms
We discriminate between measurements of the ocean surface and the ice surface by identifying which echoes are specular and which are diffuse (Peacock and Laxon, 2004). Specular echoes occur when the radar burst is reflected from a smooth, mirror-like surface such as a lead or very thin ice. In these cases the power in the range window rises and falls again very rapidly, creating an echo that looks like a spike (Fig. 3a). Diffuse echoes occur when the radar burst is reflected from a rougher surface such as an ice floe or the open ocean. In these cases the power in the range window rises rapidly but gently decays, creating an echo that looks like a step (Fig. 3b).  (Warren et al., 1999). (a-d) Mean snow depth for January-March, April-June, July-September, and October-December, respectively. (e-h) Mean snow density for January-March, April-June, July-September, and October-December, respectively. Values are from the Warren climatology (Warren et al., 1999). The white polygon in (a) and (e) represents the central Arctic area within which the two-dimensional quadratic function fitted to all measurements is well constrained by in situ measurements. We identify specular and diffuse echoes by the values of their SSD and pulse peakiness. As described in Section 2, the SSD parameter uses the multi-looking capability of CryoSat-2 to provide a measure of the variation in surface backscatter with incidence angle for all stack waveforms that make-up the 20 Hz waveform, and is provided in the L1b product. The pulse peakiness is defined as: where PP is the pulse peakiness, P max is the maximum return power of the waveform, and P mean is the mean return power. Pulse peakiness is only calculated for bins where the return echo power is above the noise floor, where the noise floor is defined as the mean power in bins 10-19. The higher the peakiness, the more the echo looks like a spike. Specular echoes are defined as those with a pulse peakiness greater than 18 and a SSD less than 6.29 for SAR mode echoes and 4.62 for SARIn mode echoes. Diffuse echoes are those with a pulse peakiness less than 9 and a SSD greater than 6.29 for SAR mode echoes and 4.62 for SARIn mode echoes. Echoes with a peakiness or SSD between these values are considered complex and removed from the processing (see Fig. 4).
Next it is necessary to differentiate between radar echoes from ice floes and echoes from the ocean, as both get classified as diffuse according to the above definitions. We define ice floe regions as those with a sea ice concentration (Section 3.2.1) greater than 75%, and ocean regions as those with a concentration of 0%. Diffuse echoes from regions with an ice concentration between these values  are not trusted to come from either ice floes or open ocean, and are removed from the processing. We derived the floe concentration threshold empirically -the aim is for it to be high enough to avoid incorrectly classifying open water as an ice floe region. This could occur due to the relatively low resolution of the sea ice concentration data, which is provided on a 25 km grid.
During the sea ice melt season it becomes difficult to discriminate between measurements from leads and sea ice floes, due to melt ponds forming on the ice that cause specular return echoes. Eventually specular echoes dominate the majority of return waveforms. For example, by July and August over regions with a sea ice concentration greater than 0%, an average of 92% and 85% of echoes are classified as lead returns, respectively. In contrast, less than 1% of the waveforms are classified as ice floe returns. Therefore, we do not run the sea ice processor in the months of May-September. In the months that we process sea ice data (October-April), we discard an average of 34% of the waveforms over regions with a sea ice concentration greater than 0%, as they are not classified as a lead or ice floe return according to the criteria outlined in this section.

Definition of sea ice type
We use daily ice type data from OSI SAF (Section 3.2.2) to flag sea ice measurements as FYI, MYI, ambiguous, ice free, or unclassified. Records are removed where the ice type is ice-free, unclassified, or ambiguous. The ambiguous classification often relates to a thin band of ice between FYI and MYI regions. This band is highly mobile during the course of a month, so its removal does not cause gaps in monthly maps of sea ice thickness, and should not significantly affect monthly volume estimates. Maps of sea ice type, for the same day each year (Fig. 5), show interannual variation in the location of the FYI and MYI edges. But in all years the MYI cover is concentrated around the coast north of Greenland and the Canadian Archipelago, and often extends into the central Arctic.

Waveform retracking
For each echo, we select a specific point on the waveform leading edge to mark the location of the ocean or ice surface. This point often deviates from the reference point to which the window delay is measured, which for CryoSat-2 L1b data is the central range bin -bin 128 for SAR data and 512 for SARIn data. The purpose of this process, known as retracking, is to find the location within the full range window where the returned power comes from the surface at nadir. This value is returned as a bin number from the retracking routine and used later in the calculation of the surface elevation (see Section 4.5). We apply a different retracking routine depending on whether the surface return is specular (corresponding to a lead), or diffuse (corresponding to open ocean or a sea ice floe).

Specular echo retracking
We apply the Giles retracking method (Giles et al., 2007) to the specular echoes from leads. This method uses two functions to describe the shape of a specular echo. The leading edge of the echo is represented by a Gaussian function (where f = f 1 in Eq. (2)). The trailing edge of the echo is represented by an exponentially decaying function (where f = f 2 in Eq. (2)). These two functions are linked by a third linking function (where f = f L ). The full retracking function is: Where Yellow shading represents first year ice (FYI), red represents multi year ice (MYI), blue represents areas where sea ice is not present, grey represents an ambiguous ice type, and pink represents areas where the sea ice type has not been set in the data, which correspond to areas of land cover.
where P r t ð Þ is the power at time t, a is the maximum amplitude of the echo, t 0 is the time that P r ðtÞ ¼ a, t b is the time period for which f ¼ f L , k governs the rate of decay of the exponentially decaying function, and r is the standard deviation of the Gaussian function. We choose a 2 and a 3 such that the function, P r , and its first derivative, are smooth and continuous.
The retracking function is fit to each CryoSat-2 waveform using the Levenberg-Marquardt non-linear leastsquares method (Marquardt, 1963). The algorithm varies the fitting parameters a, r, t 0 and k to minimise the sum of the squares of the differences between the waveform samples and the fitted function, and reports it at the end of each iteration. The fit is run for a maximum of 3000 iterations and the retracking point is selected for the iteration where the sum-squared difference is at its minimum. The retracking point is the location of the maximum amplitude of the fitted function, t 0 . The retracker returns the decimal bin number, b 0 , corresponding to the retracking point, t 0 . If an acceptable fit is not achieved after 3000 iterations then the echo is removed from processing. In these cases the echoes differ in some way from the standard specular echo shape. For example, the lead dominating the return signal may be located away from the nadir of the satellite (Armitage and Davidson, 2014).

Diffuse echo retracking
Open ocean and sea ice floe waveforms are generally noisier than those originating from leads, especially in the case of SARIn data. Therefore, we smooth diffuse echoes using a three-point moving average before retracking. The surface height of the sea ice surface can be derived from the waveform leading edge, by determining the point along the edge that corresponds to the average surface. This will coincide with the two-way travel time for the midpoint of the pulse to reflect from the mean surface at nadir (Rapley et al., 1983). Therefore diffuse echoes are retracked using a threshold retracker. The tracking point is positioned where the leading edge rise of the echo first reaches 70% of the amplitude of its first peak. As the echo bin count will cross the 70% threshold somewhere between two bins, the exact tracking bin number, b 0 , is located by linear interpolation. The threshold is applied to the first peak rather than the peak of maximum amplitude, as the maximum may be caused by off-nadir ranging to leads (Armitage and Davidson, 2014). To be defined as the first peak, a peak must have an amplitude of 20% of the maximum, or higher. We retrack to the 70% power threshold, rather than the 50% threshold of conventional, pulse-limited altimetry, due to the synthetic aperture operation of CryoSat-2 over sea ice regions. In regular pulse-limited altimetry over a flat, homogenously rough, horizontal surface the surface area illuminated by a given radar pulse expands with time ( Fig. 6a) until the trailing edge of the pulse leaves the lowest reflecting point at nadir (Rapley et al., 1983). The leading edge of the return power waveform (Fig. 6b) displays a rise in power with time as the radar pulse illuminates an increasing area on the surface, and the rise in power is proportional to the illuminated area (Fu and Cazenave, 2000). The power of the return waveform then gradually decreases, as the radar pulse expands over the surface to form an annulus and a smaller surface area is illuminated. In pulse-limited altimetry the location of the mean surface height within the antenna footprint is taken to be the position where the leading edge rise of the echo first reaches 50% of the amplitude of maximum power, for surfaces with a Gaussian height distribution (returned power is proportional to illuminated area) such as sea ice Femenias et al., 1993). However, due to the synthetic, strip-like illumination pattern of each beam and the process of slant-range correction (Section 2), the point on each return waveform that corresponds to the surface does not lie at the half-power point of the leading edge (Wingham et al., 2006), but closer to 70% of the maximum.
Diffuse echoes are flagged as invalid in cases where the leading edge rise is complex, as these echoes typically produce anomalously low surface elevations. To detect echoes with a complex leading edge rise we measure the leading edge width, where the width of the leading edge is defined using: where W le is the width of the leading edge rise in number of bins, b 70 is the altimeter echo range bin number corresponding to the 70% threshold of the first peak, and b 30 is the altimeter echo range bin number corresponding to the 30% threshold of the first peak. Waveforms with leading edge widths larger than three bins are removed from processing (Fig. 7). Causes of such echoes over sea ice floes could include off-nadir reflection from leads (Armitage and Davidson, 2014), or reflection from a very rough surface (surface roughness $1.5 m and above), such as a heavily deformed or ridged ice.

Calculation of sea ice and ocean elevations
The next step is to compute the ocean surface elevations in the leads between sea ice floes, and the surface elevations of the sea ice floes. We assume that the radar pulses penetrate through any snow cover on ice floes and scatter from the snow-ice interface, which has been shown in laboratory experiments where the snow cover on sea ice is cold and dry (Beaven et al., 1995). Despite some evidence that the scattering horizon migrates as temperature rises (Willatt et al., 2010), we do not observe any bias in our thickness retrieval when compared to year-round ice draft data (Section 6.2 and Fig. 16), and so we conclude that the impact of this effect is not significant. The elevation of the lead or floe surface is computed using the following formula: where E is the elevation of the surface above the WGS84 reference ellipsoid, A is the altitude of the satellite centre of gravity (COG) above the same ellipsoid, and R 0 is the range of the satellite to the lead or floe surface. The altitude of the satellite COG is reported in the CryoSat-2 L1b data product. It is located at the centre of the single, spherical fuel tank (Wingham et al., 2006) and will not move during the satellite lifetime. R 0 is computed using the formula: where and where R n is the range of the satellite to the surface represented by the nominal tracking bin, C G represents the geophysical corrections that are applied, C R is the retracking correction, t n is the two-way travel time of the radar signal to the surface represented by the nominal tracking bin, b 0 is the bin number returned by the lead or floe retracker (Section 4.4), b n is the bin number of the nominal tracking bin, and D b is the bin width (0.2342 m). For SAR and SARIn data (Section 4.1) b n is centrally-located, at bin 128 and 512 of the range window, respectively. C G accounts for the effects of tides and atmospheric pressure on the surface elevations, and the effect of radar range delay due to propagation through the atmosphere. The corrections applied are taken from the CryoSat-2 L1b data product, and are: dry tropospheric, wet tropospheric, inverse barometric, modelled ionospheric, ocean tide, long period equilibrium tide, ocean loading tide, solid earth tide, and geocentric polar tide. Finally, C R accounts for the fact that lead and floe surfaces are located at bin b 0 (the bin number returned by the specular or diffuse echo retracker) rather than b n (the nominal tracking bin of CryoSat-2). Over sea ice, SARIn mode echoes are noisier than those from SAR mode, due to the lower burst-repetition frequency in SARIn mode (ESA/MSSL, 2013). In early iterations of our processing, sea ice elevations calculated from SARIn mode data were biased high as power drops on the waveform leading edge resulted in the retracker prematurely identifying the first peak of the waveform in the range window. By smoothing diffuse echoes using a threepoint moving average before retracking (Section 4.4.2) we are able to remove this bias. Individual SARIn elevation estimates will still be more sensitive to noise on the waveform leading edge than those from SAR mode, although this will have a minimal impact on our monthly gridded sea ice thickness product (Section 6.1), for which we average $400 individual measurements in each grid cell.

Removal of the mean sea surface
The calculation of sea ice freeboard requires us to know the instantaneous elevation of the ocean surface beneath sea ice floes, which can be obtained by interpolating between lead tie points (Section 4.8). However, at this stage in our processing the dominant signal in sea ice and ocean elevations will be due to the Earth's geoid and the mean circulation of the ocean currents. Together these make up the mean sea surface (MSS). The MSS is a static signal over the averaging period, which due to geoid gravity anomalies is of a higher frequency than time variant changes in the ocean surface beneath sea ice floes, so we remove it from all elevation estimates before continuing. For our processing we use the University College London 2013 (UCL13) MSS (Ridout, 2014). This was built using two years of Cryosat-2 Baseline-C data and improves on earlier MSS models in the region north of 81.5°N, where previously no satellite radar altimetry data was available. It is important to note however that once the MSS has been removed the remaining sea level anomaly (SLA) will still contain long wavelength errors in tides and atmospheric corrections.
After removing the MSS, we observe the occasional satellite track where ice and ocean elevations are consistently shifted by a few metres with respect to the MSS. This is physically unlikely, as the Arctic ocean dynamic topography varies on a scale of a few hundred kilometres with a maximum amplitude of around 0.5 m (Kwok and Morison, 2011). The shift can have a variety of causes, such as errors in the orbit determination or missing geophysical corrections. Shifted tracks are detected by computing the mean SLA for each track, from the individual lead elevation measurements along it. Tracks are removed where the mean SLA is lower than À0.5 m or greater than 0.5 m (Fig. 8). Before the mean is calculated, individual lead measurements with a SLA more than 20 m or less than À20 m are removed, as they are regarded as noise spikes. However, out of 35,000 satellite tracks, only 41 were filtered in this way when using Baseline-C data.

Retracker bias
As different retracking methods are used for specular and diffuse waveforms (Section 4.4), a fixed bias is introduced between elevation estimates from leads and the sea ice surface, which will affect subsequent freeboard (Section 4.8) and thickness (Section 4.9) estimates. This bias arises because the tracking points on specular and diffuse ocean waveforms are approximated by different theoretical returns. Therefore for one or both types of waveform, the approximation may differ slightly from the true tracking point, leading to a range bias between the two retracking techniques. To investigate the nature of the bias, we compared ocean elevation estimates using the two retracking methods in the Hudson Bay, which is a region of seasonal ice cover. During June and July there in no MYI cover in the Hudson Bay and the rate of ice retreat is rapid. Therefore in some places specular echo returns from very thin ice are interspersed with diffuse echo returns from areas of open water, which are similar to those from ice floes. We selected two satellite ground tracks from June and July 2011, where the SLA profile was relatively flat (no clear influence of tides or dynamic topography), and low-noise sequences of open water returns were interspersed with low-noise sequences of thin, flat ice returns. The thin ice echoes were then retracked using the specular echo retracker (Section 4.4.1) and open water echoes using the diffuse echo retracker (Section 4.4.2), to produce SLA profiles for the two tracks (e.g. Fig. 9). SLAs estimated using the diffuse echo retracker are positively biased with respect to those estimated using the specular echo retracker. To compute the bias we performed a straight-line fit through the SLAs from the specular echo retracker and calculated the difference between the line and each SLA output from the diffuse echo retracker. The mean value of the difference over both ground tracks, and therefore the bias, was 16.26 cm. We deduct this value from all elevation estimates returned by the diffuse echo retracker. It could be beneficial to repeat the analysis of retracker bias in different regions or at a different time of year, but it requires very specific conditions that are not common in regions of Arctic sea ice cover.

Calculation of sea ice freeboard
Sea ice freeboard is the elevation of the ice surface above the ocean. Before the calculation of sea ice freeboard, all lead measurements with a SLA of more than ±3 m are removed from processing. SLA values outside of this range are likely to be caused by noise in the retracking step. For example, there might be a second peak in the specular waveform that is selected by the retracking algorithm, rather than the first peak. The majority of SLAs lie within the smaller range of ±0.5 m. Of the 57 million lead measurements collected between November 2010 and April 2017, 99.8% of the resulting SLAs lay within ±3 m and 97.1% within ±0.5 m (Fig. 10).
We calculate sea ice freeboard for each waveform classed as containing an ice floe, or ice floes, following surface type discrimination (Section 4.2) and waveform filtering (Section 4.4.2). Between November 2010 and April 2017 this was 84 million waveforms. To calculate sea ice freeboard, we subtract the ocean surface elevation beneath the ice floe from the floe elevation. The ocean surface elevation beneath sea ice floes is calculated by interpolating the ocean surface elevation between leads, by using linear regression to perform a fit to the lead elevations extending 100 km either side of each floe location. The 100 km interpolation scale was chosen to ensure that ocean elevations were sufficiently interpolated without over-smoothing. At least one lead must be present on either side of the floe for the interpolation to be valid. This reduced the number of waveforms to 77 million.
We apply a correction to the calculated freeboard to account for the reduced propagation speed of light through the snow cover on sea ice floes. The corrected freeboard is given by: where f c is the corrected sea ice freeboard, f i is the original sea ice freeboard, h s is the snow depth (Section 3.2.4), c v is Therefore, corrected freeboard is simply expressed as: Corrected freeboard values outside the range À0.3 m to 3.0 m are removed from processing. These limits were selected by analysing a histogram of all freeboards between November 2010 and April 2017 (Fig. 11). Large negative freeboard values are likely to be caused by errors in the floe retracking and are removed, but slightly negative freeboards are permitted to allow for random noise in the returns from thin ice floes and ensure that the average freeboard is not biased high. Less than 1% of freeboard values are removed according to these criteria. The upper limit of 3 m was chosen as a means to remove extreme outliers, as the diffuse echo retracker occasionally returns anomalously large freeboard values. These values are likely a consequence of complex echoes whose leading edge width did not meet the removal criteria at the earlier stage of processing (Section 4.4.2), or of poor reconstruction of the ocean surface elevation through interpolation. Poor reconstruction of the ocean surface becomes more of an issue close to land, where fewer lead measurements are available. However, outliers of this magnitude are uncommon, and no freeboard measurements exceeded 4 m from November 2010 to April 2017. Over this period the mean freeboard was 20.50 cm with a standard deviation of 17.66 cm.

Calculation of sea ice thickness
Sea ice freeboard is converted to sea ice thickness by assuming that the ice floes are floating in hydrostatic equilibrium (Laxon et al., 2003) (Fig. 12). This means that sea ice thickness can be calculated using:  where h i is the sea ice thickness, f c is the corrected sea ice freeboard, q w is seawater density (1023.9 kg m À3 ) (Wadhams et al., 1992), h s is the snow depth, q s is the snow density (Section 3.2.4), and q i is sea ice density. We use a fixed sea ice density of 916.7 kg m À3 for FYI and 882.0 kg m À3 for MYI (Alexandrov et al., 2010), to account for the higher fraction of air-filled pores in MYI (Wadhams, 2000).

Calculation of sea ice volume
To calculate sea ice volume, each individual ice thickness measurement is assigned a corresponding ice concentration value from NSIDC concentration data, which is provided on a 25 km square grid (Section 3.2.1). The concentration assigned is the value from the NSIDC concentration grid cell that is nearest to the location of the thickness measurement, converted to a fractional ice concentration. Sea ice volume is calculated monthly, so we average all CryoSat-2 thicknesses and their corresponding concentrations on to a 0.5°longitude by 0.1°latitude grid for each month. The grid cells are considered empty if they contain less than five ice thickness measurements, and are filled using a nearest neighbour interpolation with a maximum search radius of 300 km.
Next we compute a sea ice extent mask, which is needed to define the sea ice edge, using NSIDC ice concentration data from the 15th day of each month. We use ice concentration from the 15th day of each month rather than a monthly average as rapid increases in concentration, particularly towards the end of the month, can bias concentration high at the ice edge. Each 0.5°by 0.1°grid cell is assigned the NSIDC concentration value that is nearest to its centre coordinates. If the concentration value from the 15th day is above 15%, then the cell falls within the ice extent mask and a value of one is set. A value of zero indicates that the cell is not within the extent mask. As a small number of grid cells will encompass land we produce an additional dataset that contains the fraction of ocean in each cell. This is done using the Generic Mapping Tools (GMT) function grdlandmask (Wessel et al., 2013;Wessel and Smith, 1991). The final 'cell volume' for each measurement location on the grid is the product of the sea ice thickness, the fractional sea ice concentration, the cell area, the sea ice extent mask, and the fraction of the cell believed to be ocean. The sum of all filled cell volumes gives the total, Arctic-wide, sea ice volume.
We also compute sea ice volume for FYI and MYI separately, using the ice type discrimination described in Section 4.3. In each grid cell the thicknesses of FYI and MYI are summed, and then used to compute the cell fraction of FYI and MYI, as follows: where F f is the cell fraction of FYI, F m is the cell fraction of MYI, h f is the total FYI thickness falling in the cell, and h m is the total MYI thickness falling in the cell. The FYI and MYI fractions in empty cells are filled using a nearest neighbour interpolation with a maximum search radius of 300 km. We use the FYI and MYI fractions to calculate FYI and MYI sea ice volume, by multiplying the total cell volume by the fractions.

Contributing factors
We estimate monthly errors in sea ice thickness and volume by considering the contributions due to uncertainties in snow depth (4.0-6.2 cm from Warren et al. (1999)), snow density (60.0-81.6 kg m À3 cm from Warren et al. (1999)), sea ice density (7.6 kg m À3 calculated below from Sever Expedition data), sea ice extent (20,000-30,000 km 2 according to NSIDC), sea ice concentration (5% according to NSIDC), and sea ice freeboard ($9 cm calculated below from CryoSat-2 freeboard observations). Uncertainties in seawater density are neglected because they have a negligi-  ble impact on uncertainties in sea ice thickness and volume Ricker et al., 2014).
We take snow depth and snow density uncertainties from a climatology (Warren et al., 1999), which was derived from fieldwork measurements acquired between 1954 and 1991 (Section 3.2.4). The climatology provides, as an error estimate, the standard deviation of snow depth and density anomalies in each calendar month. The anomalies were defined as the snow depth or density measured at the North Pole stations operating in that year, minus the multiyear average for the latitude and longitude of the station. If more than one station was operating, the anomalies were averaged. According to the authors, these errors are likely to be an overestimate, as the anomalies were calculated relative to measurements from only a few (typically 2) stations. Therefore to some extent the errors are recording regional anomalies that are accompanied by un-sampled anomalies of opposite sign in other regions in the same month. The uncertainty associated with the snow depth on Arctic sea ice increases throughout the ice growth season (Fig. 13), as snow accumulates and the associated standard deviation of depth anomaly increases. The uncertainty associated with the snow density remains more stable.
In the absence of extensive in situ sea ice density measurements, we estimate sea ice density (FYI and MYI) and its associated uncertainties from measurements acquired during the Sever expeditions of sea ice freeboard, sea ice thickness and snow depth on sea ice (Section 3.2.3). We use a rearranged version of Eq. (14) to calculate the ice density associated with each measurement location. For the ice density calculation, we set the densities of seawater and snow as 1025 kg m À3 and 324 kg m À3 , respectively, following the method of Alexandrov et al. (2010) whose densities of FYI and MYI are applied in our sea ice processing. We considered densities falling outside the range 860-970 kg m À3 unrealistic and they were discarded. We then calculated monthly average densities and discarded averages where fewer than four measurements were available. Unlike the snow climatology, average sea ice densities are not available for all months as the Sever expedition only ran in the spring. We therefore calculated the sea ice density uncertainty as the standard deviation of all available monthly averages, of which there were 18. This results in an uncertainty of 7.6 kg m À3 . This value is likely to be an overestimate of the true uncertainty due to undersampling, as was the case with the snow depth and density uncertainties.
NSIDC estimate sea ice extent as the region where its concentration exceeds 15%, and they estimate the relative (year-to-year) error as approximately 20,000-30,000 km 2 (http://nsidc.org/arcticseaicenews/faq/#error_bars) -a small fraction (0.1-0.5%) of the total extent. NSIDC quote a figure of 5% for the uncertainty in their sea ice concentration values (http://nsidc.org/data/docs/daac/nsidc0051_ gsfc_seaice.gd.html), which is important when considering local errors. As NSIDC do not estimate the distance over which the concentration uncertainty is correlated, we assume, conservatively, that the uncertainty is correlated over the entire Northern Hemisphere for each month.
Individual freeboard measurements have a standard deviation averaging 9 cm Arctic-wide. We calculated this by computing the standard deviation of all individual freeboard measurements within a 25 km radius at 25 km increments in each month from November 2010 to April 2017, then averaging over all months. The standard deviation arises through (a) uncertainties in the floe height measurement due to speckle in the radar echoes, which decorrelates from one measurement to the next, and (b) uncertainties in sea surface height, which may be correlated in space due to the interpolation scheme based on a linear regression of measurements along 200 km sections of each satellite ground track (Section 4.8). We assume that the principle source of uncertainty on an individual radar altimeter measurement of sea ice freeboard and consequent thickness is the speckle on the echo (Laxon et al., 2013;Wingham et al., 2006), which introduces noise in maps of sea ice freeboard and thickness. However, we output sea ice thickness data on to a user-friendly 5 km square Polar Stereographic grid by averaging all thickness measurements within a 25 km radius of the centre of each grid cell, with all points receiving equal weighting. Thickness uncertainties are calculated on the same grid, meaning that speckle error is reduced to a point where it no longer dominates uncertainty estimates. We chose to grid thickness data using a radius of 25 km as it is sufficient to reduce the gaps between the ground tracks at lower latitudes as well as suitably reduce the speckle error. Reducing the size of the averaging radius below 25 km does not reveal more detail in the maps, but does increase the noise. We then calculate the uncertainty associated with each grid cell thickness. To estimate the contribution of sea surface height uncertainty to freeboard uncertainty, we examined the variability of sea surface heights over the 200 km interpola- Fig. 13. Monthly mean snow depth from the Warren Climatology (Warren et al., 1999). The climatology is provided with error estimates, which are computed as the standard deviation of snow depth anomalies in each calendar month. tion length scale from November 2010 to April 2017, and their standard deviation at orbit crossing points was 4 cm. As a conservative estimate, we assume that this variability remains correlated within the 200 km window of each freeboard calculation, and include it as an additional source of uncertainty in our gridded thickness product. The freeboard error is then a combination of the 4 cm error due to spatially correlated errors in the interpolation of sea surface heights and that due to spatially uncorrelated errors in floe height measurement due to radar speckle. As the standard deviation on an individual freeboard measurement averages 9 cm Arctic-wide, and can be calculated through propagation of errors as the root-sum-square combination of the two sources of error (Ku, 1966), this leaves 8 cm for the floe height error. Although the number of floes heights averaged in each grid cell will be sufficient to reduce the 8 cm speckle error down to a negligible value, the 4 cm error in sea surface height will be reduced in the averaging only by the square root of the number of individual passes crossing the averaging window. When gridding monthly data this is typically 4 or more passes, resulting in a 2 cm freeboard uncertainty. This scales to $20 cm thickness, or 11% of a typical growth season thickness of 1.8 m  for gridded monthly thicknesses.

Method
To calculate uncertainties in sea ice volume, we calculate the monthly rate of change of volume with respect to each parameter that has an associated error. For most parameters, we adjust its value six times, at even increments, and re-compute volume each time. The computed rates of change are then multiplied by the error in each parameter in question to estimate their partial contributions to the total volume error (Table 1). Taking snow depth as an example, we compute the volume time series seven times, changing the snow depth on each freeboard measurement by À6 cm, À4 cm, À2 cm, 0 cm, +2 cm, +4 cm and +6 cm, to compute the monthly rate of change of volume per centimetre change in snow depth. This rate is multiplied by the monthly estimate of the snow depth error to estimate the contribution to error in sea ice volume.
To estimate the rate of change of sea ice volume with respect to sea ice extent, we recomputed sea ice volume using ice extent masks that use concentration data from the 10th day and the 20th day of each month, as well as the standard 15th day (Section 4.10). From these additional estimates, it is possible to compute the monthly rate of change of sea ice volume with respect to ice extent and hence assess the impact of this on volume error (Table 1). At 0.25% or less, the error in sea ice volume associated with year-to-year uncertainties in sea ice extent is insignificant. At sub-annual timescales, it is important to consider seasonal biases in sea ice extent when charting variability. During the period of sea ice freeze up, sea ice extent could be consistently underestimated by as much as 1 million km 2 (http://nsidc.org/arcticseaicenews/faq/#error_bars). Although the influence of this uncertainty on the volume error is not insignificant (Table 1), it does not affect yearto-year comparisons, and so we do not include it in our error budget, which is designed to illuminate uncertainties in inter-annual trends.
The contribution of sea ice concentration uncertainty to the total sea ice volume uncertainty is complicated, because we use the concentration data at two stages of our sea ice processing -to discriminate between radar echoes returning from ice floes and open water, and to weight the volume calculation according to the density of leads within the sea ice pack. Therefore, we do not calculate the monthly rate of change of sea ice volume with respect to sea ice concentration. Instead, we estimate the uncertainty in volume due to a 5% error in concentration. To do this the volume time series is recomputed twice for each month. In the first case we lower the sea ice concentration at every location by 5% and remove any ice floes where the concentration falls below the threshold of 75%. In the second case we raise the sea ice concentration at every location by 5%, but cap concentration at 100%. The monthly volume error is estimated as half the difference between these two recomputed volume time series.
Finally, we combine the monthly contributions to the volume error for all significant error sources in a rootsum-square manner to arrive at an estimate of the total monthly sea ice volume error, using: Table 1 The sea ice volume error budget. The October and April error columns give a value for the Arctic-wide error, with respect to the mean value, for each significant error source. The October volume error and April volume error columns show the contribution of each source to the total estimated sea ice volume error. These are then combined in a root-sum-square manner to give an estimate of the total monthly sea ice volume error.

Factor
October where r V is the uncertainty in sea ice volume in a given month, V is sea ice volume, h s is Arctic-wide snow depth, r hs is the uncertainty in snow depth, q s is Arctic-wide snow density, r q s is the uncertainty in snow density, q i is Arcticwide ice density, r q i is the uncertainty in sea ice density, e i is sea ice extent, r ei is the uncertainty in sea ice extent, and r V c is the uncertainty in sea ice volume due to uncertainty in sea ice concentration. We do not include a term for the contribution of sea ice freeboard uncertainty in Eq. (17) as there are typically more than 1 million floe height measurements and 10,000 200 km arc segments included in each monthly volume calculation, so the impact of both of these errors will be negligible on the monthly volume uncertainty.
Year-to-year, uncertainties in Arctic-wide sea ice volume are typically about 13.5%, with small variations from month to month (Table 1). Estimating local errors in sea ice thickness is complicated due to a lack of knowledge of the distances over which the contributing factors de-correlate. The main factors for which this information is important and lacking are snow depth, snow density, and sea ice density. In the sea ice volume error budget, their uncertainty is estimated over large scales as the standard deviation of monthlyaveraged sparse field observations collected across the $9 million km 2 central Arctic region. However, we assume that these factors, and their variability, are influenced by synoptic-scale meteorology, so we estimate the length scale over which they are correlated to be comparable to that of a typical polar vortex -around 2000 km in diameter (http://www.cpc.ncep.noaa.gov/products/stratosphere/polar/polar.shtml). Taking snow depth as an example, over areas that are large in comparison to this correlation scale, the variability of spatially averaged snowfall fluctuations will diminish in the ratio 1= ffiffi ffi n p . Here, n is the effective number of independent values of accumulation sampled, and is calculated as n $ A=ðp2000 2 Þ, where A is the area of the central Arctic (where field observations were collected) in square kilometres. If n < 1, we set it equal to 1. For the 9 million km 2 central Arctic region, over which the large scale sea ice volume and thickness uncertainty is estimated to be 13.5%, n $ 3, leading to an uncertainty of 23%. Using this approach, and accounting additionally for short-scale correlated errors in freeboard associated with interpolating sea surface heights, we estimate that the uncertainty in sea ice thickness increases to 25% at the 5 km scale of gridded monthly thicknesses. This is only a first attempt to characterise local uncertainty in sea ice thickness, and more detailed observations of snow depth, snow density, and sea ice density are required to establish the extent to which their variability impacts on the retrieval accuracy. However, a 25% local error in gridded monthly estimates of Arctic sea ice thick-ness derived from CryoSat-2 observations corresponds to an uncertainty of 45 cm for a typical thickness of 1.8 m. This uncertainty is consistent with the spread of differences relative to independent estimates acquired from airborne and ocean-based platforms, of 34-66 cm (Section 6.2).

Arctic sea ice thickness results
Using Baseline-C CryoSat-2 data, we are able to produce monthly sea ice thickness estimates for seven complete seasons of sea ice growth (October-April) from 2010-2017, except for one month in October 2010. The thickness data are output on to a 5 km square Polar Stereographic grid (Section 5.1). We do not run the sea ice processor in the months of May-September as melt ponds make it difficult to discriminate between radar returns from leads and from sea ice floes (Section 4.2). Sea ice thickness maps for autumn (October/November) 2016 (Fig. 14a) and spring (March/ April) 2017 (Fig. 14b) show that the thickest ice is concentrated around the coast north of Greenland and Ellesmere Island, and extends into the central Arctic as the growth season progresses. Over the CryoSat-2 period, the average autumn sea ice thickness in the Northern Hemisphere was 1.28 ± 0.13 m, increasing to 1.94 ± 0.09 m in spring.

Evaluation of Arctic sea ice thickness
To evaluate the accuracy of our CryoSat-2 sea ice thickness estimates, we have previously compared our Baseline-B results to independent estimates of sea ice thickness and draft acquired from airborne and ocean-based platforms. We used 772,090 estimates of sea ice thicknesses derived from NASA Operation IceBridge (OIB) airborne radar and laser altimeter measurements , 430 estimates of sea ice plus snow thickness derived by  airborne laser altimeter and electromagnetic (EM) sounding measurements (Haas et al., 2009), and 80 million estimates of sea ice draft derived from upward looking sonar (ULS) measurements as part of the Beaufort Gyre Exploration Program (BGEP) based at the Woods Hole Oceanographic Institution (http://www.whoi.edu/beaufortgyre). The OIB measurements were collected over the period 2011-2013 (L4 IDCSI4 product) and 2014 (quicklook product), the CryoVEx data over , and the BGEP data over 2010 Both the OIB and Cryo-VEx measurements were only collected during spring but they survey a variety of sea ice thickness and type ( Fig. 15a and b, respectively), and the OIB data cover a significant fraction of the western Arctic. The BGEP observations are available year round, but sample a restricted distribution of ice as the moorings are fixed (Fig. 15c). Like CryoSat-2 estimates of sea ice thickness, the OIB and Cryo-VEx estimates of sea ice thickness and sea ice plus snow thickness, respectively, are derived products. It is necessary to use these derived campaign products for evaluation due to the limited spatial and temporal coverage of in situ measurements of sea ice freeboard and thickness. The campaign products have improved spatial resolution compared to satellite data.
To make our CryoSat-2 sea ice thickness estimates equivalent to CryoVEX estimates of ice plus snow thickness, we applied the mean climatological snow depth (Section 3.2.4) to each CryoSat-2 thickness value. CryoSat-2 thicknesses were then compared to the OIB and CryoVEx data by gridding all datasets onto a 0.4°latitude by 4°longitude grid, which resulted in 1110 distinct OIB values and 64 CryoVEx values. To compare CryoSat-2 estimates to BGEP data, we calculated sea ice draft from CryoSat-2 freeboard measurements using: where d i is the equivalent draft from CryoSat-2, f c is the corrected sea ice freeboard measured by CryoSat-2, q i is the sea ice density, h s is the snow depth, q s is the snow density and q w is seawater density (Section 4.9). Monthly averages of all CryoSat-2 draft estimates were then taken within 100 km of each BGEP mooring and compared with monthly averages of ice draft obtained by each mooring, resulting in 58 distinct values. Overall, our CryoSat-2 measurements agree with the OIB, CryoVEx, and BGEP measurements of sea ice thickness, sea ice plus snow thickness, and sea ice draft to within 0.5, 21.0, and 10.0 cm on average, respectively (Fig. 16a-c).
To assess the overall bias in the CryoSat-2 observations we computed sea ice thickness from the BGEP estimates of sea ice draft and the CryoVEx estimates of snow plus ice thickness, and compared these to CryoSat-2 thicknesses along with the OIB thickness estimates (Fig. 16d). We computed BGEP equivalent sea ice thicknesses, h i , using: where d i is the sea ice draft measured by BGEP ULS buoys, q w is seawater density, h s is the snow depth, q s is the snow density and q i is the sea ice density (values given in Section 4.9). To convert CryoVEx ice plus snow thickness to ice thickness, we removed the mean climatological snow depth (Section 3.2.4) from each gridded CryoVEx measurement. When combined, the average difference between the OIB, CryoVEx, and BGEP estimates of ice thickness and those derived from CryoSat-2 is 2 mm. Given that the standard deviations in spring and autumn ice thickness derived from CryoSat-2 are 28 cm and 19 cm, respectively (Section 6.1) we conclude that there is no significant bias in the satellite data. The standard deviations between the CryoSat-2 and the OIB, CryoVEx, and BGEP estimates (66, 55, and 34 cm, respectively) are comparable to the estimated accuracy of the independent measurements (40 cm, 10 cm, and 10 cm, respectively) (Farrell et al., 2012;Haas et al., 2009;Melling et al., 1995) and the satellite (13 cm) (Wingham et al., 2006). The absolute differences between ice thickness estimates derived from the satellite and independent observations may arise through uncertainties in either dataset. We intend to repeat this initial evaluation, which was done with CryoSat-2 Baseline-B sea ice thickness estimates, with our Baseline-C estimates.

Arctic sea ice volume results
There have been clear seasonal variations in the amount of Northern Hemisphere sea ice over the CryoSat-2 period (Fig. 17). For all years, the total volume of sea ice increases each month over a given growth season from October until March or April. In most years the total volume drops slightly from March to April due to the onset of summer melt. This is generally also true for the volume of FYI and MYI, although punctuated with more variability from month-to-month. The uncertainty associated with FYI, MYI and total ice volume increases through the growth season due to an increase in the snow depth uncertainty (Section 5.1 and Fig. 13).
The volume of Arctic sea ice cover has also undergone large inter-annual fluctuations (Fig. 17). We observe the  (Fig. 14), with ice being 17% thicker, on average, than the six-year mean. This retention of thick ice over the summer melt season was associated with a 5% drop in the number of days on which melting occurred, compared with the previous three years . This coincided with conditions more typical of the late 1990s. The sharp increase in sea ice volume after just one cool summer demonstrates the ability of Arctic sea ice to respond rapidly to a changing environment. The volume of the autumn FYI cover is much less variable than MYI. In spring, inter-annual variations in hemisphere-wide volume are less significant than in autumn for all ice types. For example, the autumn 2013 increase in sea ice volume was followed by a 6% increase in spring 2014 volume compared to the previous three-year average, but this was not significant.

The relationship between Arctic sea ice extent, thickness and volume
Increases in the extent and mean thickness of the Arctic sea ice cover both contribute to the increase in ice volume over the growth season ( Fig. 18a-c, respectively). For all years, the total ice extent increases month-by-month over the growth season from October to March, and drops slightly in April around the onset of summer melt. The pattern is the same for the extent of FYI, although the MYI extent experiences a slight decrease over all but one (2012-2013) growth season (Fig. 18a). This is because the extent of MYI in the Arctic is dependent on two key processes. The first is the ageing of FYI to MYI and the second is ice export out of the high Arctic to more southerly regions Kwok, 2004). The ageing of sea ice is restricted to the start of the sea ice growth season, whereas sea ice export continues throughout the year. Therefore, the MYI extent experiences a continuous decrease after October/November time (Ye et al., 2016). The mean thickness (Fig. 18b) of total and FYI increases from March to April each year. The same is true for MYI, except in 2012 and 2013 where mean thickness remained constant from March to April. This could be as Arctic sea ice continues to thicken through congelation growth after it has reached its maximum annual extent (Wadhams, 2000), although dynamics also influence the extent and mean thickness of the sea ice pack (Kwok and Cunningham, 2015). Net thinning of the sea ice cover can occur due to the divergence of sea ice, but this in turn provides new areas of open water for FYI growth. Net thickening of the ice pack can occur through ice convergence, without any thermodynamic growth (Kwok and Cunningham, 2015). Despite total and FYI thickening from March to April, their volume follows that of their extent -it increases each month over a given growth season from October to March, and drops slightly in April. The volume of MYI also increases over the growth season, although there is more variation from month to month. This suggests that seasonally, the total extent and volume of sea ice in the Arctic are dominated by fluctuations in the FYI extent and volume, and start to decrease at the onset of melt.
Inter-annually, changes in the MYI extent, thickness and volume have a significant impact on the total extent, thickness and volume of Arctic sea ice. For example, we observed the minimum total ice extent (Fig. 18a) over the CryoSat-2 period in autumn 2012, which coincided with the minimum MYI extent. Despite this record low in ice extent, the mean thickness (Fig. 18b) and volume (Fig. 18c) of total and MYI reached their minimum in autumn 2011. In spring 2014, although the total ice extent was the second lowest over the CryoSat-2 period for that time of year (after 2017), we observed the highest total ice volume. This was associated with a peak in total ice mean thickness driven by record highs in MYI extent, mean thickness and volume. These findings also demonstrate that changes in ice extent do not always result in proportionate changes in ice volume, and so thickness information is required to assess the true state of the sea ice pack.

Conclusions
By providing unparalleled coverage of the Northern Hemisphere, data from ESA's CryoSat-2 mission have allowed us to produce hemisphere-wide sea ice thickness and volume estimates since November 2010. The provision of sea ice thickness data from satellite is crucial in understanding how the ice pack as a whole is changing. Here we have provided an end-to-end, comprehensive description of the processing steps that we use at CPOM to obtain estimates of Arctic sea ice thickness and volume from CryoSat-2 data, along with a detailed analysis of the uncertainties associated with our retrieval and an evaluation of our sea ice thickness product. In theory, the method presented could be used to retrieve sea ice thicknesses in the Southern Hemisphere, and it provides the foundation to develop sea ice processing systems for other Polar orbiting satellite radar altimeters.
Ideally, our uncertainty analysis would provide an error for each point measurement of sea ice thickness rather than grid cell values. This is currently hampered by a lack of knowledge regarding the correlation length scales and temporal variations of contributing factors such as snow depth and density, and sea ice density. The main contributors to the uncertainty in our retrieval are snow depth and snow density -each month they contribute an average of 10% and 6% to the total estimated sea ice volume error, respectively. These contributions are likely to be an overestimate. Therefore, a crucial next step in Arctic sea ice research is to develop improved estimates of snow loading on sea ice. The uncertainty contributions of radar propagation speed and depth through the snow cover should also be investigated. Despite these unknowns, our CryoSat-2 sea ice thickness estimates have been shown to agree with independent estimates of sea ice thickness and draft acquired from airborne and ocean-based platforms to within 2 mm, on average -a difference that is much smaller than the accuracy of either dataset (10-40 cm). However, this initial evaluation with CryoSat-2 Baseline-B sea ice thickness estimates should be repeated with our Baseline-C estimates. It could be further improved by the provision of more independent data that are near coincident in space and time with CryoSat-2 measurements, such as CryoVEx and OIB underflights of CryoSat-2, and by accounting for the sea ice drift that may occur between the location and time of each measurement.
To further refine CryoSat-2 estimates and uncertainties of sea ice thickness and volume, we require an improved understanding of the radar propagation through the snow on Arctic sea ice and of return characteristics such as the dominant scattering horizon. These could be obtained through, for example, in situ snow and ground penetrating radar (GPR) measurements (Beckers et al., 2015;Willatt et al., 2010Willatt et al., , 2011 or dual-frequency radar altimeter studies. Using dual-frequency radar it has been shown that the effective scattering horizon of a Ka-band satellite radar altimeter is higher in the snow pack than that of the Kuband CryoSat-2 radar (Armitage and Ridout, 2015;Guerreiro et al., 2016). However, these studies are hindered by the limited spatial and temporal coincidence of the two satellite datasets and the difference in their footprint size. Dual-frequency radar returns could also provide a muchneeded insight into the interaction of radar signals with the more complex snow and ice regimes in the Antarctic (Massom and Lubin, 2006;Massom et al., 2001;Maksym and Jeffries, 2000;Willatt et al., 2010;Schwegmann et al., 2016). The study highlights the potential benefit of future satellite missions with a multi-instrument payload. To understand how sea ice behaves basin-wide on decadal timescales, the continuation of satellite monitoring of the ice is crucial, as is the capability to continue and improve ice thickness retrievals.