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

We develop a physical model capable of simulating the mean echo power of CryoSat-2 SAR- and SARIn-mode waveforms over sea-ice-covered regions. The model simulations are used to show the importance of variations in the radar backscatter coefficient with incidence angle and surface roughness for the retrieval of surface elevation of both sea ice floes and leads. The physical model is used to fit CryoSat-2 waveforms to enable retrieval of surface elevation through the use of lookup tables and a bounded trust region Newton least-squares fitting approach. The use of a model to fit returns from sea ice regions offers advantages over currently used threshold retracking methods, which are here shown to be sensitive to the combined effect of bandwidth-limited range resolution and surface roughness variations. Laxon et al. (2013) have compared ice thickness results from CryoSat-2 and IceBridge, and found good agreement; however consistent assumptions about the snow depth and density of sea ice were not used in the comparisons. To address this issue, we directly compare ice freeboard and thickness retrievals from the waveform-fitting and threshold tracker methods of CryoSat-2 to Operation IceBridge data using a consistent set of parameterizations. The purpose of the comparison is to highlight the physical basis between differences in the retracking methods. For three IceBridge campaign periods from March 2011 to March 2013, mean differences (CryoSat-2 – IceBridge) of 0.144 and 1.351 m are found between the freeboard and thickness retrievals, respectively, using a 50% sea ice floe threshold retracker, while mean differences of 0.019 and 0.182 m are found when using the waveform-fitting method. This suggests the waveform-fitting technique is capable of better reconciling the sea ice thickness data record from laser and radar altimetry data sets through the usage of consistent physical assumptions.


Introduction
Remote sensing records of Arctic sea ice thickness now span five decades and have shown nearly a 2-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 interrelated 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 on 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 N. T. Kurtz et al.: CryoSat-2 sea ice freeboard 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 data coverage of the submarine cruise tracks (Rothrock et al., 2008).Recent advances in satellite altimetry capabilities have enabled the deduction of sea ice thickness and volume over the largerscale Arctic Ocean basin on monthly timescales extending from the beginning of the growth season in October to the beginning of the melt season in May.Laxon et al. (2003) 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 • .The 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][2004][2005][2006][2007][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 decadal-scale 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 to 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 from 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 pro-cesses 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 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 snow-ice freeboard and 35 cm of snow depth, the sea ice thickness estimate differs by 1.1 m within this range of ice densities.The uncertainty in sea ice thickness from the selection of ice density, as well as natural variability in ice density and snow depth, has implications for the uncertainty in temporal trends in Arctic sea ice thickness and volume estimates (Zygmuntowska et al., 2014).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.In the context of this study, we define freeboard as the height of the sea ice layer above the sea surface.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 freeboard and thickness.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. Section 6 estimates the errors which can be present in the waveform-fitting method when backscattering from the snow layer and sea ice volume are considered.The results are then summarized and future improvements to the retrieval method are described in Sect.7.

Data sets
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 for March 2011, 2012, 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 SAR processing of CryoSat-2 utilizes an unfocussed aperture synthesis technique which utilizes Doppler beam formation to reduce the footprint size in comparison to a beam-limited altimeter.The effective footprint size after postprocessing is pulse-limited at ∼ 1650 m in the across-track direction and pulse-Doppler-limited to be ∼ 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.The bandwidth-limited range resolution is 3.125 ns (0.469 m range resolution in vacuo); however the range sampling is 1.563 ns (0.234 m), which is done in order to avoid aliasing in the Fourier domain after the power envelope is taken of the signal (Jensen, 1999).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 powerdetected 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 to best simulate the returns from SAR mode.
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 a 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 each individual data point used in this study.Additionally, a bilinear interpolation of the EGM08 geoid (Pavlis et al., 2008) has also been estimated for each measurement and subtracted from the elevation.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.
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).Available data from the final processing of the 2013 campaign suggest an additional uncertainty of 0.016 ± 0.024 cm in the snow depth data and 0.002±0.069cm uncertainty in the laser freeboard from the use of the quick-look data.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, 2013;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 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 an exponential autocorrelation.Where appropriate, we note in the text where assumptions have been made and attempt to justify them.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.7.

Physical model for CryoSat-2 echoes
Here we provide the theoretical basis for modeling the mean echo power 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, ERS, 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 phaseweighted 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 (τ ), 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; 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.This will occur in practice when the snow pack does not contain wet or saline layers which will attenuate the signal, and the surface backscatter coefficient of the sea ice layer is much higher than the sea ice volume backscatter as well as the surface and volume backscatter from the snow layer.The errors in this assumption are explored in Sect.6.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 where σ c = 2σ/c, which is the standard deviation of the surface height 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 of the snow-ice interface will take over the radar footprint.Only limited data are available to determine the spatial statistics of Arctic sea ice.Using laser altimetry data, Rivas et al. (2006) show that a Gaussian height distribution is mainly valid for smooth ice, whereas for rough ice an exponential autocorrelation and Lorentzian power spectral density more accurately characterizes the surface roughness.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, and 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 the polar and azimuth angles, respectively, 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., 2014), which should be taken into account if a more physically exact characterization of the impulse response is desired.
Following Wingham et al. (2004), Eq. ( 4) can be extended for application with CryoSat-2 through the addition of a synthetic beam gain term.The synthetic beam gain used in the processor which constructs the level 1B waveforms is defined as d (cos ζ − sin ξ k ), where d (ζ, ξ k ) is a function defined by the discrete Fourier transform of a Hamming window, which is the window used when stacking the data.d (ζ, ξ k ) is a function of the angle between the direction of a scattering element and the satellite velocity vector, ζ , and the look angle of synthetic 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.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).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. (2004) 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 be from a smoothly undulating surface.In this case, 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 this type of scattering is where h m is the root-mean-square (rms) height deviation, l is its length scale, and k 0 is the carrier wave number.The radar backscatter coefficient can then be expressed following Hagfors (1970) as where R 0 is the Fresnel reflection coefficient for normal incidence.The scattering model uses the Helmholtz-Kirchoff diffraction formula, Hagfors (1970) states that the assumption of this model is that phase modulations are taken to be deep, kσ > 1.Additional assumptions of this type of model are more explicitly stated in Ulaby et al. (1986) to be that (1) the correlation length is larger than the electromagnetic wavelength, l > λ, and that (2) the radius of curvature of the surface is large with respect to the wavelength, l 2 /(2σ √ π/6 > λ.Rivas et al. (2006) calculated typical values for Arctic sea ice of l=0.6 m for young ice and l = 3 m for deformed ice, and σ = 0.03 m for young ice and σ = 0.1 m for deformed ice.Using these results the assumptions of the model are thus generally valid.
For the small incidence angles encountered by CryoSat-2, we assume cosγ ≈ 1.Using this relation the radar backscatter can be expressed following Hagfors (1970) as We relate this to variations in σ 0 (γ , ω) by considering only variations in the polar angle direction and taking φ=γ and cosγ ≈ 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. In choosing Eq. ( 9) 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 as heavily ridged ice, this assumption will likely introduce additional uncertainty in the results since the height distribution will likely not be Gaussian.
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. (2004) and Wingham (1995).The expression used here for the impulse response of the multi-looked echo power shape follows the expression described in Galin et al. (2013) with the additional use of Eq. ( 9) to include the variation of backscatter with incidence angle.Consequently, the full expression for the impulse response is written as Table 1 summarizes the parameters and symbols used in the above equation.The expression after the summation in the second, third, and fourth lines of Eq. ( 10) corresponds to the contribution of the elliptical antenna pattern, which is taken from Wingham and Wallis (2010).The fifth line of Eq. ( 10) corresponds to the variation of backscatter with incidence angle.The sixth 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 (10) 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 surfaces consisting of a sea ice lead and floe within the footprint.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 (10) shows that as α becomes large, it dominates the decay of the trailing edge of the waveform.

Lead returns
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).Figure 3 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 over smooth leads (σ ∼ 0; p (τ ) = δ(τ )) with a suitably large value of α that the shape 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 shape.This is also illustrated in the black color 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 surface returns
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 pulse-limited altimeters, was shown previously by Wingham et al. (2004).One point to note is the 1.563 ns sampling rate 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 leastsquares fitting procedure to estimate the mean scattering point and mean surface roughness within CryoSat-2 echoes from varying surface types.Hereinafter, we refer to this waveform-fitting method as the CS2WfF method.This leastsquares 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 lookup 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 1 to 10 s 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.

Fitting routine
In order to speed up calculation and enable fitting of individual waveforms, we calculate a lookup 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 lookup table, a least-squares fitting routine using a bounded trust region Newton method (MAT-LAB 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) α, and (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 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) through the use of the pulse peakiness and stack standard deviation parameters.First, the pulse peakiness parameter is calculated following Armitage and Davidson (2014) as 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 α 0 100 < α 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 Figure 5. Ratio of tail to peak power, <P_r_l> / max(P_r) where <P_r_l> is the average power of points within 10 ns following the point of peak power.These results are taken directly from the physical model in Sect.3.This is used to provide an initial guess in the fitting of waveforms over leads.The x axis is a logarithmic scale to better show the variation over the large dynamic range of α.
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.In other words, 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 α with a maximum uncertainty of 3 cm, but the uncertainties 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 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.Waveforms that do not meet this requirement are not fitted and no elevation is retrieved; for floe points this occurs in approximately 10 % of the data for the March time periods.In total, approximately 60 % of the CryoSat-2 waveforms are used for elevation retrievals during the March campaigns; this rate is largely determined by the pulse peakiness and stack standard deviation requirements.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 ∼ 80 % of the peak power, the fitted results shown in Fig. 8 demonstrate that 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 1.563 ns sampling rate 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 Fig. 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 interrelated 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 and Laxon (2004) and Laxon (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 electromagnetic 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).In terms of observational data, the 2008 CryoVEx field experiment described by Willatt et al. (2011) showed that when cold, dry snow is present, 80 % of Kuband radar returns were closer to the snow-ice interface than the air-snow interface.However, Willatt et al. (2011) also show that during the CryoVEx 2006 experiment, when warm surface temperatures and complex snow stratigraphy were present, only 25 % of Ku-band radar returns were closer to the snow-ice interface.The assumption of the dominant radar return being from the snow-ice interface needs to be considered on a regional and seasonal basis.Sea ice freeboard is defined here to be the height of the ice layer above sea level 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 as well as the static EGM08 geoid outlined in Sect. 2. We then apply the retracking correction, which is taken from the CS2WfF 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.25 km elevation grids are created, one grid containing the floe elevations and another grid containing the lead elevations.
Freeboard is calculated by subtracting the lead elevation grid from the floe elevation grid (Eq.15).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 125 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 CS2WfF method.A map of the surface roughness (excluding sea ice leads) is shown in Fig. 10.The sea ice floe roughness also corresponds well to what may be expected 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 be added to account for variations of the speed of light within the snow pack on sea ice.This is given as f b = f b radar + h c , where the correction factor, h c , is given as Figure 10.Map of the mean gridded surface roughness, σ , from CryoSat-2 excluding sea ice lead points.Surface roughness from low ice concentration areas near the sea ice edge have also been included; these areas have a high surface roughness due to the presence of ocean waves.
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 f b 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 , and ρ 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 assumptions which are discussed in Kurtz et al. (2013b) to be consistent with the IceBridge 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. ( 18) can be written as

Comparison of CryoSat-2 freeboard and thickness data
In this section we compare CryoSat-2 retrieved freeboard data using the CS2WfF 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 CS2WfF 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. (2012), but was not done in the ELTF freeboard retrievals.Additional differences include (but are not limited to) the exact definition of the first peak, which was not explicitly defined; Laxon et al. (2013) also used a mean sea surface height data set built from a full year of CryoSat-2 observations of ocean elevation, whereas in this study we use the EGM08 geoid and exclude a correction for the dynamic topography.The removal of the time-varying and static sea surface height parameters affect the absolute uncertainties in the retrieved freeboard, but they affect the waveform and ELTF methods equally since they are applied consistently.Therefore, the comparisons done in this study are similar, but not exact, reproductions of methodologies.The purpose of the comparison is to highlight the physical basis between differences in the retracking methods.The h c snow speed-of-light correction was also not applied in the comparison between threshold and waveform retrackers, but this will not affect the comparison because it is equivalent for both data sets.The mean difference in sea ice lead elevations retrieved by the CS2WfF method described in this study and the empirical tracker described in Giles et al. (2007) is 2.8 cm (Giles et al., 2007, tracker -CS2WfF method), and the correlation is 0.7.The most significant difference between the freeboard retrieval method of Laxon et al. (2013) and the CS2WfF 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 CS2WfF method.The mean freeboard differences are 11.9, 12.7, and 11.5 cm for the March 2011-March 2013 periods.This corresponds to mean ice thickness differences of 1.12, 1.19, and 1.08 m using Eq. ( 19).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. ( 19), this corresponds to a theoretical mean sea ice thickness of snowfree ice of 2.9 m, which will be higher once one considers the contribution of snow.The CS2WfF 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-oflight correction was not applied in the comparison between threshold and waveform retrackers, 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 CS2WfF methods, we now compare both methods to independent data collected from three measurement campaigns of the Operation Ice-Bridge 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. ( 19).
Table 2 summarizes the comparison between the Ice-Bridge observations, the CS2WfF method, and the ELTF method; all comparisons are done with the gridded data.For the CS2WfF method, the mean freeboard difference (CS2-IceBridge) ranges from 1 to 3 cm, while the ice thickness difference ranges from 11 to 23 cm.Histograms of the freeboard differences are shown in Fig. 13.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 offnadir 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 to 15.9 cm, which corresponds to ice thickness differences of 112-149 cm; this is significantly higher than the CS2WfF method.For the three IceBridge campaign periods from March 2011 to 2013, mean differences (CS2 -IceBridge) of 0.144 and 1.351 m are found between the freeboard and thickness retrievals, respectively, using a 50 % sea ice floe threshold retracker, while mean differences of 0.019 and 0.182 m are found when using the CS2WfF 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 CS2WfF 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 CS2WfF method gives a mean difference which compares much better to the IceBridge ice freeboard data using explicit geophysical arguments.The correlation coefficients between the IceBridge observations and ELTF method are also lower than those found by Laxon et al. (2013) for the 2011 and 2012 time periods.The reason for the large difference is not clear at this time, but some differences are due to the previously mentioned changes between the ELTF method and Laxon et al. (2013) data set as well as the fact that the correlations in Laxon et al. (2013) are for ice thickness rather than freeboard as was done here.
The rms difference between the IceBridge data and CryoSat-2 freeboard retrievals ranges from 7.4 to 11.1 cm for the CS2WfF 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., 2013b) for the compared grid points is 5.9, 7.6, 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 static and dynamic sea surface height uncertainties, and is generally a function of distance to the nearest lead.Assuming the uncertainties between the IceBridge and CryoSat-2 data sets are uncorrelated, the observed standard deviation of differences between the two measurements, σ diff , is due to the combined uncertainty of the individual components.From stan-dard error propagation, this is written as  uncertainty in the CryoSat-2 sea ice freeboard can then be calculated as 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, 6.6, and 3.8 cm for the respective 2011-2013 campaigns.The difference in uncertainty for the different years could be due to the different proportion of ice types sampled in the three campaigns shown in Fig. 2. The correlation between the CS2WfF method and Ice-Bridge data varies substantially between the campaigns, but the correlation between the CS2WfF 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 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 that the expected correlations for two identical data sets with estimated uncertainties equivalent for the 2011, 2012, and 2013 campaigns are 0.46±0.05,0.55±0.03,and 0.60±0.04,respectively.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 nontemporally 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.It can also mean that the uncertainty in the IceBridge and/or CryoSat-2 freeboard estimates are higher than estimated.Sampling differences between the three campaigns may also be a factor.A more detailed comparison between time coincident IceBridge data flights which underflew CryoSat-2 will be the subject of a future study.

Estimation of errors due to radar penetration
In the retrieval of surface elevation, it was assumed that scattering only occurs from the sea ice surface.This assumption will introduce negligible error only when the energy returned from the sea ice surface is much higher than the energy from the snow surface and the interior of the snow and sea ice layers.Here, we estimate the error that this assumption introduces into the retrieval of sea ice elevation and freeboard 1233 through consideration of the surface and volume scattering contributions of snow and sea ice.
Following Wingham et al. (2004) after inclusion of a depth-dependent backscatter term, the received radar echo can be written as where v (τ ) is the scattering cross section per unit volume as a function of delay time.In the following we assume that, since the incidence angles are small, only the two-way vertical path within the volume is considered; we also neglect the contribution of multiple scattering since this will be small.Here, we follow the approach outlined by Arthern et al. (2001) to define v (τ ) in terms of physical parameters, which satisfies the relation surf-snow +σ 0 vol-snow +σ 0 surf-ice +σ 0 vol-ice = σ 0 , where σ 0 is the total backscatter from the surface.The parameters σ 0 surf-snow and σ 0 surf-ice are the surface backscatter coefficients of snow and ice.σ 0 vol-snow and σ 0 vol-ice are the integrated volume backscatter of the snow and ice.The location of the sea ice layer is defined as τ = 0 since this is the desired location of our tracking point.For sea ice with a snow layer, v (τ ) is written as where Equation ( 22) accounts for the attenuation of the signal in the snow and ice layers as well as the loss of power at the transition between layers with different dielectric properties.The general assumptions which lead to Eq. ( 22) are from geometric considerations which are detailed in Wingham (1995) and Arthern et al. (2001), but have been modified here to model a snow layer on top of sea ice.Here, k e-snow and k e-ice are the two-way extinction coefficients of snow and ice, respectively.We use values of k e-snow = 0.1 m −1 and k e-ice = 5 m −1 , which correspond to penetration depths of 10 m for snow and 0.1 m for sea ice (Ulaby et al., 1986).c ice and c snow are the speed of light in ice and snow, where c ice = c n ice and c snow = c n snow .We take n snow = 1.281, which corresponds to snow with a density of 320 kg m −3 (Tiuri et al., 1984), and n ice = 1.732 (Ulaby et al., 1986).k t-ice and k t-snow are the transmission coefficients between the snowice interface and the air-snow interface, respectively.They are determined from the Fresnel equations for normal incidence using the chosen values of n ice and n snow and taking the index of refraction of air to be 1.Both k t-ice and k t-snow are approximately equal to 0.98.Lastly, σ vol-snow and σ vol-ice are the radar backscatter per unit volume for snow and ice.
In this case study, we take the snow backscattering properties to be constant and vary the ice backscattering properties and snow depth.For the constant parameters, we take σ 0 surf-snow = 5 dB and σ vol-snow = 2 dB, which is in the range of values observed from snow on the Antarctic Ice Sheet (Arthern et al., 2001).We take σ vol-ice = −20 dB; this is estimated from the asymptotic value of laboratory experiments of radar backscatter as a function of viewing angle (Beaven et al., 1995).For the snow depth, we simulate a range of values from 0.02 to 0.5 m.We take the snow surface layer to be flat, and take the surface roughness of the ice layer, σ c , to be half of the snow depth to simulate the increased roughness and snow depth which typically occurs over multiyear ice.Lastly, we use three cases with σ 0 surf-ice =8 dB, 15 dB, and 20 dB.These cases are meant to simulate the effect of radar penetration into the snow when the surface backscattering between the snow and ice layers are comparable; an intermediate case; and when the backscattering from the ice layer is much higher than the snow, as was assumed in the retrieval method.
Figure 14 shows the impact of scattering from the snow layer on the CryoSat-2 waveform when the snow depth is 0.3 m and the backscatter from the sea ice and snow surfaces are comparable at σ 0 surf-ice = 8 dB.Due to the sampling resolution of CryoSat-2, the impact of a strong backscatter from the snow layer is subtle and is seen mainly on the waveform leading edge.The fitted waveform retrieves a surface elevation which is 0.06 m too high; this is due to the increased backscatter from the snow layer which is not taken into account in the CS2WfF retrieval method.Figure 15 shows the errors which are present in the CS2WfF method for a variable snow depth and ice surface backscatter; these tracking errors are directly translatable to the freeboard error which is incurred.The errors are always positive; that is, the retrieved freeboard from the CS2WfF method is biased high due to the presence of scattering from the snow layer and ice volume.For the optimum case when the backscatter coefficient of ice is much higher than that of snow, σ 0 surf-ice = 20 dB, the overall error is low and less than ∼ 0.03 m even for deep snow.This is the case which was assumed in our retrieval method.For the intermediate case when σ 0 surf-ice = 15 dB the error is low at ∼ 0.02 m or less for snow depth values typical of first-year ice (h s < 0.2 m); for the higher snow depth values typical of multiyear ice, the error is on the order of several centimeters.For the case when the backscatter from the snow layer is comparable to the ice layer, the freeboard error increases as a function of snow depth.For snow depth values typical of first-year ice, the error is about 0.04 m or less, whereas for a maximum snow depth of 0.4 m taken from multiyear ice in the March-April time period of the Warren et al. (1999) climatology, the error is 0.08 m.
The errors presented in these case studies are meant to estimate the freeboard errors which could be present in our re-trieval method.The overall error should be less than ∼ 0.04 m for first-year ice and ∼ 0.08 m for multiyear ice, but will be largely dependent on the snow depth and backscattering properties of the sea ice and snow.From Eq. ( 8) this can be seen to be due to the surface roughness and correlation length of the scattering layers.Attenuation of the radar signal in the snow layer was small in these case studies; this is appropriate for dry snow, which is typically encountered in the Arctic.The loss of power at the snow-air and snow-ice interface is also small owing to the high values of k t-ice and k t-snow , so unless the snow layer strongly attenuates the radar signal due to the presence of wet or salty layers, a significant portion of the transmitted energy will always penetrate to the sea ice interface.For sea ice, a higher volume backscatter coefficient would also change the echo shape significantly; the value chosen in this study led to a negligible contribution to the waveform shape.In order to better quantify the error which is present in actual CryoSat-2 data, the backscattering properties of the ice and snow layers need to be better quantified.The use of Eq. ( 22) could also be incorporated into the CS2WfF method to attempt to account for variable scattering from the snow and ice layers, but this is beyond the scope of the present 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, 15.9, and 11.9 cm and mean sea ice thickness differences of 144.2, 149.3, and 111.9 cm when consistent estimates of snow depth and sea ice density are used.Note that these ice thickness differences are much larger than those reported by Laxon et al. (2013), who compared ice thickness values which were retrieved using different parameterizations for snow depth and sea ice density.The mean freeboard differences for the CS2WfF method were 2.2, 2.5, and 1.1 cm, and the mean sea ice thickness differences were 20.6, 23.3, 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 modeling.The difference is also due to variations in surface roughness and the angular dependence of the backscattering coefficient.A mean bias of 1.9 cm was found in the CS2WfF method freeboard retrievals compared to the Ice-Bridge data; this bias is consistent with the estimated range  (2014).A maximum correlation of 0.57 was found between the IceBridge freeboard and thickness data and the CS2WfF 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 a Lorentzian power spectral density more accurately characterizes the surface roughness over rough ice (Rivas et al. 2006).( 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 additional uncertainty if the backscattered energy from the snow layer is comparable to the surface backscatter from the sea ice layer.The overall freeboard bias due to the lack of consideration of backscattering from the snow layer was estimated to be less than 0.04 m for first-year ice, and less than 0.08 m for multiyear ice; the bias encountered in practice will depend on the backscattering characteristics of the snow and sea ice layers.
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 may be used to reduce geoid errors, which are known to be prevalent in the Arctic Ocean (McAdoo et al., 2013); however whether this approach introduces a higher error due to dynamic sea surface height variations remains to be explored.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.

Figure 1 .
Figure 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 .
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, and (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 .

Figure 7 .
Figure7.Ratio of tail to peak power, <P_r_f> / max(P_r) where <P_r_f> is the average power of points located within 90 to 120 ns following the point of peak power.These results are taken directly from the physical model in Sect.3.This is used to provide an initial guess in the fitting of waveforms of sea ice floes.The x axis is a logarithmic scale to better show the variation over the large dynamic range of α.

Figure 8 .
Figure 8. Example CryoSat-2 waveform fits for sea ice floes.The fitted waveform is shown in blue and the CryoSat-2 data are represented by black dots.(a) and (b) are fits of waveforms which demonstrate the good agreement between the observations and the model.(c) and (d) are fits of waveforms containing multiple peaks in the trailing edge due to the presence of strong off-nadir reflections from smooth ice and/or leads.

Figure 12 .
Figure 12.Freeboard difference between the CS2WfF method and the ELTF method.

Figure 13 .
Figure 13.Histograms of the gridded freeboard difference between IceBridge and the CS2WfF method.

Figure 14 .
Figure 14.Simulated waveform for σ 0 surf-ice = 8 dB and a snow depth of 30 cm.The simulated waveform is shown in black and the sub-sampled points at the sampling resolution of CryoSat-2 are represented by the black dots.The blue line shows the best-fit waveform, which assumes backscattered power originates from the sea ice interface only.

Figure 15 .
Figure 15.Plot of the error which is present in the retrieval of freeboard for the CS2WfF method under the assumption that only backscattered power from the sea ice surface is present.The error is a function of snow depth; three cases are shown for varying strengths of the sea ice surface backscatter.
from 16 to 28 March 2011, and 14 March to 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 utilize new processing techniques to minimize freeboard biases

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

Kurtz et al.: CryoSat-2 sea ice freeboard 1235 bias
due to off-nadir ranging of lead points shown by Armitage and Davidson