Interferometric swath processing of Cryosat data for glacial ice topography

We have derived digital elevation models (DEMs) over the western part of the Devon Ice Cap in Nunavut, Canada, using “swath processing” of interferometric data collected by Cryosat between February 2011 and January 2012. With the standard ESA (European Space Agency) SARIn (synthetic aperture radar interferometry) level 2 (L2) data product, the interferometric mode is used to map the cross-track position and elevation of the “point-of-closestapproach” (POCA) in sloping glacial terrain. However, in this work we explore the extent to which the phase of the returns in the intermediate L1b product can also be used to map the heights of time-delayed footprints beyond the POCA. We show that there is a range of average cross-track slopes (∼ 0.5 to∼ 2) for which the returns will be dominated by those beneath the satellite in the main beam of the antenna so that the resulting interferometric phase allows mapping of heights in the delayed range window beyond the POCA. In this way a swath of elevation data is mapped, allowing the creation of DEMs from a sequence of L1b SARIn Cryosat data takes. Comparison of the Devon results with airborne scanning laser data showed a mean difference of order 1 m with a standard deviation of about 1 m. The limitations of swath processing, which generates almost 2 orders of magnitude more data than traditional radar altimetry, are explored through simulation, and the strengths and weaknesses of the technique are discussed.


Introduction
Surface elevation data for ice caps and glaciers derived from satellite altimetry have been used to produce both surface elevation DEMs (digital elevation models; e.g.Bamber et al., 2009) and measurements of surface elevation change (e.g.Zwally and Giovinetto, 2011).The estimation of surface elevation change forms the basis of one of three remote sensing methodologies for calculating ice sheet and glacier mass balance information (Shepherd et al., 2012;Gardelle et al., 2012).Continued monitoring of glacial ice is important for understanding and predicting the potential impact of global warming on sea level rise.
In the past, satellite radar altimetry has not provided reliable elevation data over small ice caps, like the Devon Ice Cap in Canada, and had limited utility over the sloping edges of the large ice caps.Cryosat (CS) was developed by the European Space Agency (ESA) and launched in 2010 to provide radar altimetry data specifically designed to provide improved elevation data over the edges of ice sheets and smaller ice caps, and for measuring sea ice parameters (Wingham et al., 2006).
The CS SAR/Interferometric Radar Altimeter, SIRAL, operates in a number of different modes (Wingham et al., 2006), including two which use coherent along-track "delay Doppler" processing (Raney, 1998) to create a relatively small along-track footprint (∼ 300 m).The interferometric SARIn ((synthetic aperture radar interferometry)) mode uses coherent along-track processing of the returns received from two antennas such that interferometric phase difference can be used to resolve and geocode the initial returns in the crosstrack plane (Jensen, 1999).Since launch, the SARIn mode has been switched on extensively over the sloping ice at the edges of the Greenland and Antarctic ice sheets, and over all the other glaciated areas on Earth (Parrinello et al., 2013).
This work used intermediate level (L1b) SARIn product files in which the delay Doppler processing step has been completed to waveform data, but the processing to terrain height has not.The waveform data include the unique "point of closest approach" (POCA) followed in delay time by the sum of returns from both sides of the POCA.With a crosstrack slope the POCA will be displaced from the sub-satellite track and its position and height can be calculated using the interferometric phase (Jensen, 1999).With increasing time delay beyond the POCA the composite footprint will include a region beneath the satellite with strong antenna illumination and a region upslope from the POCA with weaker antenna illumination.This is illustrated in Fig. 1.If the upslope "range-ambiguous" region is illuminated very weakly in relation to the "main-beam" region, then the interferometric phase will define the cross-track look angle to the main-beam footprint and can be used with the delay time to map the position and heights of the main-beam footprints.In contrast to the level 2 (L2) SARIn data product, containing only the POCA height and position, the current analysis can generate an across-track swath of elevations under suitable slope conditions.This processing approach was demonstrated with ASIRAS, the ESA-developed airborne simulator for Cryosat (Hawley et al., 2009).
In Sect. 2 the problem of interferometric processing of data in which the range window includes returns from two distinct areas is addressed.The pre-launch measured antenna patterns are used in Sect. 3 to estimate the relative power in the two range-ambiguous zones.The Devon Ice Cap was selected as one of the Cryosat calibration-validation sites as it has been studied for many years (reviewed in Boon et al., 2010).Section 4 includes an outline of the swath processing steps for the L1b data collected over Devon during 2011.Two important limitations of the swath processing approach are discussed in Sect. 5.In Sect.6 a product simulation approach is described that shows that the DEM produced by combining the swath processed results can be used to check and potentially refine the elevation data.The problems, strengths and potential of the swath processing technique are discussed in the final section.

Interferometry with multiple footprints
Interferometry for civilian Earth observation began in the 1980s (e.g.Zebker and Goldstein, 1986) and evolved rapidly as an important technique for terrain mapping and the measurement of relative terrain displacement (Rosen et al., Cryosat is oriented such that the baseline formed by the two receiving antennas is oriented perpendicular to both the along-track direction and to the normal to the WGS84 ellipsoid.Each crosstrack scan line will contain a unique estimate of the POCA on the ice surface, followed in delay time by the sum of returns from both sides of the POCA.With a suitable geometry the returns from the ambiguous region can be much weaker than those from the main beam underneath the satellite which allows the mapping of position and height of main-beam footprints.The dashed magenta line illustrates diagrammatically the dB variation in cross-track two-way gain.Note that the angles have been exaggerated for clarity; the first side lobe in the pattern occurs at ∼ 2 • from nadir and represents a two-way power level of −40 dB with respect to the power at nadir.

2000)
. In cross-track radar interferometry for mapping, the differential phase between the returns from two antennas is used to calculate the cross-track look angle which, together with the platform position, velocity and range, allows an accurate mapping of footprint height and position.However, when swath processing is used with CS L1b SARIn data we must consider the multiple range-ambiguous regions explicitly.In normal SIRAL operation the left hand antenna (viewed in the direction of motion) is used for transmission and both the right and the left antennas are used for reception by the two SIRAL receivers.
With the delay Doppler processing adopted for Cryosat 64 fore-aft "beams" are generated, registered, and summed to reduce the speckle effect (Wingham et al., 2006;Galin et al., 2013).A comprehensive model of the power and phase of the returns has been developed by Wingham et al. (2004Wingham et al. ( , 2006) ) and used to predict waveform response, particularly for the region close to the POCA.This model was also used in CS SARIn mode calibration (Galin et al., 2013).In the following simplified approach we explore just the implications of the fact that the time-delayed returns beyond POCA The Cryosphere, 7, 1857-1867, 2013 can originate from two totally separate footprints.As delay time increases beyond the POCA time one of the two rangeambiguous footprints in sloping terrain will move towards the sub-satellite track and the other will move away.We refer to the former as the "main" footprints and the later as the "ambiguous" footprints (Fig. 1).Note that the two parts of the composite footprint, the main and ambiguous, may be many kilometres apart but still contribute energy to the same receiver range window.
The processed receiver output connected to the left hand antenna V Li for the ith fore-aft beam at zero Doppler range r L and at an arbitrary point in time (retaining only the phase term related to the radar range between the phase centres of the antennas and footprints) can be expressed as a combination of the main and ambiguous components and is + G L (θ a , ϕ i ) A a σ 0 a (χ a ) exp(j kr La ).G L (θ, ϕ) is the one-way power gain of the left hand antenna at an across-track look angle of θ and an along-track look angle of ϕ.A m and A a are the surface areas of the main and ambiguous surface footprints defined by the pulse compression range resolution, the along-track delay-Doppler resolution and the angle the beam makes with the local surface slope.σ 0 m,a (χ ) is the normalized back scattering coefficient for the main or ambiguous footprints at a surface incidence angle of χ, j = √ −1 and k is the wave number.The range variation in any data set is very small in relation to the satellite altitude so that the r 5/2 variation in return power (Raney, 1998) can be ignored.While the delay timing does guarantee that the main and ambiguous footprints are essentially at the same range, the phase terms exp(j kr Lm ) and exp(j kr La ) representing the range between the phase centre of the left antenna and the phase centres of the main and ambiguous footprints respectively are unknown and independent of one another.
Similarly the receiver output V Ri connected to the right hand antenna for one Doppler beam and range window is proportional to ) exp(j kr Ra ).The interferogram I is formed through the look summation of the product where the * symbol indicates the complex conjugate operation.The phase of the interferogram then contains information on the path length differences, which can allow a geometric solution for the position of the footprint (as in traditional cross-track interferometry).There are four terms in the V Ri .V * Li product.Two of these, the product of the first terms of Eqs.(1) and ( 2) and the product of the second terms will remain coherent as delay time increases beyond the POCA position as these terms reflect the combination of signals from the right and left antennas for the main and ambiguous footprints separately.The main or ambiguous footprint size is such that when it is considered as an antenna the resulting "celestial footprint" subtended at the satellite is very much larger than the baseline between the two antenna phase centres and the "baseline decorrelation" condition for coherent interferometry is well satisfied (Rosen et al., 2000).This is not the case for the two cross terms involving a combination of the first and second terms of Eqs. ( 1) and (2).In this case the interferogram terms are formed from a combination of the main and ambiguous footprints, which are moving apart as delay time increases.Consequently as the separation between the main and ambiguous footprints increases above ∼ 10 km the cross-track celestial footprint created at the satellite will approach and eventually be less than the actual baseline so that the critical baseline condition (Rosen et al., 2000) will not be satisfied.These cross terms are the reason that the coherence often decreases beyond POCA.Note that the subsequent increase in coherence is due to the fact that the differential phase for the 64 individual beams in the cross terms becomes random so that when coherently summed their contribution to the interferogram diminishes.
Consequently, when the two footprints are far enough apart the interferogram I can be approximated by the sum of the two terms from the main and ambiguous beams.
Where the two-way gain R (θ, ϕ) and r m,a are the path differences from the either the main or ambiguous footprints to the two antennas.Note that the magnitude of the interferogram is now in units of power and not backscattered field amplitude.If the antenna gains illuminating the ambiguous footprints are much less than those illuminating the main beam (< −30 dB) then the phase of the interferogram will be approximately that given by the first term in Eq. ( 4).With the current SIRAL transmit-receive configuration, and when the second term in Eq. ( 4) can be ignored, the interferogram phase ψ given in the L1b product file can be related to the cross-track look angle θ to the main footprint by where B is the interferometric baseline, β is the roll angle of the baseline, and θ is positive when the look angle is to the right of the sub-satellite track when viewed in the direction of motion.Knowing the satellite position, the range and the cross-track off-nadir look direction it is then possible to www.the-cryosphere.net/7/1857/2013/The Cryosphere, 7, 1857-1867, 2013 map the surface footprint position and height as in traditional cross-track interferometry.

Suppression of ambiguous returns
We estimate the relative contribution from the main and ambiguous zones as a function of cross-track slope by using prelaunch measurements of the antenna patterns (Saab Space AB, 2007).Figure 2 (top) illustrates the relationship between the circular iso-range contours and the fore Doppler beams for a locally spherical Earth model.The positions of every second Doppler beam are indicated by the dotted magenta lines extending in a grid across Fig. 2 (top).The range and angle coordinates are presented for the situation in which there is no surface slope.However, the relative positions and angles would still apply when there is a surface tilt although in this case the 0.0 point on the plot would be the position of the POCA.In "delay Doppler" processing the slant range variation for across-track footprints is compensated such that the 64 fore-aft beams for one across-track line can be registered and "stacked" to reduce the speckle effect.The spacecraft is accurately attitude controlled to point the peak antenna gain downwards perpendicular to the WGS84 ellipsoid so that, with an across-track slope, the peak antenna gain does not coincide with the POCA.An effective two-way cross-track antenna pattern G 2-way (θ ) appropriate for interferometry was calculated using summation of the pre-launch two-dimensional patterns at the anticipated Doppler beam fore-aft angles (Wingham et al., 2006).
Figure 2 (bottom) illustrates five across-track track cuts of the two-way pattern at five values of the along-track angle ϕ from 0 • to the nominal maximum angle of ∼ 0.76 • as dotted black lines.The solid blue line represents the G 2-way (θ ) sum of the 64 fore-aft two-way pattern cuts.Note that if the POCA is directly beneath the satellite then all 512 samples in range would be illuminated by the main lobe of the twoway pattern.However, with an across-track slope of 1 • , the POCA position is displaced such that the initial returns are subject to much lower two-way gain, and subsequent returns have very different antenna gains between the left and right sides of the POCA.Experience with L1b data over Devon shows that the initial returns from glacial terrain are normally still strong enough for correct operation of the tracking loop, which defines the range window, even at across-track slopes of ∼ 1.5 • .Figure 3 illustrates the relative power of the ambiguous beam in relation to that in the main beam for across-track terrain slopes between 0.5 and 2 • .A circular Earth was used for this simulation with the initial SARIn mode sampling  scheme in which each scan line contains 512 samples of the wave form, coherence and phase separated in time by 3.125 ns. Figure 3 shows that there is a range of slopes and delays for which we anticipate that the amplitude of the ambiguous return would be small enough that the interferogram phase could be used as a first estimate of the angle to the main beam.It is important to emphasize that this is a guide to what is possible; real glacial ice terrain has a relatively smooth but complex surface topography so that the acrosstrack dimension of a main-beam footprint may be quite different from that of the equivalent footprint in the ambiguous beam, which would affect the main to ambiguous ratio.However, the implication of Fig. 3 is that there are a range of given in Eq. ( 6).Note that the cross-track angle extent in Fig. 2 (bottom) is much wider than that in Fig. 2 (top).cross-track slopes and delays from POCA that would potentially support swath processing.

Processing outline and results
We illustrate the steps in the processing by using the descending pass acquired over the Devon Ice Cap on 23 February 2011.In Fig. 4 the sub-satellite descending tracks for the period between February 2011 and January 2012 are shown as dotted black lines and the 23 February pass is shown as a solid black line.The white lines indicate the positions of the airborne laser data acquired during the 2011 CRYOVEX campaign, which have been the primary reference data for the evaluation of the Cryosat elevations.
Cryosat data are read from the ESA-provided SARIn "DBL" files (version 3.9) containing a sequence of crosstrack waveform records separated by ∼ 300 m along-track, each of which contains 512 samples of the power, coherence and phase.These values span a total range window of 240 m which will translate into a ground range swath typically up to ∼ 5 km.In the latest ESA data release the sampling has been increased by a factor of 2 so that the total swath window for 512 samples has been reduced to 120 m.In the L1b file each line record has already been averaged, or "multi-looked", by summing up to 64 registered fore-aft beams.Notwithstanding this averaging, the interferogram output still suffers from both residual speckle and thermal noise.Phase noise has been reduced from that in the original L1b file by recreating the interferogram, smoothing the real and imaginary components on a line-by-line basis and extracting the phase from the smoothed interferogram.This process effectively reduces the inherent slant range resolution from ∼ 0.5 to ∼ 2.5 m.The resulting across-track resolution then depends on the local incidence angle, e.g. at a local incidence angle of 1 • the crosstrack surface resolution is ∼ 143 m.Note that even with the smoothing, the cross-track footprint size for swath processing is normally much smaller than a typical POCA footprint.The resulting power, coherence and phase are illustrated in Fig. 5 as a function of the delay window.
Data were quality controlled interactively and regions in each line with high phase noise or low coherence (< ∼ 0.8) were removed from the analysis.The phase is unwrapped on a line-by-line basis by working in both directions from a central, high coherence region.Profiles of the power, coherence and phase are shown in Fig. 6 at the position of the white arrow in Fig. 5.The part of the line which was used for mapping can be seen as the red section of the phase (3 plot).The smoothed, unwrapped differential phase (shown in red) and information on the interferometric baseline were used to estimate the unit vectors in the "Cryosat reference frame" (Wingham et al., 2006) pointing towards each footprint in the cross-track swath.Using the data provided on satellite position and delay times, the positions and heights of footprints were calculated in a geodetic frame assuming that the differential phase reflected returns from the main beam alone and any contribution from range-ambiguous regions was negligible.The resampling from the irregularly spaced across-track sampling to a regular lat-long grid was completed in two stages: the initial processing leads to a sequence of latitude, longitude and elevation values for that part of each waveform or cross-track scan line which satisfies the coherence requirement.Although the results are irregularly spaced, they span an almost straight line in the latitude-longitude domain.Coefficients for a second order polynomial are obtained from a least squares fit of the latitude to the longitude values.Latitude and elevation data can then be interpolated at a regular sampling in longitude.For the average latitude of the Devon Ice Cap the 0.004 • sampling in longitude corresponds to a spacing of 115 m, which is larger than the normal separation of the original irregular longitude values in swath processing so that the oversampling created by the phase filtering is reduced.The second stage is then the resampling of the elevation from the irregular to regular values in latitude leading to elevations interpolated to a lat-long grid with a sampling of 0.001 • in latitude (equivalent to a spacing of ∼ 112 m).This resampling greatly simplifies the merging of data from differ-1 Fig. 6.Top: profile of the return power for the scan line at the position of the white arrow in the illustration of return power in the upper part of Fig. 5.Note the leading edge is relatively weak, indicative of the low two-way antenna gain at the POCA cross-track angle.The strong delayed returns reflect the fact that the delayed main-beam returns have a much stronger two-way antenna gain.Second from top: profile of the smoothed coherence for the same scan line.Note that the dip in the coherence at ∼ 30 m in the delay window corresponds to the situation when the baseline decorrelation affects the "cross terms" discussed in the text.The high coherence for the range window between ∼ 60 and 180 m implies that the interferogram can be represented by Eq. ( 5) and that the returns are dominated by those from the main footprints.Third from top: the original phase from the L1b file is shown for the scan line in blue.The region of the scan line for which processing was done is shown in red.Note that the phase has been low pass filtered, the boxed region in this figure has been blown up in the bottom figure and illustrates the degree of smoothing; the red line is the unwrapped smoothed phase and the blue dots are the original values.Although it is slope dependent the sampling and resolution in the across-track direction is normally better than in the along-track direction, even with this smoothing.
ent passes but the original geocoded but irregularly sampled elevation data would still be useful to many users.Figure 7 illustrates the resampled elevations as colour and shows that, in this case, results could be obtained over a cross-track distance of approximately 5 km.The results were checked against the reference airborne laser scanner (ALS) data flown over the Devon Ice Cap during the 2011 CRY-OVEX campaign.Because the laser footprint is so small (∼ 0.5 m 2 ) in relation to the Cryosat swath processed footprint (10 4 -10 5 m 2 , depending on the local incidence angle) the ALS data was used to create a reference DEM whose values were weighted averages centred on the lat-long grid used for the Cryosat geocoding.The position of the E-W ALS line is shown in Fig. 7 by the white arrows.The mean laser height minus CS swath processed height for this data set was −0.28 m and the standard deviation of the difference was 0.44 m.If the phase filtering stage was omitted the higher phase noise resulted in position errors for the footprint centres and the standard deviation of the difference increased to 0.74 m.In both cases ∼ 100 height differences were used.
Results from 24 ascending and 25 descending passes were combined separately to create two DEMs.These data were acquired throughout the period from February 2011 to January 2012, usually two or three acquisitions per month for both ascending and descending passes.The DEM illustrated in Fig. 8 was created from a combination of the results from the 25 descending passes and Fig. 9 illustrates the comparison between the CS and airborne laser altimeter derived reference data.The E-W elevation profile at 73.339 N is shown  in the upper plot and the histogram of height differences between the CS and ALS data is illustrated in the middle.The mean difference is 0.49 m and the standard deviation of the height differences is 0.75 m.The lower plot in Fig. 10 illustrates the E-W slope (approximately across-track) and shows that the western slopes of the Devon Ice Cap do match those predicted by Fig. 3 as being appropriate for swath processing.

Impact of interferometric phase and roll angle errors on swath processed height errors
Accurate mapping with CS in the cross-track direction depends on a precise knowledge of the cross-track look angle from nadir, which is estimated from the differential phase and the baseline roll angle.Variable solar illumination leads to some small bending on the satellite such that the star tracker data, which provides the satellite attitude data, may not reflect the exact roll angle of the interferometric baseline (Galin et al., 2013).At POCA, the beam intersects the surface at right angles (local incidence angle = 0 • ) so that while a small roll angle error leads to a cross-track position error, to first order it does not lead to a height error.In essence one gets the correct elevation for the wrong position.However, with swath processing the radar beam will not intersect the cross-track surface at an incidence angle of 0 • so that any roll angle error would result in mapping and height errors.For example, a roll angle error of 1.10 −4 radians (∼ .006resulting bias in height would be ∼ 1.3 m.In summary, the first order error in derived height ε h depends on the error in the cross-track off-nadir angle and the local incidence angle χ local , and is given by In Eq. ( 7) we assume that the off-nadir angles are small enough that sin θ = θ and that any error in range r is negligible.The phase error ε ph can be due to phase noise, phase calibration drift and also through the assumption that the interferogram phase reflects just the main footprints when in fact there may be a small contribution from the ambiguous region.The roll angle error ε roll is not constant and may reflect the changing thermal environment (Galin et al., 2013).As χ local is essentially zero for the POCA this equation implies that swath processing will not achieve the height accuracies possible with the standard geophysical L2 POCA height estimates.This is certainly true for the range of cross-track angles for which the SARIn mode was expected to work (nominally < ∼ 0.6 • , note the phase will "wrap" at an off-nadir angle of ∼ 0.54 • complicating the operational algorithm).However, our experience with the L2 results over Devon, is that when either the across-or along-track angle exceeds ∼ 0.5 •  the leading edge of the POCA return is quite variable and the "retracker" required to identify the POCA position and interpolate the POCA phase may not generate reliable results.In essence, the operational SARIn mode was not designed for larger off-nadir angles.This shows that swath processing can complement the L2 results by providing heights in an acrosstrack off-nadir angle regime, which is outside the range for which we expect good L2 POCA results.
To minimize a possible roll angle error the results from each CS pass were compared to reference heights derived from the ALS flown over the Devon Ice Cap during the CRY-OVEX campaigns.The baseline roll angle was adjusted such that the derived E-W slope better matched the reference data slope.For example, a roll angle offset of 0.00023 radians was used to process the 23 February data leading to a difference standard deviation of 0.43 m.However, if no roll offset was used there was a tilt in the difference results in the cross-track direction and the standard deviation increased to 0.88 m.In general the roll angle correction was typically of order 2.10 −4 radians but the uncertainty in this value does lead to potential bias errors in the absolute swath processed height estimates.In fact we cannot be certain that the roll angle correction is due purely to a roll angle error; there could be a small ambiguous zone component to the interferometric phase.If the procedure to minimize the tilt between the reference data and the CS results was omitted in processing all the data, i.e. the satellite roll angle was used without any correction, the average difference between the resulting Cryosat DEM and the reference data increased to ∼ 2.6 m but the standard deviation about the mean only increased from 0.75 to 0.83 m.This area is still under active investigation.
The impact of an unknown ambiguous zone contribution is not easy to anticipate a priori but the use of data simulation using the derived DEM does help evaluate the results; this is discussed further in the section below.The worst case height error occurs when the phase of the ambiguous zone contribution is some multiple of π/2 with respect to the phase of the main-beam contribution.If the ambiguous power is 30 dB below that in the main beam then the worst case error in the look angle is ∼ 9.5 10 −5 radians leading to a height error of 38 cm if the local incidence angle in the main beam is 0.5 • .

Data simulation
Cross-track profiles of elevation as a function of position can be estimated from the CS generated DEMs for each scan line.By using the satellite position and timing information this allows the estimation of look angle, cross-track slope, local incidence angle and footprint size for each of the cross-track footprints.Note that this can be done for both the "main" and "ambiguous" beams separately.Also, the relative power in the main and ambiguous beams can be estimated using the footprint size, the local incidence angle and the antenna gain patterns.The relative power for the 23 February data is illustrated in the upper panel of Fig. 10 and the simulated power for the main beam is illustrated in the lower panel.This was based on the size of the footprint in the main beam and a relative backscatter model linear in local incidence angle (−2.5 dB times the local incidence in degrees) and clearly captures the main features of the actual data.In the same way, the relative received power in the range-ambiguous beam can be compared to that in the main beam directly beneath the satellite.This shows that for much of the received swath the contribution of the upslope range-ambiguous region is very small (−30 to −40 dB with respect to the main beam).Also, with the knowledge of the look angles it is possible to estimate the phase as well as the magnitude of the contribution from the ambiguous region so that regions with potentially significant errors can be removed.However, with the crosstrack slopes on the western flank of the Devon Ice Cap, the minimum coherence criteria used already eliminates much of the areas affected by the ambiguous regions upslope.A re-examination of the processed data shows that while the simulation approach does improve results through the elimination of some poor regions the problem of small bias errors can remain.

Discussion and conclusions
The increased data volume possible with swath processing of Cryosat SARIn L1b data over sloping glacial ice allows the creation of DEMs where previously radar altimetry provided only a sparsely sampled along-track elevation profile often with some uncertainty as to the position of the footprint (see e.g.Hurkmans et al., 2012).Even at latitude ∼ 75 • N the available data from one year over the Devon Ice Cap was sufficient to allow the creation of DEMs from both ascending and descending passes.At higher latitudes the sub-satellite tracks converge and the task becomes easier.
However, there are problems that do not exist with the conventional processing to the L2 SARIn products.Firstly, the time-delayed footprints beyond POCA contain contributions from at least two footprints and useful results can only be generated when the returns from one contiguous footprint are much stronger (> ∼ 30 dB) than any other footprint in the same time-delay window.This restricts the method to areas with an average cross-track slope of ∼ 0.5 to ∼ 2 • .Secondly, swath processing places stringent requirements on the baseline roll angle and phase calibration, which are difficult to satisfy.In particular, the uncertainty in the baseline roll angle will lead to biased height errors.Careful comparison of the CS DEM results with our various surface elevation data shows that spatially varying bias errors of order 1 m can exist.Further, we find that there appears to be a bias again of order 1 m between the DEM derived from the descending passes and that derived from the ascending passes.The reason for this difference remains uncertain.However, one advantage of swath mode processing is that the complexity associated with interpolating the surface return position (the "retracker") is avoided.
Once a DEM of the area has been completed it is possible to use the simulation approach to estimate the contribution from the various footprints and to further refine the results through the elimination of any data, which may have significant errors.Our results over Devon also show that results are worse when there is an average along-track slope > ∼ 0.5 • , or when there is significant topographic variation.For example it was much more difficult to obtain useful results on the eastern flank of the Devon Ice Cap which has a much more varied topography.
While the results to date are not good enough to map the seasonal time rate of change of Ku band radar altimeter elevation due to snow accumulation, surface melt, firn compaction, etc., further work will explore the limits of the technique and hopefully improve the absolute accuracy.Nevertheless, the demonstration that multiple passes can be combined to create DEMs of sloping glacial ice terrain is of value.It is anticipated that this approach could be used over at least parts of the major outlet ice streams and glaciers in Greenland and Antarctica covered by CS SARIn data.Improved surface elevation and slope data is important for a better understanding of ice dynamics: it will enable better L. Gray et al.: Interferometric swath processing of Cryosat data mapping of the slope in the direction of movement (related to the driving stress) and also the mapping of surface elevation "dips" in ice streams and glaciers such as the elevation results derived from 19 Cryosat descending passes for the subglacial Lake Mercer (SLM) and shown in Gray et al. (2013).Periodic filling and draining of SLM and many other subglacial lakes was documented by Fricker et al. (2007) andSmith et al. (2009) with ICESat data after the initial observation of the associated vertical surface movement through 3-D InSAR work (Gray et al., 2005).The increased data volume possible with CS swath processing combined with the improved spatial resolution will help in mapping the surface subsidence or inflation associated with the periodic filling or draining of subglacial lakes.In this regard McMillan et al. (2013) have already used Cryosat data to measure the ice surface uplift subsequent to the drainage of the Cook E2 subglacial lake originally identified by Smith et al. (2009).Although they used the Cryosat level 2 elevation results, the improved resolution and ability to map the cross-track position of the footprint helped make this possible.
While these new results on DEM creation from Cryosat L1b SARIn data are encouraging, further work is required to better refine and exploit this potentially valuable data source.In the current work, airborne laser data was used to refine results, this is clearly not generally possible.However, other approaches are possible; low resolution DEMs for Antarctica or Greenland (e.g.Bamber et al., 2009) could be used as an initial input for simulation and then refinement of swath processing results.Also for appropriate cross-track slopes the POCA estimates from either the ESA L2 or L1b files could be used as reference data to help refine the baseline roll angle or the potential contribution from the range-ambiguous zone in swath processing.Work is underway to explore these options.The ultimate height accuracy possible with swath mode Cryosat processing remains to be determined but the Devon results imply that height errors with a standard deviation of less than 1m are possible, although there can be bias errors of order 1 m.

Figure 1 .
Figure1.Cryosat is oriented such that the baseline formed by the 2 receive antennas is oriented perpendicular to both the along-track direction and to the normal to the WGS84 ellipsoid.Each cross-track scan line will contain a unique estimate of the 'point-of-closestapproach (POCA) on the ice surface, followed in delay time by the sum of returns from both sides of the POCA.With a suitable geometry the returns from the ambiguous region can be much weaker than those from the main beam underneath the satellite which allows the mapping of position and height of main beam footprints.The dashed magenta line illustrates Fig.1.Cryosat is oriented such that the baseline formed by the two receiving antennas is oriented perpendicular to both the along-track direction and to the normal to the WGS84 ellipsoid.Each crosstrack scan line will contain a unique estimate of the POCA on the ice surface, followed in delay time by the sum of returns from both sides of the POCA.With a suitable geometry the returns from the ambiguous region can be much weaker than those from the main beam underneath the satellite which allows the mapping of position and height of main-beam footprints.The dashed magenta line illustrates diagrammatically the dB variation in cross-track two-way gain.Note that the angles have been exaggerated for clarity; the first side lobe in the pattern occurs at ∼ 2 • from nadir and represents a two-way power level of −40 dB with respect to the power at nadir.
Figure 2 (top): The surface area illuminated by a burst of 64 pulses extends fore-aft and sideto-side of the sub-satellite point.A locally circular Earth model has been used to illustrate the relative positions of some of the circular iso-range contours and half of the 32 fore Doppler beam positions.The central (dark blue) semi-circular contours illustrate the first 9 iso-range contours after the POCA.Subsequently, the blue semi-circles indicate the positions of every 25 th iso-range contour.The 512 samples in each waveform or scan line are separated in the L1b files used in this work by 0.47 m in slant range.The cross-track footprint dimension change rapidly from the large size at POCA (~1.5 km) to ~20 m for the 400 th footprint beyond the POCA.The positions of the Doppler beams extending across-track are shown by dotted magenta lines, for clarity every second beam is shown.Both distance and angular coordinates in degrees are given for the across-and along-track dimensions with respect to the POCA.

Fig. 2 .
Fig. 2. (Top): the surface area illuminated by a burst of 64 pulses extends fore-aft and side-to-side of the sub-satellite point.A locally circular Earth model has been used to illustrate the relative positions of some of the circular iso-range contours and half of the 32 fore Doppler beam positions.The central (dark blue) semicircular contours illustrate the first nine iso-range contours after the POCA.Subsequently, the blue semi-circles indicate the positions of every 25th iso-range contour.The 512 samples in each waveform or scan line are separated in the L1b files used in this work by 0.47 m in slant range.The cross-track footprint dimension change rapidly from the large size at POCA (∼ 1.5 km) to ∼ 20 m for the 400th footprint beyond the POCA.The positions of the Doppler beams extending across-track are shown by dotted magenta lines, for clarity every second beam is shown.Both distance and angular coordinates in degrees are given for the across-and along-track dimensions with respect to the POCA.(Bottom): the relative two-way antenna gain is shown for cross-track cuts at five different fore angles corresponding to the central and the 8th, 16th, 24th and 32nd beams.These are the dashed lines while the solid blue line corresponds to the effective two-way pattern G 2-way (θ) given in Eq. (6).Note that the cross-track angle extent in Fig. 2 (bottom) is much wider than that in Fig. 2 (top).

Figure 3 .Fig. 3 .
Figure 3. Ratio of the interferogram power in the ambiguous beam to that in the main beam for cross-track slopes between 0.5° and 2.0° based on a locally spherical Earth model.The minimum in the suppression of the ambiguous power occurs when the ambiguous footprints are illuminated by the low antenna gain region between the main lobe and the first side-lobe at ~ 1.5° from nadir.

Figure 4 .Fig. 4 .
Figure 4.The western part of the Devon Ice Cap with topography illustrated in colour.The sub-satellite position of the Feb. 24 th descending pass is shown as a black line and the positions of the airborne laser scanning altimeter data is shown in white.The dotted black lines indicate the sub-satellite tracks of the 25 descending passes which were used to create the DEM illustrated in Fig. 8.The position of the Devon Ice Cap in the Canadian Arctic is shown in the insert.In this figure the topography is derived from the Canadian digital elevation database merged with the CS swath processed results shown in Figure 8.

Figure 5 .Fig. 5 .
Figure 5. Illustration of the return power for the central part of the Feb. 23 pass (top), the filtered coherence (middle), and the differential phase (bottom), all as a function of the delay window expressed as a slant range distance.Note that the phase has not been unwrapped at

Figure 7 .Fig. 7 .
Figure 7. Illustration of the swath processed geocoded elevations for the Feb. 23 rd descending pass.The oblique white line indicates the sub-satellite track and the position of the maximum antenna gain.The 2 white arrows indicate the position of the E-W pass of the airborne ALS data which were used as a height reference.

Figure 8 .
Figure 8. Results from 25 descending passes have been combined to create a DEM of the western slopes of the Devon Ice Cap.Colour has been used to indicate height in meters.The uneven edges of the DEM at the north and south edges are due to the relatively large alongtrack slope and the resulting poor solution often indicated by lower data coherence and higher phase noise.

Fig. 8 .
Fig. 8. Results from 25 descending passes have been combined to create a DEM of the western slopes of the Devon Ice Cap.Colour has been used to indicate height in metres.The uneven edges of the DEM at the north and south edges are due to the relatively large along-track slope and the resulting poor solution often indicated by lower data coherence and higher phase noise.

Figure 9 .Fig. 9 .
Figure 9.The upper plot illustrates the E-W elevation at 73.339 N for both the Cryosat DEM data and the ALS laser heights.The middle plot illustrates the histogram of the height differences and the lower plot shows the E-W slope derived from the Cryosat results.

Figure 10 .
Figure 10.Comparison of the relative received power (upper image) with the simulated power for the main beam (lower image) for the Feb 23 descending pass data.In both cases the dynamic range from blue to red is 20 dB and in each case the dB scale is with respect to the maximum value.The simulation is based on the footprint size, a near-nadir backscatter model and the two-way antenna gain pattern given by Eqn. 6.

Fig. 10 .
Fig. 10.Comparison of the relative received power (upper image)with the simulated power for the main beam (lower image) for the 23 February descending pass data.In both cases the dynamic range from blue to red is 20 dB and in each case the dB scale is with respect to the maximum value.The simulation is based on the footprint size, a near-nadir backscatter model and the two-way antenna gain pattern given by Eq. (6).