An improved CryoSat-2 sea ice freeboard and thickness retrieval algorithm through the use of waveform fitting

Introduction Conclusions References


Introduction
Remote sensing records of Arctic sea ice thickness now span five decades and have shown nearly a two-fold decrease in mean winter thickness (Kwok and Rothrock, 2009) while observations over the past three decades have shown a 17.2 % decade −1 decline in the areal coverage of multiyear ice (Comiso, 2012).The inter-related decline in sea ice thickness and multiyear ice coverage is tied to declining trends in ice age and survivability (Maslanik et al., 2007;Maslanik et al., 2011).These changes have significant impacts to the climate, with a notable aspect of declining sea ice cover being linked to the observed higher than global average increase in Arctic surface air temperatures, a phenomenon known as Arctic amplification (Serreze et al., 2009).This occurs due to the increase of energy transferred from the atmosphere to the ocean as sea ice volume decreases (Kurtz et al., 2011;Rigor et al., 2002), which enhances warming and moistening of the lower troposphere (Boé et al., 2009;Screen et al., 2013).Changes in Arctic sea ice have also led to growing interest in determining predictability of the response of the sea ice cover to a changing climate.These interests range from efforts to improve short-term seasonal predictions (Lindsay et al., 2012;Eicken, 2013) to long-term predictions of when an ice-free summer may occur (Wang and Overland, 2012), and if ice-free summers can be sustained over the long-term (Tietsche et al., 2011).A key factor which links these disparate study areas is the need for continuous large-scale sea ice thickness observations to link physical processes to changes in sea ice and climate.
The earliest historical remote sensing record of Arctic sea ice thickness is composed of declassified submarine sonar observations extending back to 1958 (Rothrock et al., 1999).The submarine sonar sea ice thickness record is composed of numerous profiles within the central Arctic Ocean which need to be statistically analyzed to separate spatial, annual, and interannual variability within the limited regional data (Rothrock et al., 2008).Recent advances in satellite altimetry capabilities have enabled the deduction of sea ice thickness and volume over the larger-scale Arctic Ocean basin on monthly time-scales extending from the beginning of the growth season in October to the beginning of the melt season in May.Laxon et al. (2013) produced the first results of Arctic sea ice thickness from ERS-1 and ERS-2 satellite radar altimetry measurements spanning October 1993 to March 2001 up to the latitudinal limit of 81.5 ERS-1/2 radar altimetry record has also been extended using data from the Envisat satellite altimeter which showed large-scale thinning following the then record 2007 sea ice minimum (Giles et al., 2008).For the period spanning 2003-2008, data from the ICESat satellite laser altimetry mission provided a record of sea ice volume with increased coverage up to the latitudinal limit of 86 • .ICESat observed further decline in the thickness and volume of the Arctic sea ice cover in agreement with the radar altimetry results (Kwok et al., 2009).Presently, ESA's CryoSat-2 mission (Wingham et al., 2006), launched in 2010, is producing a continuous time series of radar altimeter measurements up to a latitudinal limit of 88 • , providing unparalleled coverage of the Arctic sea ice cover.Laxon et al. (2013) produced the first estimates of sea ice thickness and volume derived from CryoSat-2 data and validated the data with multiple in-situ data sets.The CryoSat-2 results were combined with ICESat estimates to produce the first decadalscale record of basin-wide Arctic sea ice volume from satellite altimetry.Data were also compared to estimates from the Pan-Arctic Ice-Ocean Modeling and Assimilation System (PIOMAS) model, which has shown volume loss of nearly 3×10 3 km 3 decade −1 from 1979-2010 (Schweiger et al., 2011) and similar trends from 1979 to the present.The combined ICESat and CryoSat-2 time series of sea ice volume change provides a useful tool to assess the PIOMAS data set which shows a loss of sea ice volume over a much longer time period.
With the advent of sea ice volume records from different satellite altimetry data sources comes the need to reconcile the assumptions used in the retrieval processes to produce a continuous time series and quantify uncertainties.Differences in sea ice thickness estimates from altimetry data arise in particular to the use of different density values and snow depth estimates which are used in the retrieval of sea ice thickness.
These quantities are due to environmental processes and should be applied in a consistent manner in the retrieval of sea ice thickness regardless of which instrument is used.In the case of sea ice density, previous studies have utilized a wide range of values, which will result in large differences between data sets if the same physical Figures assumptions are used.For example, in the study by Kwok et al. (2009) an ice-density of 925 kg m −3 was used, Kurtz et al. (2011) used a value of 915 kg m −3 , while Laxon et al. (2013) used an estimate of 917 kg m −3 for first year ice and 882 kg m −3 for multiyear ice.In these studies, the range of sea ice density values for multiyear ice is particularly large at 43 kg m −3 .For a typical multiyear sea ice floe with 60 cm of snowice freeboard and 35 cm of snow depth, the sea ice thickness estimate differs by 1.1 m within this range of ice densities.Despite the large-scale mean agreement of the sea ice thickness data sets described in previous studies, this discrepancy in physical assumptions points to the source of the differences as being due to potential biases in the freeboard and snow depth data sets used.This large discrepancy underscores the need to establish a set of consistent physical constants for use in the retrieval of sea ice thickness from satellite radar and laser altimetry data.
The focus of this study is to develop a new method for the retrieval of sea ice freeboard from CryoSat-2 data.We demonstrate that this method is consistent with independent measurements from airborne laser and radar altimetry data sets from NASA's Operation IceBridge mission to retrieve sea ice thickness which eliminates the need to utilize different ice density and snow depth values as an effective bias correction.The study is organized as follows: Sect. 2 describes the data sets used in the study.Section 3 describes the physical model which is used to simulate CryoSat-2 returns from Arctic sea ice covered regions.A procedure to utilize the model to fit CryoSat-2 waveforms for the retrieval of surface elevation is developed in Sect. 4. The new retrieval procedure is used to estimate Arctic sea ice freeboard and the results are compared to a threshold tracking method and independent freeboard observations from airborne data in Sect. 5.The results are then summarized and future improvements to the retrieval method are described in Sect.6. Figures

Back Close
Full resolution in vacuo).The satellite operates in SARIn mode over a spatially limited section of the Arctic Ocean, SARIn mode utilizes dual receive antennas to obtain phase information which can be used to detect the angle of off-nadir reflections.The focus of this study is to describe retrieval methods which can be used for the power detected waveforms.Thus, in order to maintain consistency in the retrieval algorithms developed here, phase information is not used and the SARIn data are truncated from 512 to 128 range bins.
The window delay field in the level 1B data provides the one-way travel time from the center range gate to the satellite's center of mass.We use this data to retrieve an elevation above the WGS84 ellipsoid by multiplying by the speed of light in vacuum, applying geophysical and retracking corrections, and subtracting these from the satellite center of mass altitude.Corrections for the elevation are given in the data products for the wet and dry tropospheric delay time, ionospheric delay, oscillator drift, inverse barometer effect, dynamic atmospheric correction, ocean equilibrium tide, long period ocean tide, load tide, solid earth tide, and pole tide.These corrections have been applied to the data used in this study.Retracking the mean scattering surface within the radar waveform is the focus of this study which allows for the sea ice surface elevation to be determined.Figures

Back Close
Full Data from NASA's Operation IceBridge airborne mission are used for comparison with monthly mean CryoSat-2 data for three campaign periods spanning March 2011 to March 2013.We use data from the IceBridge Sea Ice Freeboard, Snow Depth, and Thickness products (Kurtz et al., 2012) (Kurtz, 2013), but it is possible that the uncertainties in this data set are higher than in the final archival product (Kurtz et al., 2013a).The data consist of measurements from the Airborne Topographic Mapper laser altimeter (Krabill, 2010), Digital Mapping System camera (Dominguez, 2010), and the University of Kansas' 2-8 GHz snow radar (Leuschen, 2010;Panzer et al., 2013).Data from the individual instruments have been synthesized to provide sea ice freeboard, thickness, and snow depth at a 40 m spatial sampling resolution along all available flight lines using the methodology described in Kurtz et al. (2013b).Uncertainty estimates are also provided with the data products, which are estimated from the number of sea surface height tie points, distance to the local sea surface tie points, and the estimated covariance of the sea surface height for each flight.In this study we restrict data usage to where the uncertainty in the laser altimeter derived sea ice freeboard is less than 0.1 m.A map of the IceBridge derived ice freeboard used in the study is shown in Fig. 2.

CryoSat-2 multi-look echo phenomenology
In this section the behavior of the CryoSat-2 waveforms over surface types encountered in sea ice covered regions of the Arctic are simulated through the use of a physical model.The model shows the theoretical variation of the echo tracking point needed for the retrieval of surface elevation from the different surface types encountered.Before describing the model used in the simulation of CryoSat-2 returns, we acknowledge that due to the inherent complexity of scattering from sea ice covered regions assumptions Figures need to be made to simplify the problem to attain a tractable solution.In particular, we treat the scattering from sea ice as a surface problem.We furthermore assume the height deviations within the radar footprint are Gaussian and have a spatial exponential autocorrelation.Where appropriate, we note in the text where assumptions have been made and attempt to justify them, though there is of course a limit to the accuracy of the assumptions used.Given the assumptions made in the scattering model, it must still be treated as empirical, the validity of which is thus based on the degree to which it is able to model the phenomena of CryoSat-2 returns.Towards this end, potential improvements to the physics of the model are discussed in Sect.6.

Physical model for CryoSat-2 echoes
Here we provide the theoretical basis for modeling the mean echo power cross product from CryoSat-2 SAR and SARIn mode waveforms over sea ice covered regions of the Arctic.CryoSat-2 differs from previous generation pulse-limited radar altimeters (e.g.Envisat, Jason-1/2) largely in two ways: (1) the radar altimeter of CryoSat-2 consists of two antennas which have been narrowed in the across-track dimension to fit within the launcher fairing, thus it has an elliptical rather than circular antenna pattern which alters the impulse response (Wingham and Wallis, 2010).
(2) Unfocussed aperture synthesis is employed to reduce the along-track footprint size of the surface return.The level 1b data products are the result of a beam formation process which sums phase weighted and slant range corrected echoes taken from different look angles (see Wingham et al., 2006 for details).
The received radar echo, Ψ(τ), from a uniformly backscattering planar surface can be expressed as (e.g.Brown, 1977;Raney, 1998;Wingham et al., 2004) the convolution of the compressed transmit pulse after signal processing, P t (τ), the surface height probability density function, p(τ), and the "rough surface" impulse response (a factor of the surface geometry and antenna pattern), I(τ), (1) Figures

Back Close
Full where ⊗ represents convolution and τ represents the time delay relative to the time of the first surface arrival (τ = 0 is thus the point of closest arrival, which is here considered to correspond to the mean scattering surface).The use of Eq. (1) assumes that only surface scattering from the snow-ice interface is present (i.e. the surface is assumed to be perfectly conducting), surface scattering from the snow-air interface and volume scattering from within the snow and ice layers are neglected.This assumption is justified when the dominant reflection of energy occurs from the snow-ice interface i.e. when the density of the snow pack is relatively low (compared to ice), the reflection coefficient of the snow-air interface is much less than that of the snow-ice interface, and volume scattering within the ice layer is low.This will occur in practice when the snow pack does not contain wet or icy layers, and the roughnesses are not sufficiently different so as to cause a significant increase in the ratio of the Fresnel reflectivities of the snow-air and snow-ice interfaces.The CryoSat-2 compressed transmit pulse is well-represented by a sinc function described in Galin et al. (2013) as where p 0 is the peak power of the compressed pulse and B w is the received bandwidth (B w = 320 MHz).
The surface height probability density function is not known a priori, rather it must be determined through analysis of the waveform shape as will be shown in Sect. 4.Here we assume that p(τ) follows a Gaussian distribution given by of the snow-ice interface will take over the radar footprint.Future research in this area is needed to determine whether a different height distribution assumption can be used to improve the accuracy of the retrievals.Following Brown, (1977), the impulse response used to determine the power for a conventional altimeter is where λ is the center wavelength, g is the antenna pattern, σ 0 is the radar backscattering coefficient, t is the time from the instant of transmission, r is the range from the radar to the elemental scattering area dA on the surface.The angular component of the antenna gain pattern, (Θ, Ω), is measured relative to the antenna boresight, while the angular component of the backscatter dependence, (Γ, χ ) is relative to the surface normal.To simplify the problem, we here make the approximation that the satellite pitch and roll are zero and that the surface normal is parallel to the nadir direction.This allows for the angular components to be written as (Θ, Ω) = (Γ, χ ) = (γ, ω) where γ and ω are respectively the polar and azimuth angles subtended at the altimeter between the antenna boresight and scattering element.The standard deviations of the pitch and roll values over Arctic sea ice regions are small at 0.006 • and 0.01 • , respectively, with some of the observed variability due to noise in the star tracker measurements (Galin et al., 2013).The recorded mean pitch and roll over Arctic sea ice regions is less than 0.01 • , however, there is a known bias in the recorded pitch and roll values due to an error in the star tracker rotation matrices (Galin et al., 2013;Galin et al., 2014) which should be taken into account if a more physically exact characterization of the impulse response is desired.Following Wingham et al. (2006), Eq. ( 4) can be extended for application with beam k from nadir, ξ k .As described in Galin et al. (2013), the impulse response must also be summed over the different look angles used in the beam formation process.For the study by Galin et al. (2013) the look angle was defined in terms of the higher angular sampling in the burst which allows for the retrieval of the echo power and thus σ 0 .Since we are only concerned with the echo power shape, in the context of this study we define ξ k = k • 0.0238 which refers to the look angles from the stack data (the definitions of "burst" and "stack" are described in Wingham et al., 2006).Similarly, d 0 is the FFT of a Hamming window which is used in the formation of the mean echo cross product and D 0 is the FFT gain of the synthetic aperture minus the Hamming window loss.Using these definitions, the impulse response can then be written as We expand on the models of the CryoSat-2 impulse response described in Wingham et al. (2006) and Galin et al. (2013) by including a backscatter coefficient which varies with incidence angle, this will be shown to be necessary for modeling of the CryoSat-2 waveforms over sea ice.Over the range of incidence angles (up to 0.76 • ) encountered by CryoSat-2 for sea ice regions, the type of scattering is here assumed to remain within the specular scattering regime.For specular scattering, only surface facets which are tilted normal to the direction of the incident radiation contribute to the backscattering (Hagfors, 1964;Valenzuela, 1977).Hagfors, (1964) showed that for smoothly undulating surfaces with an exponential autocorrelation of height features the theoretical received power, Φ, with respect to incidence angle, φ, for a plane wave undergoing specular scattering is

Conclusions References
Tables Figures

Back Close
Full where h m is the rms height deviation, l is its length scale, and k 0 is the carrier wavenumber.We relate this to variations in σ 0 (γ, ω) by considering a scenario which varies only in the incidence angle by taking a ratio of the received power for nadir scattering and at an angle φ = γ (which considers only variations in the polar angle direction) which gives For the small incidence angles encountered by CryoSat-2, we assume cos 4 γ ≈ 1 which gives the approximate variation of backscatter with incidence angle to be where σ 0 (0 • ) is the backscattering coefficient at nadir and α = is a dimensionless variable that quantifies the efficiency of backscattering from a surface as a function of incidence angle.α is not known a priori and is determined from estimates of the waveform shape as described in Sect. 4. Note this treatment does not quantify σ 0 (0 and its dependence on the rms height deviation and length scale, it is here treated as an unknown constant which affects the echo amplitude, but not shape.Thus, the value of σ 0 (0 • ) is not retrieved in the context of this study.In choosing Eq. ( 8) a Gaussian height distribution with an exponential autocorrelation of height features is assumed since it allows for an explicit relationship between the surface height deviation and surface slope to be realized (Hagfors, 1964), which then allows for the simple mathematical characterization of the power directionality dependence described in Eq. ( 6).Physically, we may expect this assumption to be valid when the surface consists of a largely homogenous and isotropic field of height features.However, in areas such 732 Introduction

Conclusions References
Tables Figures

Back Close
Full as heavily ridged ice this assumption will likely introduce additional uncertainty in the results since the height distribution may not be Gaussian and the autocorrelation may have a different spatial dependence.Equation ( 5) can be reduced to a line integral around an isorange circle following the approximations to the scattering and geometry described in Wingham et al. (2006) and Wingham, (1995).The expression used here for the impulse response of the multilooked echo power shape follows the expression described in Galin et al. (2013) with the additional use of Eq. ( 8) to include the variation of backscatter with incidence angle.Consequently, the full expression for the impulse response is written as Tables Figures

Back Close
Full the contribution of the elliptical antenna pattern, it is taken from Wingham and Wallis (2010).The fourth line of Eq. ( 9) corresponds to the variation of backscatter with incidence angle.The fifth line corresponds to the gain of each synthetic beam and the application of a Hamming window.The form of the equation accounts for the slant range correction of each synthetic beam which is employed in the CryoSat level 1B data processor.

Waveform simulations
Equations ( 1), ( 2), (3), and ( 9) describe the physical model used to simulate CryoSat-2 waveforms over sea ice.As there is no closed form solution to these equations they must be calculated numerically.Here we detail the theoretical behavior of the waveform shape over both sea ice leads and floes by using this model.The free parameters in the model simulations shown here are σ and α. σ is varied from 0 to 0.4 m which represents the expected range from a smooth lead to ridged sea ice over the CryoSat-2 footprint.α is varied from 0 to 5 × 10 7 which represents the range from open ocean returns and very rough sea ice (where α ≈ 0; there is little to no backscatter dependence with incidence angle) to a perfectly smooth lead where backscatter from the nadir point dominates the echo.The leading edge of the waveform is affected by both α and the surface roughness.Increasing surface roughness increases the width of the echo, particularly from the rise time edge to the peak, it has much less impact on the trailing edge (Wingham et al., 2004) which is largely affected by the τ −1 2 behavior of the area of the range cells (Wingham et al., 2006).Equation ( 9) shows that as α becomes large, it dominates the decay of the trailing edge of the waveform.

Sea ice leads
The effect of variation in σ 0 with incidence angle is to decrease the effective illuminated area on the surface.This is apparent over sea ice leads where returns from geometrically small leads dominate the echo from radar altimeters (Drinkwater, 1991).shows simulations of CryoSat-2 waveforms for σ = 0.02 m and the observed range of α over sea ice leads (shown in Sect.5).It can be seen that placement of the tracking point to determine the surface elevation for sea ice leads is sensitive to α, because it determines the contribution of the off-nadir beams used in the retrieval.The mean scattering surface can be seen to correspond to the echo maxima for α 5 × 10 7 , and progressively moves toward the waveform leading edge as α decreases.
Mathematically, it can be shown that over smooth leads (σ ∼ 0; p(τ) = δ(τ)) with a suitably large value of α the impulse response function goes to a delta function, I(τ) = δ(τ, ξ k ).Note however, that Eq. ( 6) assumes incoherent reflections whereas scattering from a surface with σ = 0 will be coherent, which will affect the pulse amplitude, but not shape, which is the focus of this study.The received echo shape for a perfectly smooth lead will be which is simply a copy of the transmit pulse.This is also illustrated in the lowest amplitude waveform in Fig. 3 and can be seen in select CryoSat-2 waveforms over leads (an example of which is shown in Sect.4).As α decreases, returns from off-nadir are incorporated and the trailing edge of the waveform becomes longer.With the inclusion of more returns from off-nadir the mean scattering surface shifts leftward from the maximum peak power.Quantitatively, for σ = 0.02 m the echo peak corresponds to τ = 0.000 ns for α = 5 × 10 7 , and τ = 0.203 ns (0.030 m) for α = 5 × 10 5 , this range of 3 cm is the maximum sensitivity of elevation retrievals from leads due to α variations.

Sea ice floes
Simulated CryoSat-2 echoes from sea ice floes are shown in Fig. 4. For a Gaussian surface height distribution, Fig. 4 shows the mean scattering surface occurs when the leading edge reaches approximately ∼ 85-95 % of the peak value, with some variation of this threshold due to α and σ variations.The result that the retracking point for SAR echoes is near the peak, rather than at the half power point as is found in conventional Figures

Back Close
Full pulse-limited altimeters, was shown previously by Wingham et al. (2004).One point to note is the effective 1.563 ns sampling resolution of the instrument may not allow for the peak power to be well-determined for waveforms with low σ and high α, this will impact threshold algorithms which rely on a peak power ratio.For the 50 % threshold tracker used by Laxon et al. (2013), the simulations show biases due largely to variations in σ and less significant biases due to α variations.The biases range from −2.969 ns (leading to an elevation bias of +0.445 m) for σ = 0.4 m and α = 10 3 to −0.531 ns (+0.08 m) for σ = 0 and α = 10 5 .The variation of the τ = 0 point for different threshold values shown in Fig. 4b demonstrates that the freeboard for threshold tracking methods will likely be biased.However, the basin-wide bias encountered in an operational setting can not be accurately quantified from the simulations since it will be dependent on the combined effect of the surface roughness, the surface height distribution within the footprint, and the finite range resolution of the instrument.

Surface elevation retrieval algorithm
In this section, the physical model is combined with a least squares fitting procedure to estimate the mean scattering point and mean surface roughness within CryoSat-2 echoes from varying surface types.This least squares fitting procedure is analogous to routines which fit physical models to waveforms over ocean returns to retrieve surface elevation and other parameters such as significant wave height.Since Eq. ( 1) as developed here does not have a closed form solution, we describe the procedures which are used to fit the waveforms from numerical solutions.We show that through the use of look-up tables, the computation time of a least squares fitting routine is sufficient to fit the waveforms without the need for a closed form solution.Our fitting routine can fit a single CryoSat-2 L1B SAR/SARIn waveform on the order of one to ten seconds using a standard desktop computer, and a single month of CryoSat-2 data over Arctic sea ice can be processed in ∼ 10 days.Thus, the retracking method using the best model fit is practical from a processing standpoint.Figures

Back Close
Full

Fitting routine
In order to speed up calculation and enable fitting of individual waveforms, we calculate a look-up table of L(τ) = P t (τ) ⊗ I(τ) for a discrete set of cases encountered over Arctic sea ice and placed the data on an irregular grid.A fitting routine with pre-computed interpolation coefficients is then used to linearly interpolate between these discrete cases and quickly provide a solution for the function and its first and second derivatives (using the method of finite differences) for any queried point within the parameter space.
After creation of the look-up table, a least squares fitting routine using a bounded trust region Newton method (MATLAB function lsqcurvefit; described in Coleman and Li, 1996) is used to minimize the difference between the model fit and each CryoSat-2 echo power waveform, P r .A bounded trust region Newton method was chosen because the method is globally convergent, relatively independent of the problem size, and few iterations are needed to converge to a solution (Coleman and Li, 1994;Coleman and Li, 1996).Four free parameters are used in the fitting routine which is characterized by the equations where P m is the modeled waveform, P r is the observed echo power, and τ i corresponds to the observed echo power at point i of the waveform.The four free parameters are: (1) A f , the amplitude scale factor, (2) t the echo time shift factor, (3) α, (4) σ.
Given the dynamic range of the input parameters, and the fact that the solution which minimizes the square of the differences may not be physically correct, we specify an initial guess for each waveform and provide upper and lower bounds for the unknown parameters which are dictated by the physical system.The initial guess for A f is taken to be equal to the waveform peak power for all cases.For all other parameters, the Introduction

Conclusions References
Tables Figures

Back Close
Full methods for initial guess and upper and lower bounds are provided in the specific cases outlined below.

Leads
CryoSat-2 data over leads are identified in a similar manner to Laxon et al. (2013).First, a pulse peakiness parameter is calculated as Following Laxon et al. (2013), leads are defined as having a PP > 0.18 and a stack standard deviation < 4.An initial guess of σ = 0.02 m is used and the bounds are taken to be 0 ≤ σ ≤ 0.1 m. t i is first estimated to be the point of maximum power.An initial guess for α is estimated from the theoretical waveform peak to tail ratio which is taken from Sect. 3 and shown in Fig. 5.The tail is defined as the mean power of the six range bins (10 ns) following the point of peak power.The bounds for α are taken to be with a maximum uncertainty of 3 cm, but the uncertainty in surface elevation caused by errors in the choice of α in the fitting routing are likely small since the returns from leads can be seen to be very well-represented by the physical model.

Sea ice floes
For sea ice floes we follow Laxon et al. (2013) and define waveforms from floes as having a PP < 0.09 and a stack standard deviation greater than 4 (3 for SARIn mode regions).The initial guess for t i is taken from Laxon et al. (2013) as the first point where the waveform power reaches 50 % of the power of the first peak.The first peak is defined as the first local maximum on the waveform leading edge with a power value greater than 50 % of the point of highest power in the waveform.The waveform is only used when the power of the first peak is greater than 80 % of the highest power value in the waveform.The upper and lower bounds for t i are taken to be ±6 ns (±0.9 m) from the initial guess point.The initial guess for α is determined by the trailing edge of the waveform in a similar manner to that of leads.Figure 7 shows the tail to peak ratio which is used for sea ice floes.For sea ice floes, the tail is taken from the mean power of the set of measurements between 90 to 120 ns (58-78 range bins) after the point of peak power.The upper bound for α is taken to be 100 times the initial guess for α and the lower bound is taken to be the initial guess for α divided by 100.The trailing edge of the waveform is used for sea ice floes since the larger off-nadir angles experienced at larger delay times tends to eliminate the more "peaky" aspects from flat targets such as off-nadir leads.The initial guess for σ is set to 0.1 m, with a range of possible values from 0 to 1 m, if the initial guess for α is less than 8000 (which occurs over the open ocean and very rough sea ice floes) then the upper bound for σ is set to 6 m.
Example fits of the physical model to sea ice floes are shown in Fig. 8.The model fits the CryoSat-2 return waveform very well for both smooth (Fig. 8a) where σ = 0.05 m) and rough ice (Fig. 8b where σ = 0.34 m), which provides confidence in the ability of fitting model to be used to retrieve surface elevation over sea ice floes.We note that while Sect. 3 estimated the tracking point to be where the waveform leading edge reached 739 Introduction

Conclusions References
Tables Figures

Back Close
Full ∼ 80 % of the peak power, the fitted results shown in Fig. 8 demonstrate the finite range resolution of the instrument changes this value in practice such that a choice of an 80 % threshold would not yield accurate results in all cases.This is illustrated in Fig. 4b for the σ = 0 cases where the waveform leading edge, peak power point, and trailing edge are all located within the effective 1.563 ns sampling resolution of CryoSat-2, thus the peak power may not be adequately determined within a given waveform due to sampling limitations.
Model fits for areas with a mixture of smooth and rough surface types are shown in Figs.8c and d.The physical model developed in Sect. 3 assumed a surface with uniform characteristics which leads to the observed variations from the model fit.When the surface is not largely homogeneous within the CryoSat-2 footprint, a mixed return will result due to the different backscattering properties within the footprint.This is due to the inter-related variations in σ 0 , α, and σ which will combine to create a signal which has multiple peaks, unlike the single-peak smooth theoretical echoes seen in the model.The use of the pulse peakiness parameter to distinguish between sea ice floes and leads is discussed in Peacock andLaxon (2004), andLaxon (1994).In this study, we used the pulse peakiness and stack standard deviation thresholds used by Laxon et al. (2013) to minimize errors caused by mixed returns.The fitted returns show that if the smooth areas within the radar footprint have a large enough off-nadir angle so as to make the secondary peaks distinguishable from the main peak, then they do not largely impact the fitting routine since the location of the mean scattering surface is on the waveform leading edge.

CryoSat-2 derived sea ice properties
In this section we discuss the procedure for retrieving sea ice properties, including freeboard, roughness, and thickness from the CryoSat-2 data set.The freeboard and thickness results are compared to a threshold tracker method for sea ice floes and to independent measurements from NASA's Operation IceBridge campaign.

Sea ice property retrievals
For the frequency range used by CryoSat-2, the surface return from sea ice covered regions is assumed to be from the snow-ice interface as has been shown to be the dominant reflecting surface in laboratory experiments (Beaven et al., 1995), and field experiments when cold, dry snow is present (Willatt et al., 2011).Sea ice freeboard is thus defined here to be the height of the ice layer a.s.l. and is calculated as where h floe is the sea ice floe elevation and h ssh is the sea surface elevation.For all CryoSat-2 waveforms, we first removed the time varying sea surface height parameters outlined in Sect. 2. We then apply the retracking correction which is taken from the waveform fitting model used in Sect. 4. We calculate a monthly mean freeboard by gridding all sea ice floe and lead data points to a 25 km polar stereographic grid, with each grid point required to contain five floe elevations and five sea ice lead elevations to be flagged as containing a valid freeboard retrieval.After this initial gridding we then smooth the data by taking the average value for all points within ±2 grid points.This effectively reduces the spatial resolution to ∼ 100 km.A map of the mean gridded CryoSat-2 freeboard retrievals is shown in Fig. 9. Since the radar measures the ice freeboard which is the dominant factor in the retrieval of sea ice thickness, the spatial distribution of freeboard heights is expected to be similar to that of the ice thickness.
The map shows a spatial pattern which is consistent with past observations (Bourke and Garret, 1987;Kurtz et al., 2011;Kwok et al., 2009;Laxon et al., 2013) with high freeboards in the multiyear ice regions north of the Canadian Archipelago and Greenland, and lower freeboards in the first year ice regions of the Arctic Ocean and outlying seas.
The roughness of the scattering surface, σ, can also be retrieved from the waveform fitting method.A map of the surface roughness (excluding sea ice leads) is shown in from known dynamics and circulation patterns in the Arctic Ocean with the roughest ice corresponding to the multiyear ice area north of Greenland and the Canadian Archipelago.Gridded data points for σ and log 10 α are highly correlated with a correlation coefficient of −0.8, this demonstrates that, as expected, an increasing ice surface roughness corresponds to a lower angular variation in the radar backscatter coefficient.
In order to retrieve the ice freeboard needed for sea ice thickness retrievals, a geophysical correction to the CryoSat-2 freeboard must also added to account for variations of the speed of light within the snow pack on sea ice.This is given as fb = fb radar + h c , where the correction factor, h c , is given as where h s is the snow depth and c snow is the speed of light within the snow pack.c snow is parameterized following Tiuri et al. (1984) to be where ρ s is the density of snow with units of g cm −3 .For the comparison with IceBridge data discussed in the next section, this geophysical correction adds a mean value of 4.9 cm to fb radar to attain the true ice freeboard.Sea ice thickness, h i , can be retrieved from the CryoSat-2 data set through the assumption of hydrostatic balance where ρ w , ρ i , ρ s are the respective densities of sea water, ice, and snow.Thus, the retrieval of sea ice thickness requires an independent snow depth data set as well as assumptions of the density properties of the surface.In this study we use density Figures assumptions which are discussed in Kurtz et al. (2013b) to be consistent with the Ice-Bridge data.The density of sea ice is taken to be 915 kg m −3 , the density of snow is taken to be 320 kg m −3 , and the density of sea water is taken to be 1024 kg m −3 .Using these values Eq. ( 16) can be written as

Comparison of CryoSat-2 freeboard and thickness data
In this section we compare CryoSat-2 retrieved freeboard data using the waveform fitting method and an empirical lead and threshold floe tracker.We then provide an independent comparison to Operation IceBridge data.

Comparison to sea ice floe threshold tracker
In order to illustrate differences between the new freeboard retrieval method developed in this study, we compare freeboard retrievals from the waveform fitting retracker to a similar freeboard retrieval method outlined in Laxon et al. (2013).The method of Laxon et al. (2013) uses a 50 % threshold tracker in the retrieval of sea ice floe elevations, and retracks sea ice lead returns using an empirical fit function described in Giles et al. (2007).Our reproduction of a similar method is hereinafter referred to as the ELTF (empirical lead and threshold floe) retracker.We note that several differences are present between the freeboard retrieval used by Laxon et al. (2013) and the ELTF method used here.The primary difference is that Laxon et al. ( 2013) subtracted a bias from the sea ice lead elevations by taking the difference between returns from the ocean when sea ice is not present and returns from leads in the nearby ice pack.This was done following Giles et al. (2007), but was not done in the ELTF freeboard retrievals.Additional differences include (but are not limited to) the exact definition of the first peak, as well as the use of a mean sea surface height data set in place of the EGM08 geoid.Therefore, the comparisons done in this study are similar, but not

TCD Figures Back Close
Full exact reproductions of methodologies.The purpose of the comparison is to highlight the physical basis between differences in the retracking methods.The mean difference in sea ice lead elevations retrieved by the waveform fitting method described in this study and the empirical tracker described in Giles et al. (2007) is 2.8 cm (Giles et al., 2007 tracker-waveform fitting method), and the correlation is 0.7.The most significant difference difference between the freeboard retrieval method of Laxon et al. (2013) and the waveform fitting method is the use of the 50 % threshold tracker in the retrieval of sea ice floe elevations, which as illustrated in Sect. 3 is expected to be biased high from theoretical arguments since the selected threshold should be closer to the waveform peak.Figure 11 shows the retrieved freeboard using the ELTF method, and Fig. 12 shows the difference with the waveform fitting method.The mean freeboard differences are 11.9 cm, 12.7 cm, and 11.5 cm for the March 2011-2013 periods.This corresponds to mean ice thickness differences of 1.12 m, 1.19 m, and 1.08 m using Eq. ( 17).The differences shown in Fig. 12 also show significant spatial and interannual differences between the methods.The mean freeboard using the ELTF method for March 2011 is 31.3cm, using Eq. ( 17) this corresponds to a minimum mean sea ice thickness of 2.9 m, which will be higher once one considers the contribution of snow.The waveform fitting method gives a minimum sea ice thickness of 1.8 m, which is much closer to the mean thickness of first year ice which is now the dominant ice type in the Arctic (Comiso, 2012).The h c snow speed of light correction was not applied in the comparison between threshold and waveform retrackers because it is equivalent for both data sets, but this will slightly increase the mean minimum thickness.Thus, the higher freeboard values retrieved by the threshold method are likely biased high, which is in agreement with the theoretical arguments presented in Sect.3.

Comparison to IceBridge data
In order to compare the ELTF and waveform fitting methods, we now compare both methods to independent data collected from three measurement campaigns of the 744 Introduction

Conclusions References
Tables Figures

Back Close
Full Operation IceBridge mission.The mean IceBridge snow depth has been subtracted from the laser altimeter freeboard to determine the ice freeboard, and the data have then been gridded to the same 25 km polar stereographic grid as the CryoSat-2 data.A grid point is defined as containing valid data for comparison when there are greater than 200 IceBridge measurements and a valid gridded CryoSat-2 measurement.Since snow depth information is available from the IceBridge data set, we add the h c correction factor to the CryoSat-2 retrieved freeboards and also estimate sea ice thickness using Eq. ( 17).
Table 2 summarizes the comparison between the IceBridge observations, the waveform fitting method, and the ELTF method.For the waveform fitting method, the mean freeboard difference (CS2-IceBridge) ranges from 1-3 cm while the ice thickness difference ranges from 11-23 cm.The slightly higher freeboard retrieved by CryoSat-2 is consistent with the results of Armitage and Davidson, (2014) who estimate that the sea surface height will be biased low by ∼ 2 cm due to off-nadir ranging to leads when a minimum pulse peakiness of 0.18 is used as a threshold for lead classification.The mean freeboard difference (CS2-IceBridge) for the ELTF method ranges from 11.9-15.9cm which corresponds to ice thickness differences of 112-149 cm, this is significantly higher than the waveform fitting method.As shown in Sect.3, this is likely due to the choice of the 50 % threshold which was shown to be too low in comparison to theoretical estimates which show the tracking point should be closer to the peak power value.Surface roughness and the finite sampling resolution of the radar also plays a role as well.We note that Laxon et al. (2013) did not add a correction for the speed of light within the snow pack and also subtracted a constant value from the sea surface elevation due to the use of different fitting models between open ocean and leads.In the waveform fitting retrieval scheme illustrated in this study, no such bias in the sea surface height needs to be removed because the same model is used to fit waveforms from open ocean, sea ice floes, and leads.The addition of the snow speed of light correction will also apply equally to each method.Thus, the waveform fitting method Introduction

Conclusions References
Tables Figures

Back Close
Full gives a mean difference which compares much better to the IceBridge ice freeboard data using explicit geophysical arguments.The root-mean-square difference between the IceBridge data and CryoSat-2 freeboard retrievals ranges from 7.4-11.1 cm for the waveform fitting method and is higher at 14.1-19.8cm for the ELTF method.The mean estimated IceBridge freeboard uncertainty (taken from the data products using the method described in Kurtz et al., 2013) for the compared grid points is 5.9 cm, 7.6 cm, and 6.3 cm for the respective 2011-2013 campaigns, the uncertainty in the sea surface height is due to a combination of instrumental uncertainties and is also a function of distance to the nearest lead.For the waveform fitting method, the estimated uncertainty in the CryoSat-2 sea ice freeboard, σ cs2−fb , can be calculated as where σ diff is the observed standard deviation of differences and σ IceBridge is the uncertainty in the IceBridge ice freeboard.One complication with this estimate is that σ IceBridge is a set of values, rather than a constant number.However, using the mean value of σ IceBridge for each campaign gives an estimate for the CryoSat-2 freeboard uncertainty of 9.2 cm, 6.6 cm, and 3.8 cm for the respective 2011-2013 campaigns.
The correlation between the waveform fitting method and IceBridge data varies substantially between the campaigns, but the correlation between the waveform fitting method and IceBridge data is higher for all campaigns than the ELTF method.While the low correlations to the IceBridge data may be some cause for concern, we note that the reasonable spatial distribution of sea ice freeboard shown in Fig. 9 and estimated uncertainties for the IceBridge and CryoSat-2 freeboard retrievals place this into a context which can be understood.Given the uncertainties present in both the IceBridge and CryoSat-2 data and the small dynamic range of the freeboard values, it is possible a high correlation value can not be attained from the comparison.To test this hypothesis, a Monte Carlo simulation was conducted using the CryoSat-2 freeboard values.Two sets of numbers were constructed by adding a random distribution Introduction

Conclusions References
Tables Figures

Back Close
Full with zero mean and a standard deviation equivalent to the estimated mean uncertainties for the CryoSat-2 and IceBridge data.The correlation was computed for each set of numbers and the simulation was run 1000 times.The simulation shows the expected correlations for two identical data sets with estimated uncertainties equivalent for the respective 2011, 2012, and 2013 campaigns are 0.46±0.05,0.55±0.03,and 0.60±0.04.
The March 2013 data set shows a correlation which is consistent with the estimated uncertainties for the data.However, it is not clear why only this campaign shows a lower RMS difference and a correlation which is in line with expectations.A possible explanation is that additional variability due to the use of non-temporally coincident data sets added additional uncertainty to the estimates shown here, and that the lower bound uncertainty of 4 cm is correct for the CryoSat-2 freeboard retrievals.A more detailed comparison between time coincident IceBridge data flights which underflew CryoSat-2 will be the subject of a future study.

Conclusions
A new method to fit CryoSat-2 level 1B waveforms using an empirical model was developed.This waveform fitting procedure was used to retrieve sea ice freeboard from CryoSat-2 over Arctic sea ice.Through comparison with Operation IceBridge data for the 2011-2013 campaigns, this study has shown that fitting of the CryoSat-2 level 1B waveforms using a physical model can be used to obtain improved results over the empirical lead and threshold tracker (ELTF) methods which are similar to those used by Laxon et al. (2013).The ELTF method was found to have respective mean freeboard differences (CryoSat-2 -IceBridge) of 15.4 cm, 15.9 cm, and 11.9 cm and mean sea ice thickness differences of 144.2 cm, 149.3 cm, and 111.9 cm.The mean freeboard differences for the waveform fitting method were 2.2 cm, 2.5 cm, and 1.1 cm, and the mean sea ice thickness differences were 20.6 cm, 23.3 cm, and 10.6 cm.The larger RMS and mean differences in the ELTF tracker method were found to be largely due to the choice of the 50 % threshold, which was shown to be too low based on theoretical Introduction

Conclusions References
Tables Figures

Back Close
Full modeling.The difference is also due to variations in surface roughness and the angular dependence of the backscattering coefficient.A bias of 1.9 cm was found in the waveform fitting method freeboard retrievals compared to the IceBridge data, this bias is consistent with the estimated range bias due to off-nadir ranging of lead points shown by Armitage and Davidson (2014).A maximum correlation of 0.57 was found between the IceBridge freeboard and thickness data and the waveform fitting method, this correlation is consistent with an estimated uncertainty of 4 cm in the retrieved CryoSat-2 freeboard for a 100 km gridded data point.Despite having a physical basis, and having a small bias compared to the airborne observations, the model used to fit the CryoSat-2 waveforms is still essentially empirical.In order to move towards a more physically exact model a number of points need to be taken into account, which are largely due to the considerable variability of surface types (and their associated backscattering properties) which can be found within the radar footprint: (1) it was assumed that the distribution of surface heights within the footprint can be approximated with a Gaussian function, though it is possible that the presence of ridges will lead to a more skewed distribution.(2) The presence of ridges may also lead to an electromagnetic bias if the scattering from the ridge peaks is different than the surrounding ice, this is similar to a known phenomenon which has been observed in open ocean returns wherein wave troughs have a higher backscatter than wave crests (Yaplee et al., 1971).(3) The model assumes the antenna boresight is always at nadir and the surface normal is parallel to the nadir direction.However, recently discovered pitch and roll biases within the CryoSat-2 data mean that the antenna boresight is slightly off-nadir which should be taken into account.(4) Mixed returns containing more than one surface type are not dealt with in the model, and the fitting procedure only works when the mixture of surface types allows for sufficient separation between peaks to fit the return.The usage of additional statistics such as goodness of fit estimates may be used to further reduce errors caused by mixed returns in future studies.(5) Surface scattering from the snow-air interface, and volume scattering within the snow pack and sea ice were not considered.In particular, this may cause Introduction Full additional uncertainty if the snow cover contains dense layers or if the roughness of the snow-air interface is much lower than that of the snow-ice interface.Further improvements to the retrieval of surface elevation, sea ice freeboard, and sea ice thickness can also be done which should reduce the uncertainty in the measurements.Modeling of the phase and further analysis of the SARIn data areas in the Arctic Ocean may lead to ways to identify off-nadir sea ice leads and reduce the observed ∼ 2 cm bias which this causes in the sea surface height and freeboard data sets.In a similar manner to what was described in Laxon et al. (2013), the retrieval of a large volume of sea surface height estimates will allow for the construction of a high resolution mean sea surface height data set which can be used to reduce geoid errors which are known to be prevalent in the Arctic Ocean (McAdoo et al., 2013).Lastly, an evaluation of the IceBridge snow depth measurements needs to be done to improve basin-wide snow depth on sea ice estimates.This has been done for a single season of data (Kurtz and Farrell, 2011) compared to the snow depth climatology of Warren et al. (1999), and for passive microwave retrievals of snow depth on sea ice for first year ice (Brucker and Markus, 2013).The focus of a future study will be to utilize existing observations to improve estimates of snow depth on sea ice to be used in the retrieval of sea ice thickness from the CryoSat-2 time series.
Overall, this study has further demonstrated the capabilities of CryoSat-2 for the retrieval of sea ice freeboard and thickness.The advantage of the retrieval processes used in this study is that they are compatible with the laser altimetry record and show that the two records can be reconciled to produce a more complete time series of sea ice volume change.This has distinct advantages for the expected launch of the ICESat-2 laser altimeter mission in 2017.The lifetime of CryoSat-2 is expected to overlap with the ICESat-2 mission, as is the new Sentinel-3 radar altimeter mission.The combined satellite radar and laser altimetry data provided by these missions will thus provide unmatched information on the state of the Arctic sea ice cover.Introduction

Conclusions References
Tables Figures

Back Close
Full  Full  Full Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | The primary data set used in this study comes from ESA's CryoSat-2 satellite(Wingham et al., 2006).Data are taken from the baseline B Level 1B SAR and SARIn mode data products forMarch 2011March  , 2012March  , and 2013.Example CryoSat-2 SAR mode waveforms and terminology employed in the description of the waveform features are shown in Fig.1.CryoSat-2 is a radar altimeter which operates at a center frequency of 13.575 GHz and has a receive bandwidth of 320 MHz.The effective footprint size after post-processing is ∼ 1650 m in the across track direction and ∼ 380 m in the along track direction.The power detected echoes contain 128 range bins in SAR mode and 512 range bins in SARIn mode, each range bin is sampled at 1.563 ns (0.234 m range Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | σ c = 2σ/c which is the surface roughness in the time domain and c is the speed of light in vacuo.A Gaussian height distribution was chosen since it is dependent on only a single parameter, and it is not presently known what form the height distribution Discussion Paper | Discussion Paper | Discussion Paper | CryoSat-2 through the addition of a synthetic beam gain term, where d 0 (cos ζ − sin ξ k ) is the synthetic beam gain, which is a function of the angle between the direction of a scattering element and the satellite velocity vector, ζ , and the look angle of synthetic 730 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | 0 < 100α 0 where α 0 is the initial guess.Example fits to CryoSat-2 waveforms over leads are shown in Fig. 6.The figure shows the behavior of the CryoSat-2 waveform for increasing values of α.As shown in Sect.3.2.1, over smooth leads (σ ∼ 0) with a large value for α (α 5×10 7 ), the received waveform is simply a copy of the transmit pulse which may be slightly broadened by the small surface roughness within the lead.An example CryoSat-2 waveform showing this behavior can be seen in Fig. 6d.For all lead cases, the tracking point for the mean scattering surface is near the maximum peak of the return.For α < 5×10 7 , returns from off-nadir begin to broaden the waveform and shift the mean scattering surface leftward from the maximum peak due to the inclusion of off-nadir look angles i.e. it is determined by the combined effect of the impulse response for each off-nadir look angle and the slant range correction used in the data processor.As shown in Sect.3.1.1,the tracking point for the mean scattering surface over leads is thus sensitive to the choice of α Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Fig. 10.The sea ice floe roughness also corresponds well to what may be expected Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | teams and air crews for long hours in the field and at home collecting and processing the data and the National Snow and Ice Data Center for archiving and publishing the data.This work is funded by NASA's Airborne Science and Cryospheric Sciences Programs.Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

Fig. 1 .Fig. 2 .Fig. 3 .
Fig. 1.Example CryoSat-2 waveforms.(a) Example waveform for a sea ice floe.The waveform contains an off-nadir reflection from a surface with high backscatter which results in a secondary peak.(b) Example specular waveform from a lead.

Figure 4 .Fig. 4 .Fig. 5 .Fig. 6 .Fig. 7 .Figure 8 .Fig. 8 .Fig. 9 .Fig. 10 .Fig. 11 .Fig. 12 .
Figure 4. Simulated CryoSat-2 echoes over sea ice floes for a) the typical range of delay times provided in the level 1B product over sea ice, b) a zoomed in plot showing the behavior near the echo delay time of 0. Solid lines correspond to α = 10 3 and dashed lines correspond to α = 10 5 .39 Fig. 4. Simulated CryoSat-2 echoes over sea ice floes for (a) the typical range of delay times provided in the level 1B product over sea ice, (b) a zoomed in plot showing the behavior near the echo delay time of 0. Solid lines correspond to α = 10 3 and dashed lines correspond to α = 10 5 .
from 16-28 March 2011, and 14 March-10 April 2012.Data from the quick look data set have been used for 20-27 March 2013 since final data from the campaign are not yet available.The 2013 quick look data utilizes new processing techniques to minimize freeboard biases

Table 1 .
Summary of parameters and symbols used in the CryoSat-2 model.

Table 2 .
CryoSat-2 freeboard and thickness retrievals compared to IceBridge airborne data.Values pertaining to sea ice thickness are in parentheses.