Application of toPEX Altimetry for Solid Earth deformation Studies

This study demonstrates the use of satellite radar altimetry to detect solid Earth deformation signals such as Glacial Isostatic Adjustment (GIA). Our study region covers moderately flat land surfaces seasonally covered by snow/ice/vegetation. The maximum solid Earth uplift of ~10 mm yr-1 is primarily due to the incomplete glacial isostatic rebound that occurs around Hudson Bay, North America. We use decadal (1992 2002) surface height measurements from TOPEX/POSEIDON radar altimetry to generate height changes time series for 12 selected locations in the study region. Due to the seasonally varying surface characteristics, we first perform radar waveform shape classification and have found that most of the waveforms are quasi-diffuse during winter/spring and specular during summer/fall. As a result, we used the NASA β-retracker for the quasidiffuse waveforms and the Offset Center of Gravity or the threshold retracker for the specular waveforms, to generate the surface height time series. The TOPEX height change time series exhibit coherent seasonal signals (higher amplitude during the winter and lower amplitude during the summer), and the estimated deformation rates agree qualitatively well with GPS vertical velocities, and with altimeter/tide gauge combined vertical velocities around the Great Lakes. The TOPEX observations also agree well with various GIA model predictions, especially with the ICE-5G (VM2) model with differences at 0.2 ± 1.4 mm yr-1, indicating that TOPEX has indeed observed solid Earth deformation signals manifested as crustal uplift over the former Laurentide Ice Sheet region.


IntroductIon
Satellite radar altimetry, which was initially designed for accurate measurements of sea surface height, ocean circulation and sea level, has been demonstrated to be applicable to non-ocean surfaces as well.It has been shown that radar altimetry is capable of measuring surface elevations and changes over ice sheets, including Antarctica and Greenland, and subsequently the ice sheet mass balance and its role in global sea level change (Zwally et al. 1989;Wingham et al. 1998).Furthermore, through the analysis of radar waveforms, the profile of backscattered power and the geophysical information related to the near-surface properties of the ice sheet such as the extinction coefficient can be derived (Davis and Zwally 1993).In addition, satellite radar altimetry has been used to measure inland water level variation for hydrologic studies (Birkett et al. 2002;Frappart et al. 2006).Whereas these studies focused on large river basins such as the Amazon, attempts have been made to measure water level change over vegetated wetlands using TOPEX/ POSEIDON radar altimetry (Ibaraki et al. 2006).
The performance of the radar altimeter over land surfaces is significantly different from that over the ocean surface primarily due to different surface roughness and elevation variability.Due to these limitations, there have been attempts to generate digital elevation models (Hilton et al. 2003) or use only the backscattering coefficients from satellite altimetry to classify surface properties (Papa et al. 2002(Papa et al. , 2003)).
Here, we attempt to detect solid Earth deformation due to Glacial Isostatic Adjustment (GIA) using the land eleva-tion measurements from the TOPEX radar altimeter, which is equipped with a circularly polarized antenna and has observed the Earth on a decadal time scale (1992 -2002).GIA produces deformation due to the viscoelastic response of the mantle to the deglaciation of the Laurentide Ice Sheet.The study area covers a relatively flat land surface where the maximum solid Earth uplift reaches ~10 mm yr -1 near Hudson Bay, the location of the former Laurentide Ice Sheet.Radar altimetry has been used to generate digital elevation models with an accuracy on the order of several meters, but was mostly limited to extract elevations along profiles.There has not been any prior attempt to compute spatiotemporal changes of land surface elevations using satellite radar altimeter measurements.In this paper, we describe the waveform classification and retracking methodologies towards generating decadal time series of land surface changes as observed by the TOPEX/POSEIDON radar altimeter.These measurements are then compared with predictions using contemporary GIA models, and other independent measurements.

dAtA
TOPEX/POSEIDON is the first dual-frequency (Kuand C-band) radar altimeter space mission designed to accurately measure global ocean topography with an accuracy of several cm (averaged over 3 seconds) and with 10 days temporal sampling.The TOPEX/POSEIDON satellite was launched on August 10, 1992, on a near-circular orbit repeating every 9.9 days.The orbit inclination is 66° that enables the global observation of the ocean within ±66° latitude bounds.In this study, data from cycles 9 through 364 of the TOPEX Geophysical Data Record (GDR) and the Sensor Data Record (SDR) are used.The Ku-band ten-perframe (10-Hz) surface height data with an along-track spatial resolution of ~660 m are available from the GDR/SDR files.Geophysical corrections such as solid Earth tide and media corrections including dry and wet troposphere delays and ionosphere delay have been applied to the data (see next section).The 10-Hz geodetic coordinates are computed using the Precise Orbit Ephemeris (POE) data provided by NASA Goddard Space Flight Center (GSFC) and tenper-frame time tags are calculated from twenty-per-frame ranges contained in the SDR.The SDR files also contain the ten-per-frame 64-sample waveforms.

MEthodology
First, we search for moderately flat (slope < 0.1°) topographic surfaces as study regions around the Hudson Bay using the ten-per-frame (high-rate) ranges from the TOPEX GDR.Excluding small lakes or potholes around the Hudson Bay, we accept the GDR data provided that each observational point contains 4 -15 ten-per-frame range mea-surements with standard deviations less than 40 cm.In this paper, we intend to show the analysis of the land surface height (LSH) change time series at each of the 12 selected observational points over relatively flat topographic regions near Hudson Bay (Fig. 3) as a demonstration of the validity of our methodology.It should be noted that we anticipate additional observational points could potentially be found in a future study which is out of the scope of this paper, e.g., for a detailed study of the Laurentia GIA phenomena, or other solid Earth deformation signals over other selected regions of the world.
Geophysical corrections such as dry troposphere, solid Earth tide and pole tide corrections have been applied using the values provided in the GDR.Although TOPEX/ POSEIDON is equipped with a three-channel (18, 23, and 37 GHz) microwave radiometer for measuring integrated water vapor contents, interpolated wet troposphere correction based on the European Center for Medium Range Weather Forecasting (ECMWF) model is used in this study since the wet radiometric observations at these frequencies are overwhelmed by the opaqueness of the land surface, rendering them not usable.The decision whether to use the (nadir) dual-frequency ionosphere correction or the DORIS (relative) ionosphere correction is based on a variance analysis assessing the residual RMS of the linear fit to the high-rate range measurements.The variance reduction study result (Table 1) indicates that one should choose either the DORIS or the dual-frequency ionosphere correction depending on the residual RMS of the fit, and that one should not assume that the DORIS ionosphere correction has to be used over non-ocean surfaces (Birkett 1995).

Waveform Shape Analysis and classification
Since the study region near Hudson Bay is not a surface of homogeneous characteristics such as permanently ice covered or inland water surface, in the first step, the classification of the waveform shapes using the so-called waveform shape discriminator is performed.It is based on the specularity and the power distribution of the waveform (Callahan 1992).It was originally developed to affirm that the waveform is obtained over ocean or large water body surface, and not corrupted by land returns.The specularity check compares the power in the 60 th waveform sample with the maximum power, and the power distribution check compares the sum of power in the waveform samples from 5 to 24 with the sum of power in samples from 25 to 60.If both of the ratios are less than empirically predetermined threshold values, which are 0.23 and 11, respectively, the waveform is classified as specular.If one of the ratios is larger than the threshold values, then it is classified as quasi-diffuse.The term quasi-diffuse is used here since it has a more rapidly falling trailing edge than the waveform obtained from openocean, but shows a broader high peak than the specular waveforms.It is found that most of the waveforms are quasi-diffuse during the winter, and as the summer begins, the number of specular waveforms starts to increase.Figure 1 shows the distribution of the waveform shapes obtained from TOPEX repeat cycles 9 to 364 (1992 -2002) for one of the observational points in the study region.
From the waveform shape distribution, we can usually predict that the snow and/or ice physical surface exhibit the quasi-diffuse waveform, since it is expected that the increase of snow/ice accumulation contributes to an increased number of scatter points, causing the waveform to have a more diffusive shape.On the other hand, the specular returns are likely to be from the surfaces with vegetation and/ or soil, which contain irregularly back-scattering points, and they should be dominated by the power signals returned from a few brightest points on the back-scattering surface, leading to a narrow-peaked waveform.Here, although we provide a speculation of the land surface types resulting in the TOPEX waveform shapes, it can be qualitatively confirmed.It will be shown later that the surface height changes agree with seasonal changes in the Northern Hemisphere as our study region primarily contains ice surface in the winter/ spring, and vegetation/soil surface in the summer/fall.The waveform classification is important as it provided more confidence in the radar waveform retracking results in the study, as there is no available in situ data to independently validate our results.

retracking
The performance of a radar altimeter over varying terrain differs significantly from that over ocean surface, and may lead the altimeter's onboard tracker to fail in precisely predicting the range.Therefore, the altimeter range measurements over non-ocean surfaces must be corrected for the deviation of the mid-point of the leading edge from the tracking gate of the onboard tracker (for TOPEX 64-sample waveform, 24.5).This procedure is known as retracking.The retracking algorithms for non-ocean surfaces can be divided into threshold or functional fit (Zwally and Brenner 2001).They have originally been developed for ice sheet altimeter data retracking for mass balance studies (Martin et al. 1983;Bamber 1994;Davis 1997), and have now been extended to coastal circulation and marine gravity studies (Anzenhofer et. al 2000;Deng and Featherstone 2006), and to hydrologic studies for large inland water bodies (Frappart et al. 2006).
Table 1.Spatially and temporally averaged RMS residuals (cm) of a linear fit to the high-rate measurements for each of the 12 study sites.The selected ionosphere corrections and the retracking algorithms used for the specular waveforms are marked as bold numbers.Since the retracking correction essentially has the same role as the acceleration correction and the corrections of the significant wave height (SWH) and attitude effects for ranges predicted by the onboard tracker that is based on the so-called α-β filter, these corrections should be added back to the surface heights before applying the retracking corrections (Brooks et al. 1997;Zanife et al. 2003).However, the SWH/attitude correction, provided in the GDR is not calculated over the non-ocean surfaces and set to a default value in the TOPEX ground processing, meaning that it is not available.Concerning the acceleration correction, which is not provided in the GDR, we estimated the height acceleration using the same algorithm as in the TOPEX ground segment.It is found that the estimated accelerations and subsequently the acceleration corrections are flagged as bad by the GDR algorithm over our study regions.It is because most of the twenty-per-frame (20-Hz) ranges, used as input for the acceleration computation, are not valid over non-ocean surfaces.Therefore, we can directly use the tenper-frame (10-Hz) surface heights contained in the GDR without nullifying any of these two corrections for ranges of the onboard tracker.
Prior to retracking, to mitigate the waveform anomalies such as zero-leakage and the offset leakage effects, we employed the sets of multiplicative and additive waveform factors (Hayne et al. 1994).These factors are originally developed for the TOPEX Side A altimeter, but also applicable to Side B (D. Hancock, personal communication, 2006).In addition, after the retracking, we checked the estimated retracked gate to see if it lies inside the Automatic Gain Control (AGC) gate (9 to 40 for 64-sample TOPEX waveform), which represents the average signal level across the leading edge of the waveform (Chelton et al. 2001).Finally, we have not accounted for the TOPEX Side A altimeter point target response (PTR) degradation problem after TOPEX cycle 140, and definitely after cycle 167 (P.Callahan, personal communication, 2007), which could decrease the accuracy of the land elevation time series.

retracking Algorithms (1) nASA V4 (β-) retracker
The β-retracker is the first retracking algorithm developed to obtain corrected ranges from SEASAT-1 radar altimeter over the Antarctic and Greenland continental ice sheets (Martin et al. 1983).This algorithm uses 5-or 9-parameter functions to fit a single-or double-ramped waveform, respectively.The double-ramped waveforms can be found when two distinct surfaces at different elevations exist within the range window.It should be noted that we retracked only single-ramped waveforms to consistently observe one distinct surface.The Ice Altimetry Group of NASA's Goddard Space Flight Center (GSFC) has developed algorithms for retracking polar ice altimetry based on Martin's functions.There have been four versions of the GSFC retracking algorithms (Zwally 1996), and we used Version 4 which employs an exponential function instead of a linear function to fit a fast-decaying trailing edge which is commonly found in our case: where (2) offset center of gravity (ocog) retracker The OCOG algorithm was developed as an empirical technique to produce ice sheet data products from ERS-1/2 radar altimetry.It calculates the center of gravity, amplitude, and width of a rectangle using full waveform samples (Bamber 1994).The squares of the sample values are used to reduce the effect of low amplitude samples in front of the leading edge (ibid).In addition, the waveform samples 45 to 50 are excluded to avoid the leakage effects of TOPEX waveforms (Hayne et al. 1994).where P(i) is the waveform sample value at the i th bin, and n a is the number of aliased sample (n a = 4 for TOPEX).Finally, the leading edge position (LEP) is given by, (3) threshold retracker The threshold retracking algorithm was developed primarily to measure ice sheet elevation change (Davis 1997).
The leading edge position is determined by locating the first waveform sample to exceed the percentage (i.e., threshold level) of the maximum waveform sample amplitude.The pre-leading edge DC level is computed by averaging the waveform sample 5 to 7, and again, the samples from 1 to 4, 45 to 50, and 61 to 64 are excluded.Davis (1997) suggests the 50% threshold for surface-scattering dominated waveforms, and 10% or 20% threshold level for volumescattering surface.In this study, we adopted the 50% threshold retracker.

combining retracking Algorithms
The waveform shape analysis described in Section 3.1 leads us to believe that there is a need for combining different retracking algorithms, which have their respective advantages and disadvantages.Although the threshold retracker is the most simple algorithm, we used the NASA V4 (β-) retracker for quasi-diffuse waveforms because it is based on the principle of fitting the waveform shape using a physical model (Brown 1977).On the other hand, the 50% threshold or OCOG retracker are used for the specular waveforms because the β-retracker fails to perform for the narrow high-peaked waveforms.The optimal combination of the retracking algorithms (β-retracker with threshold or β-retracker with OCOG) as well as the ionosphere correction is selected based on the calculated RMS residuals of a linear fit to the high-rate surface heights as a by-product of data compression to 1-Hz surface heights.Then, the spatially averaged RMS surface height residuals over each selected study site are again averaged from TOPEX cycles 9 to 364 to represent the uncertainty of candidate sets of retracking algorithms and the ionosphere correction.Table 1 shows the analysis performed to decide which retracking algorithm and which ionosphere correction (DORIS or dualfrequency) are to be used for each study site (Fig. 3).
Finally, a TOPEX time series is generated using the spatially averaged (or compressed) retracked measurements.An example of a time series over one of the study regions is illustrated in Fig. 2. It is generated by combining two time series generated from specular and quasi-diffuse measurements, respectively.Figure 2 shows that the increasing surface height during the winter season corresponds to the increasing snow/ice depth, whereas the decreasing surface height during the spring season corresponds to the melting of the accumulated snow/ice, which further validates our hypothesis of the surface classifications.

rESultS
A collinear analysis has been chosen for the relatively smooth terrain surfaces in this study.This is because the crossover technique, normally used in applications of radar altimetry over ice sheets, produces results over test areas of a limited spatial coverage.It is noted, however, that there is also a recent study to use along-track repeat altimetry over ice sheets to avoid the spurious differences between ascending and descending track heights at the crossover points.However, this only worked for different radar altimeters (ESR-1/-2) with non-circularly polarized antennae (Legresy et al. 2006), and not for TOPEX.As indicated earlier, this paper intends to demonstrate the methodology and constructs time series data at 12 observational points over a moderately flat land area for the data analysis.Slope error correction, another critical correction for land altimetry data, is not considered in this study because it is required only for absolute height determination (Anzenhofer et al. 2000) and presumably the gradient of the terrain surface would not change significantly with time.One issue that must be considered is the effect of the surface gradient in generating a time series of the surface height differences as the horizontal location of each observation point changes from cycle to cycle.Since we have selected relatively flat surfaces around Hudson Bay as study regions, we can model the surface as a plane and the surface gradient can be computed from the satellite observations, i.e., the retracked surface heights.The detailed gradient estimation algorithm can be found in Guman (1997), and a brief description is given in this paper.It should be emphasized that this approach will be valid only over reasonably flat surfaces.Our local mean land surface in a study region (say, bin) is modeled with along-and cross-track gradients as,

LSH a b dx c dy
where LSH is: mean land surface height; a: height of the plane at the bin center; b: along-track land surface gradi-Fig.2.An example of a selected time series generated by combining the specular and quasi-diffuse retracked measurements.
ent; c: cross-track land surface gradient; dx: along-track displacement of a data point from the bin center; dy: crosstrack displacement of a data point from the bin center.The parameters a, b, and c can be estimated by the least squares adjustment.However, as can be seen from Fig. 2, in addition to the spatial variation, the land surface height also exhibits temporal changes from cycle to cycle.Thus, we adopted a model with six parameters that includes linear, annual, and semi-annual variations. [ where A is: the offset or bias; B: the linear slope; c1, c2: the amplitude of the cosine term; s1, s2: the amplitude of the sine term; ω: the annual frequency.
To reduce the effect of the terrain surface gradient on the estimate of the annual/semi-annual variations, equations ( 8) and ( 9) are implemented in an iterative scheme.We first estimate a preliminary annual/semi-annual variability using equation ( 9), and estimate the preliminary land surface gradient using the annual/semi-annual variation-removed surface height.This procedure is iterated once more to yield a gradient-corrected surface height.Table 2 shows the variation (standard deviation) of the amplitude of the time series before and after the gradient correction.
Due to the varying amplitudes of the time series, the Multi-Taper Method (MTM) (Ghil et al. 2002) is used instead of the typical least squares approach to estimate the uplift rate from the generated time series.We speculate that the seasonal signals shown in Fig. 2 are due to the varying radar scatters and/or the penetration depths of two different surfaces (ice or vegetated land).The goal is to remove the seasonal signals (with varying amplitudes and possibly varying frequencies) while better preserving the linear trend signal.The MTM spectral analysis method provides a novel approach for spectral estimation of a time series, which is believed to exhibit spectra containing both continuous and singular components, and signal reconstruction from selected spectral components.This method has been widely applied to problems in geophysical signal analysis (Ghil et al. 2002;Kuo 2006).We first estimate the spectra of the time series and then reconstruct the contributions of selected components of the time series.Finally, a linear trend can be estimated by linear regression using the reconstructed time series.Figure 3 illustrates the generated time series with their estimated linear trends over the 12 selected observational points within the study region.
Figure 3 shows that the selected data points south of Hudson Bay are uplifting at a rate of approximately 8 to 13 mm yr -1 , whereas the study sites near the Lake Superior and Lake Ontario show smaller and negative uplift rate, respectively.The result qualitatively agrees with computed uplifts using various GIA models, which predict that the maximum uplift occurs around the Hudson Bay and minimum or subsidence near or south of the so-called transition zone which is located around Lake Erie and the southern Lake Michigan.Figure 4 shows the uplift rates from other techniques including the estimated vertical motions at water level gauges around the Great Lakes by combining TOPEX/ POSEIDON altimetry and long-term water level gauge records (Kuo et al. 2008), and the GPS vertical velocities from National Geodetic Survey (NGS) solution (Snay et al. 2007; R. Snay and M. Cline, personal communication, 2006).Another two GPS solutions at Churchill and Kuuj (Calais et al. 2006) are included to show the qualitative agreement between the GPS velocities and the nearby TOPEX altimetry observed land uplift rates.The Kuuj GPS site solution is less accurate because of the limited data span (~3 years).
GIA models used for qualitative comparisons with the TOPEX observations include the ICE-4G (VM2), ICE-5G (VM2) (Peltier 2002(Peltier , 2004;;data courtesy, W. Peltier 2006), the simple three-layer viscous model (hereafter TLM) used in the BIFROST GPS/GIA study (Milne et al. 2001;data courtesy, J. Mitrovica and G. Milne 2004), and the L20 model (Wang and Wu 2006;data courtesy, P. Wu 2006).The first three models are radially symmetric while the L20 is a 3-D (incompressible) model with laterally heterogeneous viscosity profiles and lithospheric thickness using realistic Earth model parameters.The TLM consists of sev-  eral parameters with different lithospheric thickness (LT), upper mantle viscosity (UMV) and lower mantle viscosity (LMV).The combinations of Earth parameters used in this study are LT = 120 km, UMV = 0.5, 1 (× 10 21 Pa·s), LMV = 1, 3, 10 (× 10 21 Pa·s) (see Table 3).It should be noted that the TLM model may not have been optimized over the Laurentia region, however, as it is shown later, the model (LT = 120 km, UMV = 1 × 10 21 Pa·s, and LMV = 3 × 10 21 Pa·s) has the same qualitative agreement (Fig. 5) with the TOPEX observations as the ICE-5G or the L20 models.The ICE-3G (Tushingham and Peltier 1991) and ICE-4G ice loading histories were adopted to generate the TLM and L20 models, respectively.
From Table 3, it can be seen that the estimated TOPEX land deformation rates show good agreement with various GIA model predictions, in particular with respect to ICE-5G (VM2).The mean and standard deviation of the differences are 0.2 and 1.4 mm yr -1 , respectively, indicating an excellent agreement.In addition, the difference between TOPEX results and the L20 laterally varying viscosity and lithospheric thickness model is 0.2 ± 2.0 mm yr -1 .The ICE-4G (VM2) model is negatively biased with respect to other GIA models and TOPEX observations (Fig. 5).It is also remarkable to  see from Fig. 5 that the altimetry-derived vertical velocities follow the pattern of variations predicted by several GIA models.However, while the comparison shows good agreement between TOPEX observations and models, it is largely a qualitative comparison as the observation samples (12 sites) are quite limited.

concluSIonS
This paper describes an efficient retracking approach to enable TOPEX radar altimeter measurements over land for observations of solid Earth deformations around Hudson Bay.The estimated vertical land motion is primarily attributed to GIA because of the well-known incomplete glacial isostatic rebound of the mantle as a result of the Laurentide Ice Sheet deglaciation since the Last Glacial Maximum, 18000 -20000 years before present.The estimated TOPEX land uplift rates agree well with the combined TOPEX/PO-SEIDON and tide gauge measurements (Kuo et al. 2008) and the GPS vertical velocities (Snay et al. 2007;Calais et al. 2006).The results further agree well with various GIA models.It was demonstrated that radar altimetry can be used for geodynamics studies, in this case, for detecting the GIA signal over land.Future studies must include improvement of the retracking algorithm to obtain more data points over potentially steeper regions and also use other satellite radar altimeters such as ENVISAT and GFO, and ICESat laser altimetry, in order to increase the spatial coverage and to expand the measurement time span.It is anticipated that this technique could potentially be applied to other regions of the world and detect other solid Earth deformation signals if the vertical motion rates are comparable in magnitude to GIA.

Fig. 1 .
Fig. 1.Waveform shape distribution for a selected data point in the Hudson Bay study region.

Fig. 3 .
Fig. 3. TOPEX cycle 343 ground tracks with retracked ellipsoidal heights (middle).Time series are generated over the study sites (red triangle) using different retrackers chosen based on the shape of the waveforms.The linear trends are estimated from the Multi-Taper Method (red straight line).

Fig. 5 .
Fig. 5. Comparison of TOPEX observed land vertical motion and GIA models at 12 study sites.

Table 2 .
Result of the surface gradient correction shown with standard deviation (cm) of the amplitude in the time series.

Table 3 .
Comparison of the TOPEX observed land uplift rates with GIA models.