chapter2.tex

%
% Chapter 2 - Method
% ------------------
% Jul 1, 2015
%

\chapter{Time series of ice-shelf height}

{\sl
\noindent
This chapter, in full, is a reprint of:\\
Developing improved decadal records of Antarctic ice-shelf height change from multiple satellite radar altimeters, F.\,S.~Paolo, H.\,A.~Fricker, L.~Padman, {\rm Remote Sens. Environ.} in revision 
}


\section{Abstract}

Antarctica's ice shelves, the floating extensions of the ice sheet, exert an important dynamic constraint on the flow of ice from the grounded ice sheet to the ocean, and hence on changes in global sea level. Thinning of an ice shelf reduces its ability to restrain the ice discharge from the grounded ice-sheet interior. However, our understanding of how ice-shelf processes couple ice-sheet changes to climate variability is still rudimentary. In part, this is due to the brevity and low resolution of surveys of ice-shelf thickness relative to the time scales on which ice-sheet mass fluctuates. We present improved procedures to construct 18-year (1994--2012) time series of ice-shelf height around the entire Antarctic continent, merging data from multiple overlapping satellite radar altimeter missions (ERS-1, ERS-2, and Envisat). We apply an averaging scheme to enhance the signal-to-noise ratio of height changes over the floating ice shelves, and extract low-order polynomial trends using a robust approach accounting for both bias and variance in the fit (regularized regression with cross-validation). We identify the main processes affecting the estimation of ice-shelf height from satellite radar altimeters. We estimate uncertainties by bootstrap resampling of the residuals of the fit, allowing us to construct formal confidence intervals (unlike the standard error-propagation approach). Our results show that, for the ice shelves, densification of the surface is a more important effect than penetration of the radar signal in biasing the estimated height changes. The 18-year record of surface height allows us to map the temporal progression and spatial structure of changes in ice-shelf height, which provides insights on how ice shelves respond to the changing atmospheric and oceanic conditions.

\section{Introduction}

The Antarctic Ice Sheet contains ice above floatation equivalent to ~58 m of global sea-level rise \parencite{Fretwell2013} and, over centennial-to-millennial time scales, plays an important role in pacing sea-level changes \parencite{Alley2005}. Satellite measurements over the past two decades show mass loss from the grounded ice sheet \parencite{Shepherd2012} including accelerating loss in the Amundsen Sea Embayment \parencite{Sutterley2014}, a region with large potential to contribute to sea level rise \parencite{Rignot2014, Joughin2014}. There is an urgent need to understand the mechanisms behind these current changes as a step toward predicting Antarctica's contribution to global sea-level change over the next century. 

Most ice mass loss from Antarctica takes place through iceberg calving and basal melting from the ice shelves, the floating extensions of the ice sheet \parencite{Depoorter2013, Joughin2012, Rignot2013}. Ice shelves restrain the discharge of grounded ice into the ocean through a buttressing effect \parencite{Joughin2011, Schoof2007}. Small perturbations in the adjacent oceanic and atmospheric conditions can have a large impact on the extent and thickness of ice shelves \parencite{Dutrieux2014, Rignot2004, Scambos2004}, which reduces their buttressing capability. The response of ice shelves to climate change is, therefore, a key component in assessing future loss of grounded ice.

Our understanding of the many processes that affect ice-shelf mass balance is too rudimentary to allow prediction of ice-sheet change under projected future climate states. There are two complementary ways forward: to develop our theoretical framework of the actual mass-loss processes so they can be better represented in models; and to empirically relate observed ice-sheet change to ocean and atmospheric variability. In a recent study \parencite[][see chap. 3]{Paolo2015} we reported changes in Antarctic ice shelf height, and inferred thickness during the 18-year period 1994--2012. The temporal and spatial resolutions of that record were ~3 months and ~30 km, respectively. This record will be used to improve knowledge of ice-shelf response to climate by comparing measured variability of ice thickness with measured or modeled changes in the ocean and atmosphere. This continuous, highly-resolved record overcomes the limitations of previous studies of ice-shelf thickness, which analyzed much shorter records and/or predominantly reported simple linear trends for large areas \parencite{Pritchard2012, Shepherd2010, Zwally2005}.

In this paper we document comprehensive methods for constructing continuous records from multiple satellite radar altimeters to obtain reliable time series of height change over the longest possible time period, providing detailed justification for the data processing approach used by \textcite{Paolo2015}. We present improved procedures for merging data from three overlapping satellite missions, enhancing the signal-to-noise ratio of changes in ice-shelf height, and extracting 18-year mean trends and acceleration. The method reveals complex patterns of ice-shelf height variability in both time and space. We introduce an alternative approach to the standard error propagation for the uncertainty analysis: bootstrapping applied to time-dependent data. 

\section{Satellite radar altimeter missions}
\label{sat-ra}

We used data from three European Space Agency (ESA) satellite radar altimeter (RA) missions: the European Remote Sensing Satellite-1 and Satellite-2 (ERS-1, 1991--1996; and ERS-2, 1995--2003); and the Environmental Satellite (Envisat, 2002--2012). Each satellite carried a standard pulse-limited altimeter with a footprint size of $\sim$3--5 km over the predominant flatter portions of the ice sheets. For reasons explained below in Section \ref{sat-ra}, we do not use the first two years of ERS-1 data, resulting in an 18-year continuous record for 1994--2012.

ERS-1 was launched in July 1991 and operated between December 1991 and June 1996. It flew at an altitude of $\sim$785 km with an inclination of 98.5 degrees (latitudinal limit of 81.5$^\circ$), with three different orbital repeat periods: 3-day (Ice Phase, December 1991 to March 1992 and December 1993 to April 1993), 35-day (Multidisciplinary Phase, April 1992 to December 1993 and March 1995 to June 1996) and 168-day (Geodetic Phase, April 1994 to March 1995). ERS-2 was launched in April 1995 and operated from May 1995 to June 2003 in an orbit with a repeat period of 35 days, following the ERS-1 35-day orbit. In June 2003 the on-board tape recorder used for the RA data failed, and the mission continued without altimetry until July 2011. The RA system on these two satellites was based on a linearly polarized antenna at 13.8 GHz Ku-band, yielding along-track measurements that are $\sim$340 m apart (sampling rate of 20 Hz, where each "sample" is the average of several radar echoes).

The ERS-1/ERS-2 orbit covered about 85\% of the Antarctic ice-shelf area, missing only the southern portions of the large Filchner-Ronne (Fig.~\ref{c2f1}) and Ross ice shelves. Over the ice shelves, ERS-1 and ERS-2 operated in both a standard `ocean mode' and a specialized `ice mode', alternating between them. The range window (the segment of return echo that is recorded) for ice mode was four times wider than for ocean mode, increasing the chances of capturing return signals over rough topographic surfaces. However the broader range window resulted in four times coarser sampling of the return waveform, leading to a less precise estimation of arrival time and, therefore, of the surface height.


\begin{figure}[!ht]
  \centering
  \includegraphics[width=.78\textwidth]{img/map_xovers_grid_v8.png}
  \caption[ERS-1 radar-altimeter data coverage over the Weddell Sea]{
  \ssp \footnotesize
  ERS-1 radar-altimeter data coverage over the Weddell Sea region of Antarctica, for a typical 3-month period with the 35-day repeat orbit. Red dots are crossover locations on the ice shelves. The $\sim$30-km grid that we used to average the crossovers is overlaid. 
  ice shelves.} 
  \label{c2f1}
\end{figure}


Envisat was launched in March 2002 into the same 35-day orbit as ERS-1 and ERS-2, providing data at the same 20 Hz sampling rate. The altimeter on this satellite (RA-2) \parencite{Roca2009}, a linearly-polarized radar, operated via a single antenna dish on the Ku-band (13.6 GHz) and S-band (3.2 GHz). This dual frequency enables the correction of height-measurement errors introduced by the ionosphere. RA-2 had an improved precision over ice surfaces compared to the ERS altimeters, and operated with three sampling modes: `fine', `medium' and `coarse', which differ in their sampling of the waveform of the return echo. Envisat stopped operating in April 2012.

We obtained Level-2 RA data as Ice Data Records (IDR) for each mission from NASA's Goddard Space Flight Center (GSFC) Ice Altimetry group (\url{http://icesat4.gsfc.nasa.gov/}). At GSFC, the RA waveforms were retracked using a multi-parameter range-retracking algorithm (the $\beta$-retracker) \parencite{Zwally2001}. The following corrections were applied by GSFC: atmospheric range corrections; instrument corrections; slope corrections; ocean and solid earth tides \parencite{Brenner1983, Zwally2001, Zwally2005}; (for ERS) removal of a 40.9 cm bias from ERS-1 heights to account for a change in instrument parameter used for ERS-2 \parencite{Femenias1996}; corrections for drifts in the ultra-stable oscillator and bias changes in the scanning point target response that are obtained from ESA; and upgraded orbits (DGM-E04 orbits for ERS) which have a radial orbit precision of 5-6 cm \parencite{Scharroo1998}.


\section{Processing satellite radar altimeter data}

Our determination of height changes over the ice shelves from multi-mission satellite RA data is based on crossover analysis \parencite[e.g.,][]{Davis2004, Wingham2009, Zwally2005}, which estimates change in surface height at intersections between time-separated ascending and descending satellite tracks. Along-track (or repeat-cycle) analysis methods that provide higher point density have recently been introduced \parencite{Flament2012, Moholdt2010, Pritchard2012}; however, crossover analysis remains the most precise technique, since differences are derived from precisely co-located height estimates on the ascending and descending tracks. For crossover analysis interpolation is performed between measurements $\sim$340 m apart in the along-track direction (where there is a well defined profile from a consistent set of measurements, i.e., one orbit); in contrast for repeat-cycle analysis each measurement is projected onto a reference ground track, with the actual orbits varying up to $\sim$2 km from each other from repeat to repeat (so interpolation is performed in the across-track direction without having a good estimate of the slope). The spacing of crossovers becomes much smaller as latitude approaches the satellite's turning point (Fig.~\ref{c2f1}).

In this section we describe the steps required to process RA data from multiple satellite missions. Steps include (see flowchart in Fig.~\ref{c2f2}): subsetting data over floating ice; data editing and additional corrections required for floating ice and for radar signal interactions with the surface layer on the ice shelves; crossover analysis; methods for merging data records from the different missions; and the averaging scheme required to achieve satisfactory data coverage and accuracy. The output is 18-year long records of ice-shelf surface height, which we then analyzed using statistical techniques to derive long-term trends, acceleration and uncertainties.


\begin{figure}[!ht]
  \centering
  \includegraphics[width=.6\textwidth]{img/flowchart_v4.pdf}
  \caption[Flowchart of our RA data processing scheme]{
  \ssp \footnotesize
  Flowchart of our RA data processing scheme, showing the steps performed to construct an 18-year record of ice-shelf surface-height change from multiple satellite altimeters. Each step is described in Sections \ref{sec:ra-processing} and \ref{sec:ra-analysis}.}
  \label{c2f2}
\end{figure}


\subsection{Initial RA data selection and editing}
\label{sec:ra-processing}

\subsubsection{Antarctic ice shelf mask}

For this study we are interested only in Antarctica's floating ice shelves; therefore, we subsetted all of our RA datasets based on a 1-km resolution Antarctic ice-sheet mask \parencite{Depoorter2013} that was constructed using a composite of InSAR \parencite{Rignot2011}, ICESat \parencite{Brunt2010, Fricker2009}, MODIS MOA \parencite{Scambos2007} and ASAID \parencite{Bindschadler2011} products. This high-resolution product allows us to exclude all grounded ice, which includes islands, ice rises and rumples that are completely surrounded by ice shelves.

We considered only the portions of ice shelves that were present and afloat through the entire 18-year observation period. We excluded regions where ice front advance or retreat occurred during the 18 years and portions of a few ice shelves that became ungrounded between 1994 and 2012.

\subsubsection{Altimeter data editing}

We excluded data if any of the following conditions were met: the return waveform had no leading edge, or had specular shape indicating surface melt ponds and melt streams \parencite{Phillips1998} as indicated by flags in the IDR; any geophysical correction was missing from the IDR; the crossover point was less than 3 km from any ice-shelf boundary (grounding line and ice front). The last criterion was designed to avoid biasing estimates (i) with high topographic gradients within the pulse-limited radar footprint, (ii) with changes in grounded ice that are an order of magnitude higher than over floating ice, and (iii) where geophysical corrections for ocean processes such as tides are unreliable, especially close to the grounding line where the ice is not in full hydrostatic equilibrium \parencite{Fricker2006}.


\subsubsection{Additional data corrections}

The ocean tide corrections applied to the IDR (from the CSR3 global tide model) are not accurate in Antarctica \parencite{King2005}; for that reason, we removed this correction from the data and applied an improved tide correction using the regional Circum-Antarctic Tidal Simulation model (CATS2008a; an updated version of the inverse tide model described by \textcite{Padman2002}). This tide model has higher resolution ($\sim$4 km) and a more accurate land mask than global models, resulting in more accurate tide prediction close to the coast. We also corrected for the ocean tidal loading (the elastic deformation of the seabed in response to the tide load) using the TPXO7.2 model \parencite{Egbert2002}.

\subsubsection{Estimating height differences}

To estimate the height differences at crossovers we first precisely locate the intersection points between each ascending and descending satellite tracks. We then find the two data points on either side of this intersection for each track pair, interpolate these to estimate the height at actual crossover location, and difference values from both tracks. When some data are missing along the satellite track, we only estimated crossover height-differences if the gap between data points on either side of the intersection point was less than 3 km. This criterion ensured that the RA pulse-limited footprints at both ends of the gap still overlapped enough so that interpolated heights along both tracks at the precise crossover location were consistent representations of the same location on the ice surface.

Our approach to developing time series of ice-shelf height with respect to an initial epoch $t_0$, $h(t-t_0)$, uses a multi-reference time scheme \parencite{Khvorostovsky2012, Li2006} applied to each satellite mission separately. This method computes crossover height differences for each pair of ascending and descending tracks, $\Delta h(t_i,t_j)$, throughout the entire record for a specific mission regardless of time separation between track pairs. For example, in the 9-year Envisat record, time separation of differences, $\Delta t=t_j-t_i$, at crossovers range from close to zero to nine years. 

For ERS-1, where data were acquired in different modes (`ice' and `ocean') at different times over the ice shelves, we analyzed the crossovers $\Delta h(t_i,t_j)$ estimated from pairs of tracks with the same operation mode (ice-mode crossed with ice-mode or ocean-mode crossed with ocean-mode). Results for both modes were practically identical. In contrast, estimates of $\Delta h(t_i,t_j)$ obtained from a mixed-mode crossover (one track in ice-mode, the other in ocean-mode), were not consistent with the same-mode pairings. We also found that ice-mode data by themselves ware insufficient to provide near-full coverage over the ice shelves for the first few years of ERS-1 operation. Hence, for ERS-1 we retained both types of same-mode crossovers (ice-ice and ocean-ocean). For Envisat, we used the fine-mode data only. 

\subsection{Averaging to enhance the signal-to-noise ratio}
\label{sec:ra-analysis}

For a given change in ice thickness, the estimated height-change signal over a floating ice shelf is about an order of magnitude smaller than the same signal over grounded ice. This is because, outside of a narrow ice-flexure zone roughly 1--10 km wide seaward of the grounding line, ice shelves are in hydrostatic equilibrium, with a density that is about 10\% less than seawater. Since ice thickness and mass are the primary variables we wish to document, this reduced signal implies that there is a need to reduce the noise in the height estimates over ice shelves, a problem that is compounded when the goal is to identify variability over short time scales.

To enhance the signal-to-noise ratio over floating ice we implemented two averaging steps to reduce the variance in the height estimates, as described below.

\subsubsection{Spatial gridding and temporal binning}

We defined blocks of height data, $h(\mathbf{x},t_i)$, consisting of all ascending and descending tracks within non-overlapping 3-month time bins and spatial grid cells of 0.75$^\circ \times 0.25^\circ$ in longitude and latitude, respectively ($\sim30 \times 30$ km at 71$^\circ$S). At shorter time intervals, e.g., one month, the data records are often discontinuous. Aggregating data over 3-month bins provides continuous records while still resolving the annual cycle. The choice of grid-cell size ($\sim$30 km) was a compromise between characteristic spatial scales of ice-shelf processes and the spatial distribution of crossovers (Fig.~\ref{c2f1}). We identify each 3-dimensional block (one time bin and spatial cell), by its central time, $t_i \rightarrow i=0,1,2,..,N$, where $N$ is the total number of time bins within each satellite record, and by its central location, $\mathbf x \rightarrow 0^\circ <$ longitude $< 360^\circ$ and $-81.5^\circ <$ latitude $< -64.5^\circ$ (the center of the grid cells).

For each pair of blocks for the same cell but for separate time bins $t_i$ and $t_j$, we developed the set of all crossover height differences, ${\Delta h(t_i,t_j)}$, found by differencing all height values from the two data blocks. During this process we excluded crossover values higher than 15 m (gross outliers). For all pairs of ${t_i,t_j}$, we therefore obtained average height changes at each location as:

\begin{equation}
  \overbar{\Delta h}(\mathbf{x}, t_i, t_j) = \frac{1}{n_\text{ad} + n_\text{da}} 
  \left\{
  n_\text{ad} \, \text{Md}\!\left[ \Delta h_\text{ad}(t_i, t_j) \right] +
  n_\text{da} \, \text{Md}\!\left[ \Delta h_\text{da}(t_i, t_j) \right]
  \right\}
  \label{c2e1}
\end{equation}

\noindent
where $\overbar{\Delta h}$ is the weighted-mean height-change estimate between times $t_i$ and $t_j$, $\Delta h_\text{ad}$ and $\Delta h_\text{da}$ are height changes formed by differencing ascending-descending and descending-ascending satellite ground tracks, respectively, $n_ad$ and $n_da$ are the number of crossovers of each type per block, and $\text{Md}[\,]$ is the median operator. We used the median instead of the mean because individual values of $\Delta h(t_i,t_j)$ can have large errors and their distribution can be non-Gaussian. It is necessary to use the average between ascending/descending and descending/ascending crossovers to remove any time-invariant biases, such as those introduced by the combination of satellite orientation, antenna polarization and ice surface anisotropy.

If we assume that true temporal variations in $h$ are small for small time separations (within the time bin itself, i.e., less than three months), then the statistics of $\Delta h(t_i,t_j)$ are a measure of noise in the height difference estimates. This noise value may be spatially dependent; for example, relatively smooth and flat regions of ice shelf may have smaller noise in $\Delta h$ than regions where the ice-shelf surface is sloping or rough.

\subsubsection{Derivation of height-change time series}

The simplest way to obtain a time series of height change for a specific spatial cell is to difference the first block of data consecutively with all subsequent blocks \parencite[e.g.,][]{Zwally1989}. This approach, however, results in all height differences in the resulting record being dependent upon the quality and data availability of the reference epoch (the first block). Instead, for each cell we evaluated all possible combinations between time bins within each satellite mission \parencite{Zwally2001, Li2006}, and further combinations between the formed differences as in \textcite{Khvorostovsky2012}. This procedure was carried out as follows.

For each grid cell we constructed a set of $N-1$ time series of average height differences with respect to each of the $N-1$ separate epochs, and then formed a full matrix (one per cell) containing one height-change record per row. For example, taking $t_0$ to be the reference epoch (any time could be chosen) we then have:

\begin{equation}
  [\Delta h] = 
  \begin{pmatrix}
    h_{0,1} & h_{0,2} & h_{0,3} & \cdots & h_{0,N} \\
    -- & h_{0,1} + h_{1,2} & h_{0,1} + h_{1,3} & \cdots & h_{0,1} + h_{1,N} \\
    h_{0,2} - h_{1,2} & -- & h_{0,2} + h_{2,3} & \cdots & h_{0,2} + h_{2,N} \\
    h_{0,3} - h_{1,3} & h_{0,3} - h_{2,3} & -- & \cdots & h_{0,3} + h_{3,N} \\
    \cdots & \cdots & \cdots & \cdots & \cdots \\
    h_{0,N} - h_{1,N} & h_{0,N} - h_{2,N} & h_{0,N} - h_{3,N} & \cdots & --
  \end{pmatrix}
  \label{c2e2}
\end{equation}

\noindent
where $h_{i,j} = \overbar{\Delta h}(\mathbf{x}, t_i, t_j)$ from Eq.~\ref{c2e1}. In the example above, every row constitutes a time series of height change with respect to $t_0$. This (final) matrix is constructed by (i) forming the upper triangle of the matrix from the calculated average differences (this is a one-sided matrix containing the height differences between all possible time combinations) (Fig.~\ref{c2f3}, top row), (ii) forming the lower triangle as the negative of the upper one (this is the transposed of the previous one-sided matrix multiplied by $-1$), and (iii) adding the two single-sided matrices to form a full matrix (Fig.~\ref{c2f3}, middle row) and then adding each element $j$ of the reference row (first one in this example, time series with respect to $t_0$) to the respective row $i = j$ (Fig.~\ref{c2f3}, bottom row). We note that any row in the matrix could be used as the reference series (i.e., any epoch); the optimal one, however, will cover the largest time span in the presence of data gaps.


\begin{figure}[!ht]
  \centering
  \includegraphics[width=.85\textwidth]{img/multiref_matrix_v4.png}
  \caption[Representation of the multi-reference time-series approach]{
  \ssp \footnotesize
  Representation of the multi-reference time-series approach. (left) Individual time series of cumulative change. (right) Diagram representing the matrix formed with the time series on the left (one time series per row). From top to bottom is depicted the process of forming single-grid-cell average time series: (top) the one-sided matrix of average differences; (middle) the two-sided matrix; (bottom) referencing all the rows (time series) to a common epoch for posterior averaging (this is the matrix in Eq.~\ref{c2e2}).} 
  \label{c2f3}
\end{figure}


Next, we weight-averaged the matrix rows by the number of crossovers used to form each average difference to produce a mean time series per cell, with reduced statistical error:

\begin{equation}
  h(\mathbf x, \tau) = \langle [\Delta h] \rangle 
  \label{c2e3}
\end{equation}

\noindent
where $\tau = t - t_0$. We propagated individual errors in the averaging procedure assuming normality \parencite{Li2006}.

To evaluate possible signal aliasing we used a small region to test aggregating data in 3-month bins that were overlapped by $1/3$ as in \textcite{Khvorostovsky2012}. The overlapping approach required a much longer computation time, and the resulting average time series showed roughly the same pattern of annual and inter-annual fluctuation as the ones derived from the non-overlapping approach, and so we did not apply this to the full dataset.

\subsection{Surface and volume scattering variation}
\label{sec:bs-variation}

Interpreting heights derived from RA data over snow-covered ice surfaces requires an understanding of the complex interaction of the radar wave with the snowpack and firn (compressed snow that has not yet formed into solid ice). Over ice sheets, part of the incident radar energy penetrates into the snowpack/firn \parencite{Ridley1988} by an amount that depends on the properties of the surface layers including temperature, density, grain-size and moisture content \parencite{Davis1993}, leading to volume scattering with reflections within sub-surface layers and ice lenses. Volume scattering increases the path length of the radiation and, therefore, its travel-time back to the satellite and the inferred height. The waveform shape in the presence of volume scattering is also different than if only surface scattering occurred. In some cases, most of the echo can be from a well-defined sub-surface horizon within the snowpack \parencite{Thomas2008}, leading to a surface-like waveform shape, but a height estimate that is still lower than the true surface.

The penetration depth (where the radar pulse decays by $1/e$ of its initial intensity) varies with density, grain-size, and water content of the snow/firn. Penetration depth can be several meters over the colder and dryer parts of the Antarctic plateau where grain sizes are smaller, but is much lower (on the order of centimeters) over the warmer and wetter ice shelves where grain sizes are larger \parencite{Davis1996}. Note that the dependence is on the bulk electromagnetic properties (and not grain size as such). To demonstrate this we used the \parencite{Davis1993} surface/volume algorithm to estimate the extinction coefficient ($k_\text{e}$; which is inversely proportional to penetration depth) for the entire Antarctic area under the ERS-1 satellite's coverage, using ERS-1 altimeter level-1 ocean-mode waveforms (Phase C) (Fig.~\ref{c2f4}). In general, $k_\text{e}$ is on the order of 4--5 times higher over the ice shelves than elsewhere on the ice sheet. This means that the penetration bias is not as large over the ice shelves as it is over the ice sheet interior.


\begin{figure}[!ht]
  \centering
  \includegraphics[width=.85\textwidth]{img/extinction_coef_v3d.png}
  \caption[Estimated radar altimeter extinction coefficient]{
  \ssp \footnotesize
  Estimated radar altimeter extinction coefficient ($k_\text{e}$) over the Antarctic ice sheet. $k_\text{e}$ was derived from ERS-1 Phase C ocean-mode return waveforms over the ice sheet and ice shelves using a surface/volume scattering algorithm \parencite{Davis1993}. Note that the $k_\text{e}$  (which is inversely related to the penetration depth) is generally higher in the ice shelves than on the plateau. This difference can mostly be explained in terms of grain sizes: on the ice shelf, individual particles are larger than those on the plateau, due to much warmer temperatures, combined with successive melt-freeze cycles (i.e., wet snow metamorphosis) \parencite{Zwally1994}; therefore $k_\text{e}$ is high. On the plateau, grain sizes depend on surface temperature and accumulation rate. In general, as elevation and distance from the coast (continentally) increases, temperatures decrease; therefore the grain sizes and $k_\text{e}$ decrease.
  }
  \label{c2f4}
\end{figure}


The return radar waveform used to estimate the surface height is the sum of surface echo (energy scattered by the surface) and volume echo (energy scattered within the snowpack). Perturbations in the return radar waveform occur due to changes in the properties of the snowpack, which alter the effective backscatter layer, potentially leading to a biased surface-height estimate. Two competing effects (with opposite signs) play a significant role in altering the shape of the return waveform: i) penetration of the signal into the firn layer increases volume scattering (increasing the waveform's leading-edge width); and ii) densification of the surface increases the air-snow dielectric contrast (decreasing the leading-edge width). We investigated corrections to minimize these effects as described in Section \ref{sec:bs-corr} below.


\subsubsection{Correction for changes in backscatter}
\label{sec:bs-corr}

Although much progress has been made in understanding the relation between backscatter fluctuation (resulting mainly from variations in near-surface properties) and altimeter-derived height \parencite{Arthern2001, Davis1993, Legresy1998, Partington1989, Remy2012, Ridley1988}, it is still an active area of investigation, and its full impact on the RA measurement error remains poorly known \parencite{Remy2012}. To attempt to minimize the impact of fluctuations in backscatter on the estimated height, an empirical adjustment using backscatter information is performed \parencite[as have been applied by,][]{Davis2004, Khvorostovsky2012, Remy2009, Wingham2006, Wingham1998, Zwally2005}.

We used the altimeter automatic gain control (AGC; proportional to the log of the gain) as a measure of backscattered power to adjust our height-change time series using corresponding series of changes in AGC, $g(\mathbf x,\tau)$ \parencite{Zwally2005}. We formed these backscatter-change time series for each grid cell in the same way as the height-change series described above. We derived a spatially varying sensitivity factor to scale the backscatter values at each grid-cell:

\begin{equation}
  S(\mathbf x) = \frac{\Delta h(\mathbf x)}{\Delta g(\mathbf x)}
\end{equation}

\noindent
where $h$ and $g$ are all the height and backscatter differences formed at location $\mathbf x$. In other words, this sensitivity is the regression slope of the correlation between height and backscatter changes that we derived independently at each grid cell and for each RA mission. We then corrected each time series as:

\begin{equation}
  h(\mathbf x,\tau)_\text{cor} = h(\mathbf x,\tau) 
    - S(\mathbf x) \, g(\mathbf x,\tau) - h_0(\mathbf x)
  \label{c2e4}
\end{equation}

\noindent
where $h_\text{cor}$ is backscatter-corrected height-change, $S$ is the sensitivity factor (altimeter dependent), $g$ is proportional to gain change, and $h_0$ is the regression intercept (not needed for correcting full independent records relative to some epoch). For the regression procedure we used a robust method, `maximum likelihood-type estimator' \parencite[][p. 43]{Huber1981}, instead of the commonly used least-squares approach, which can be sensitive to outliers (especially when few data are available).

In selecting an optimal backscatter correction scheme we tested three different approaches. We correlated series of backscatter and height changes: i) using the cross-calibrated full 18-year records, ii) using data within 2 and 3-year sliding windows (similar to \textcite{Khvorostovsky2012}) and iii) using independent single-mission records. For each approach we performed the correlation using both the original time series (absolute values) and the differenced series (derivative). We found that the third approach was the only one to perform consistently under a variety of conditions (e.g., high vs. low variance, high vs. low correlation, high vs. low heteroscedasticity). Note that using differenced series is equivalent to using high-pass filtered records and, therefore, higher correlations are expected since backscatter change is primarily a function of seasonality. However, in correlating differenced series only the short-term fluctuations in backscatter are taken into account. Hence, we corrected our height-change records by correlating series of absolute values using method (iii).

\subsubsection{Densification from changes in the firn column}

One approach to account for changing surface mass balance over an ice shelf, specifically firn compaction/densification, is to use an atmospheric-based model to predict how the firn column evolves under modeled climate conditions. For example, a laser altimeter study by \textcite{Pritchard2012} used modeling of changes in the firn column to isolate the part of the signal due to fluctuations in the snowpack. We performed a correlation analysis between our backscatter-uncorrected time series, which contain the seasonal response of the ice-shelf surface as detected by the altimeter's automatic gain control (backscatter), and the height changes derived from a firn densification model (FDM). To construct the firn height-change record we used the \textcite{Ligtenberg2011}'s FDM from which we (i) extracted firn height-change series for each location in our grid; (ii) averaged the 48-hour firn records with a 3-month moving window; and (iii) retained the discrete time-values that corresponded to our 3-month ice-shelf height records.

We found no correlation between the firn-height changes from the model and the (backscatter-uncorrected) RA-derived changes in ice-shelf height at the fine resolution required for this study (individual grid cells of $\sim$30 km). Since backscatter changes are indicative (among other processes) of the seasonal fluctuation in the firn height (see text below), some degree of correlation between the uncorrected RA records (which contain at least part of the seasonal firn-height cycle) and the FDM is expected. The firn model was able to provide, however, an estimate of the expected seasonal variation of the firn height at the larger scale (Fig.~\ref{c2f5}). Additionally, since some degree of penetration through fresh snow is expected at the frequency of the Ku-band altimeter, satellite RAs are less sensitive than laser altimeters to fluctuations in surface mass balance. For these reasons we found no justification to use an FDM to make any firn-change correction to RA-derived height when converting to thickness and mass.


\begin{figure}[!ht]
  \centering
  \vspace{.4cm}
  \includegraphics[width=.85\textwidth]{img/correlation_h_b_f_v7.png}
  \caption[Correlation between backscatter (AGC), ice-shelf height]{
  \ssp \footnotesize
  Correlation between backscatter (AGC), ice-shelf height and firn-height change derived from a firn densification model. (left column) Derivative of circum-Antarctic-wide average time series; (right column) respective correlations. From top to bottom, correlations between the derivative of backscatter and derivatives of: (a) uncorrected ice-shelf height; (b) corrected ice-shelf height (at the grid-cell level); and (c) firn-column height time series.
  } 
  \label{c2f5}
\end{figure}

\subsection{Merging individual satellite mission time series}

The analyses described in Sections \ref{sec:ra-processing} to \ref{sec:bs-variation} led to independent time series for each of the three satellite missions. These time series then needed to be concatenated and further quality-controlled to generate the long multi-mission record of ice-shelf height. At each grid cell we cross-calibrated the records from the three satellite missions by computing the offsets (weighted-average of differences) between the time series during the periods when consecutive missions overlapped (ERS-1 to ERS-2: 1995--1996 and ERS-2 to Envisat: 2002--2003). We excluded the first two years of ERS-1 data (1992--1994) from our analysis because they show an apparent strong negative anomaly that extends over several sectors of Antarctica, and this anomaly also appears in the altimeter's backscatter. In the absence of evidence supporting such a strong continental-scale signal, we assume this feature is an artifact of the data. Our final dataset therefore spans 18 years from 1994 to 2012.

We corrected a few step changes greater than 3 m (usually associated with anomalous backscatter) by leveling the segments on each side by the difference in their medians, and removed any peak greater than 3 standard deviations from the polynomial trend (see below). We rejected full 18-year records if time series from any of the three missions did not overlap due to data gaps, or if the full data set spanned less than 70\% of the 18-year time interval.

\section{Time series analysis}

Estimating the evolution of changes in ice-shelf height is an important step for fingerprinting the oceanic and atmospheric processes driving such changes. A challenging task in analysis of any climate data is the extraction of underlying trends from relatively short and noisy records, particularly when there is considerable change in variance within the record itself (heterocesdasticity). Ice shelves couple with the atmosphere (temperature, radiation balance, precipitation and wind stress), ocean (circulation, temperatures and waves), and sea ice (concentration and thickness), over a wide range of time scales \parencite{Paolo2015}. Estimated trends are, therefore, strongly dependent upon length of the time interval over which the trends are calculated.

Previous authors have adopted a variety of approaches for trend extraction such as fitting simple straight lines \parencite{Pritchard2012, Shepherd2010}, quadratics \parencite{Wingham2009}, cubics \parencite{Schenk2012}, sinusoids \parencite{Zwally2005}, and using autoregressive models \parencite{Davis2005}. We took a different approach by (i) selecting the best fit from a range of possible models (polynomial of degree 0, 1, 2 or 3), and (ii) taking into account both bias and variance in the fitting procedure (i.e., regularization).

\subsection{Extracting trends from height-change time series}

We fitted low-order polynomial models to each time series (Fig.~\ref{c2f6}), allowing the degree $m$ to vary from 0 (constant height) to 3 (cubic fit), using the lasso approach for regularized regression \parencite{Tibshirani1996}. The regularization parameter (and degree of the polynomial) was selected by ten-fold cross-validation, where the optimal fit is the one that minimizes the average `mean square prediction error' across folds (i.e., on each fold a model is fitted to 90\% of the data, testing a sequence of 100 regularization parameters; the mean square error (MSE) of each fitted model/parameter is then computed with respect to the 10\% of the data left out) \parencite{Friedman2010}. With this bias-variance-tradeoff approach, our objective is to allow a degree of curvature for estimating potential acceleration of height-change rate while avoiding the rapid oscillations introduced by ordinary least-squares polynomial fits \parencite{Paolo2015}. The polynomials are defined as (dropping the spatial variable for generalization):

\begin{equation}
  \hat h(\tau) = \sum_{m=0}^3 \beta_m \, \tau^m
  \label{c2e5}
\end{equation}

\noindent
where $\hat h$ is the polynomial model, and $\beta_m$ are the coefficients of the polynomial. We then calculated average and instantaneous rates of change as the slope of the secant and tangent lines, respectively, to the fitted polynomial by finite difference approximation. 

\begin{figure}[!ht]
  \centering
  \includegraphics[width=.566\textwidth]{img/map_single_series_v9b.png}
  \caption[Examples of single-grid-cell 18-year time series]{
  Examples of single-grid-cell 18-year time series over the Antarctic ice shelves. (top) Map showing individual grid cells with estimated 18-year mean polynomial trends over the ice shelves for one sector of Antarctica (inset map). The grid was smoothed and interpolated using a Gaussian Kernel with sigma equal to grid-cell size.  (lower panels) Comparison of polynomial versus straight-line fits for trend determination on three individual height-change time series (locations shown in map). The average polynomial rate is estimated from the slope of the line connecting the polynomial end points (Eq.~\ref{c2e6}).}
  \label{c2f6}
\end{figure}


In estimating ice surface-height trends we allow for a trend $\dot p$ in sea-surface height due to the inverse barometer effect associated with a regional atmospheric-pressure trend \parencite{Padman2003} (source: ERA-Interim \parencite{Dee2011}), and a regional sea-level trend $\dot \eta$ (source: AVISO \parencite{LeTraon1998}). The average rate of change $\Delta \hat h/\Delta \tau$ between two times $\tau_i$ and $\tau_j$ is then given by:

\begin{equation}
  \frac{\Delta \hat h}{\Delta \tau} = 
    \frac{\hat h(\tau_j) - \hat h(\tau_i)}{\tau_j - \tau_i}
    + \dot p + \dot \eta
  \label{c2e6}
\end{equation}


To estimate average acceleration we calculated the mean slope of the (numerical) derivative of the polynomial fit (second derivative). We performed the regularized regression procedure in two ways: a) on each individual grid-cell time series (over 2650 cells), which we used to construct spatial maps of average ice-shelf height changes (Fig.~\ref{c2f6}); and b) on ice-shelf-wide and region-wide average time series for estimating regional trends (Fig.~\ref{c2f7}; \parencite{Paolo2015}). 


\begin{figure}[!ht]
  \centering
  \includegraphics[width=.71\textwidth]{img/tseries_v4.png}
  \caption[Regional time series of average Antarctic ice-shelf height]{
  \ssp \footnotesize
  Regional time series of average Antarctic ice-shelf height for: (top) West Antarctic ice shelves; (middle) East Antarctic ice shelves; and (bottom) all Antarctic ice shelves. Adapted from \textcite{Paolo2015}.
  } 
  \label{c2f7}
\end{figure}


Performing regularized regression on regional-average time series, where each point is the average of many other points, is preferred over aggregating estimates of trends from fitted polynomials on individual grid cells. In the latter the regularized regression is not as well constrained due to a lower signal-to-noise ratio (Fig.~\ref{c2f8}). A consequence of this analysis approach to trend estimation is that the sum of ice-shelf/regional values does not equal the total value computed for all of Antarctica. This is because different region sizes lead to different constraints in the regression (i.e., the all-Antarctica trend comes from its own lasso polynomial estimate while each individual ice-shelf trend also comes from its own lasso fit); all estimates are consistent, however, within the formal errors.


\begin{figure}[!ht]
  \centering
  \vspace{.5cm}
  \includegraphics[width=.73\textwidth]{img/map_noisy_series_v3.png}
  \caption[Comparison of noise level between single-grid-cell]{
  \ssp \footnotesize
  Comparison of noise level between single-grid-cell and ice-shelf-average time series. (top) Example of time series derived for three individual grid cells on Larsen-C Ice Shelf. Green lines show polynomial fits. Differences between single-cell time series represent both true spatial variations (signal) and contributions from noise. (bottom) Time series generated by averaging all individual-grid-cell time series over the entire Larsen-C Ice Shelf. The polynomial is fitted to the averaged time series (it is not the average of single-cell fits).
  } 
  \label{c2f8}
\end{figure}


\subsection{Estimating uncertainties}

\subsubsection{Time series error bars}

Our averaging procedure facilitates the formal derivation of statistical error for individual height values within each mean time series. Two factors contribute to lower the uncertainty in the height estimates: i) the large number of crossovers that contribute to the final height estimate per location per time step ($\sim$10--100); and ii) the long-period records (18 years), which makes the derivation of long-term slopes robust to high-frequency fluctuations (e.g., seasonality).

Our uncertainties are meant to reflect: (a) the sampling error (number of crossovers); and (b) the backscatter-related error. Hence, we derived the time-series error bars as:

\begin{equation}
  \delta h(\mathbf x,\tau) = \sqrt{
    \text{SE}[\Delta h(\mathbf x,\tau)]^2
    + \text{SE}[h(\mathbf x)_\text{cor} - h(\mathbf x)_\text{unc}]^2
    }
  \label{c2e7}
\end{equation}

\noindent
where $\text{SE}[\,]$ is the standard error, $\Delta h(\mathbf x,\tau)$ represents all height differences at time $\tau$ and grid-cell location $\mathbf x$ (these are the crossovers being averaged at each block in Eq.~\ref{c2e2}), and $h(\mathbf x)$ is the formed height record at each grid cell, backscatter corrected (cor) and uncorrected (unc). The first term on the right-hand side represents the noise at the grid-cell level (time and spatially dependent). The second term assesses the magnitude of the backscatter-series to height-series correlation (spatially dependent) by looking at the difference between the corrected and uncorrected records. Both terms refer to the mean time series (Eq.~\ref{c2e3}~and~\ref{c2e4}).

\subsubsection{Trend and acceleration error bars}

We estimated the uncertainties in the derived trends from the polynomial fits based on the variance of the residuals. In general this provides a higher uncertainty than that derived from the individual data points. This is because the spread of the residuals arises from true natural variability in addition to observational uncertainty. Furthermore, the spread of the residuals contain valuable information on whether the polynomial fit is a good statistical model for the data. We computed the residuals as:

\begin{equation}
  \varepsilon(\tau) = h(\tau) - \hat h(\tau)
  \label{c2e8}
\end{equation}

\noindent
where $\varepsilon$ are the residuals, $h$ is the observed time series of height changes and $\hat h$ is the polynomial fit. We then used the `bootstrap' approach \parencite{Efron1993} in which the residuals of the polynomial fit are randomly resampled with replacement. That is, every subsample of the residuals contains data from the original sample that can appear multiple times, so that each subsample have the same number of elements as the original data set, but with less information about the noise (with respect to the fit). We then added back each resampled residual to the original fitted model to construct the bootstrap samples:

\begin{equation}
  h^*(\tau) = \hat h(\tau) + \varepsilon^*(\tau)
  \label{c2e9}
\end{equation}

\noindent
where $h^*$ is the bootstrap time series and $\varepsilon^*$ are the resampled residuals. We refit the polynomial model (as previously described) to each bootstrap time series and calculated the following: average rate of change, the derivative of the polynomial and the average rate of change of the derivative (average acceleration). By doing this repeatedly, we constructed an empirical distribution for each parameter of interest. We then estimated formal confidence intervals (95\%) from this distribution. We constructed 1000 bootstrap samples for each individual ice-shelf/region time series, and 500 for each grid-cell time series (a total of $\sim$1,330,000 sets of calculations).

The total uncertainty in the estimated ice-shelf height rate also includes the systematic component coming from the regional sea-level change and inverse barometer effect. However, these terms are an order of magnitude smaller than the bootstrap-derived component. The uncertainty is then given by:

\begin{equation}
  \delta \dot h = \sqrt{
    (\delta \dot h^*)^2 + (\delta \dot \eta)^2 + (\delta \dot p)^2
    }
  \label{c2e10}
\end{equation}

\noindent
where $\delta \dot h^*$ is the uncertainty in the rate of height change at the 95\% confidence level from the bootstrap distribution, $\delta \dot \eta$ and $\delta \dot p$ are the uncertainties in regional sea-level and regional atmospheric pressure trends, respectively.

Estimating uncertainties using bootstrapping is a `top-down' approach, which has the advantage of not relying on the assumption that characteristics of the noise are known (e.g., normally distributed), or requiring that the relationships between the different sources of error are specified since an algebraic solution of error propagation (`bottom-up' approach) is not required.

\section{Results and discussion}

\subsection{Comparison between backscatter and height}

Our comparison of the height-change time series with changes in backscatter for all ice-shelf grid cells around Antarctica (Fig.~\ref{c2f5}) showed that the majority of the grid cells presented a positive correlation (indicative of densification of the surface) instead of a negative one (indicative of penetration). The backscatter variation impacts significantly the estimated height-changes (correlation $r \approx 0.81$) and the respective correction greatly reduces this effect (correlation $r \approx 0.16$). The anti-correlation between backscatter variation and modeled firn-height changes (correlation $r \approx -0.73$) implies that backscatter increases with firn compaction (densification of the surface), which fluctuates with the seasons. That is, for an increase in backscatter (presumably due to densification) there is an apparent increase in height even though surface densification implies compaction of the firn layer (surface lowering).

\subsection{Trend and acceleration in ice-shelf height change}

We derived two main products from our 18-year long, high-resolution ($\sim$3 months and $\sim$30 km) record of ice-shelf height change.

\paragraph{Regionally-averaged time series.} Regional averaging increases the statistical significance of the information that can be recovered from the available data (e.g., acceleration from the curvature of the fit). We computed circum-Antarctic average time series of height change from 1994 to 2012 for West Antarctic, East Antarctic and all Antarctic ice shelves (Fig.~\ref{c2f7}). These results have revealed accelerated losses in total Antarctic ice-shelf volume at decadal time scale (these results are discussed by \textcite{Paolo2015}).

\paragraph{High-resolution spatial maps.} With analysis at the grid-cell level it is possible to map, for most of the Antarctic ice-shelf area, the spatial structure of the 18-year average rate of change, average acceleration and instantaneous rates in ice-shelf surface height (Fig.~\ref{c2f9}). Mapping higher-order terms such as average acceleration (slope of the derivative of the trend) reveals the regional contribution to the accelerated loss. The primary factor in the overall acceleration of the Antarctic rate of ice-shelf loss is the cessation (or deceleration) of earlier gains in the East Antarctic ice shelves \parencite{Paolo2015}.


\subsection{Processes affecting RA-derived ice-shelf height}

The main sources of uncertainty affecting our ability to use RA-derived height estimates (instantaneous and/or trend) over the Antarctic ice shelves are as follows (the order of magnitude given if for the uncertainty in instantaneous values, unless specified):

\begin{itemize}
  \item[i.]     {\it Surface and volume scattering variation}, $O(1)$ m. Despite considerable progress, ongoing investigation on the interaction of the radar wave with the snow/firn/ice surface and how the surface properties change over time is essential, particularly for the modern CryoSat-2 radar mission and its innovative Doppler-delay altimeter system \parencite{Wingham2006}.
  \item[ii.] 	{\it Firn-column changes}, $O(\text{0.01--1})$ m. This effect is important to separate the mass-loss-driven changes from the surface-compaction-driven changes in the RA observation. However, there is still no reliable model-derived correction to account for this effect at the small spatial scales (tens of kilometers) that we are considering.
  \item[iii.]   {\it Atmospheric-pressure variation}, $O(\text{0.01--0.1})$ m. Since ice shelves are floating, depressions on the sea-surface height forced by changes in atmospheric pressure (inverse barometer effect) will be included in the altimeter ice-shelf-height measurement. There is considerable fluctuation and significant trend in the atmospheric pressure around the Antarctic ice shelves \parencite{Padman2003}, especially on the time scales of a few years that we attempt to resolve here.
  \item[iv.]    {\it Regional sea-level trend}, $O(\text{0.001})$ m/yr. The trend in sea-level rise (which includes changes in the geoid) varies considerably from the global average ($\sim$3 mm/yr) at the regional scale \parencite{Church2004}; and this becomes particularly important when making inferences of ice-shelf height over long time periods (e.g., decades).
  \item[v.]     {\it Advection of surface topographic features}, $O(\text{1--100})$ m. The relatively large footprint size ($\sim$3--5 km) of pulse-limited RAs can lead to moving topographic features (e.g., megaripples and crevasses) being misinterpreted as interannual fluctuations in ice-shelf surface and derived changes in basal mass balance (melting).
  \item[vi.]    {\it Errors in ocean tide and load tide corrections}, $O(\text{0.01--0.1})$ m. Ocean tide height changes can be large, up to $\sim$8 m at spring tides under the southern Ronne Ice Shelf. Standard tide model corrections reduce the tide-induced height error to $<$0.1 m \parencite{King2005, King2011, Stammer2014}. Further reduction by averaging is not effective because the ocean tide is a long-wavelength process. Instead, improvements in tide models are required.
\end{itemize}


\begin{figure}[!hb]
  \centering
  \includegraphics[width=\textwidth]{img/map_4panels_v5.png}
  \caption[Maps of rate and acceleration of Antarctic ice-shelf height]{
  \ssp \footnotesize
  Maps of rate and acceleration of Antarctic ice-shelf height change. (top) 18-year average rate and acceleration of ice-shelf height change; (bottom) instantaneous rate of change in surface height for 2000 and 2012. Rates and acceleration were derived from polynomial models, such that (a) negative acceleration represents either increase loss-rate or decrease gain-rate (and vice-versa for positive acceleration), and (b) instantaneous rates reflect the decadal fluctuations.
  } 
  \label{c2f9}
\end{figure}


\subsection{Methodological limitations}

We are unable to sample near the ice-shelf grounding line. Our grid-cell size is limited by the spatial distribution of the satellite ground-track crossing points (Fig.~\ref{c2f1}). We have eliminated all crossovers within 3 km of the grounding line where data may be corrupted due to the floating-to-grounded transition (high slopes, crevasses, and flexure so that hydrostatic corrections for ocean height do not work). These sampling issues are most severe for small ice shelves that have large melt rates in the grounding zone (such as Pine Island and Dotson), and they imply that estimated ice-shelf loss could be higher than our analyses suggest.

During the observation period, some ice shelves experienced significant grounding-line retreat \parencite{Rignot2014}, without a commensurate retreat of their ice fronts, meaning that those ice shelves increased in area. This ice-shelf migration could impact the average value of volume loss if the new ice-shelf thickness is taken into account, although these area changes are small relative to the total ice-shelf area (less than~7\%). Since our objective is to map and quantify the varying impacts of the ocean and atmosphere to ice loss around the Antarctic ice shelves, we used a fixed area approach (i.e., Eulerian reference frame).

\section{Conclusions}

We have presented a method for improved analyses of height data from multiple satellite RAs and quantifying uncertainties in the resulting data products. Using this method, we constructed an 18-year time series of height changes at $\sim$3 month intervals and $\sim$30 km grid cells over Antarctica's floating ice shelves. Our data set allowed us to estimate the temporal progression and spatial structure of ice-shelf height changes in Antarctica between 1994 and 2012.

We have demonstrated the following. Substantial averaging both in time and space is required to construct reliable RA height records over floating ice shelves. Densification of the surface strongly affects the height-change estimate, and the backscatter correction significantly reduces this effect. Densification is a more important effect than penetration in biasing the height-change estimates over the ice shelves. Given the high interannual-to-decadal variability that is present, a simple straight-line fit fails to capture the underlying trends, some degree of curvature in the trend is needed. Polynomial trends allow us to obtain information on the evolution and spatial structure of changes, e.g., instantaneous rate of change (derivative of the trend), and average acceleration (slope of the derivative). Given the convoluted nature of the error sources, we propose that a top-down approach for uncertainty estimation (e.g., bootstrap applied to time-dependent data) constitutes a more accurate alternative for error analysis.

\clearpage
\section*{Acknowledgements}

This work was funded by NASA [awards NNX12AN50H 002 (93735A), NNX10A-G19G and NNX13AP60G]. This is ESR contribution XXX. We thank J. Zwally's Ice Altimetry group at the NASA Goddard Space Flight Center for distributing their RA data sets for all satellite radar altimeter missions (\url{http://icesat4.gsfc.nasa.gov}). We are grateful to G. Hyland of the Australian Antarctic Division in Hobart for running the Davis S/V algorithm used to make Fig.~\ref{c2f4}. We thank C. Davis and D. Wingham for RA processing advice. We thank A. Shepherd for helpful comments that improved the manuscript. We also thank the San Diego Supercomputing Center for their assistance in HPC. 


{\sl Chapter 2}, in full, is in revision for publication of the material as it
may appear in {\it Remote Sensing of Environment} 2015. Paolo, Fernando S.;
Fricker, Helen A.; Padman, Laurie. The dissertation author was the primary
investigator and author of this paper.