Systematic Errors Observed in CryoSat-2 Elevation Swaths on Mountain Glaciers and Their Implications

Our awareness of ice caps’ and mountain glaciers’ sensitivity to climate change has driven major advances in the application of remote sensing techniques during the past decade. Regarding ESA’s SARIn altimeter CryoSat-2, processing the full waveform to generate swaths of elevation estimates has become standard practice in regions of complex topographies. This technique provides information on areas where we would be blind otherwise. In this article, we discuss systematic errors and analyze their impact on surface elevation measurements and change rates of two test areas. In particular, we focus on periodically occurring errors in elevation swaths, caused by the superposition of coherent signals from range-ambiguous surfaces. They can lead to measurement errors in excess of 10 m, affect most measurements in mountainous regions, are difficult to exclude with established post-processing techniques, and occur repeatedly for satellite revisits introducing a 369-day periodicity—difficult to distinguish from the annual cycle. We show a correlation between derived elevation swaths and the sensor view angle and explore the influence of common data exclusion choices on higher level products. Our results indicate that these systematic errors hold a substantial share of the error budget and that the choice of thresholds impacts higher level products. We conclude that error correlations need to be considered to characterize the data accuracy. With the established data editing strategies, systematic errors prevent resolving seasonal mass changes of single mountain glacier basins and impact aggregates over larger areas or longer periods.

munities in several mountainous regions will likely be negatively impacted by retreating glaciers, as they depend on the meltwater of glaciers, e.g., for drinking water or irrigation, or are threatened by landslides and floods [1]. Also, mountain glaciers and ice caps contribute substantially to sea level rise and have collectively been the largest contributor in the 20th-century [2]. To gain insights into the dynamic response of glaciers to changing environmental conditions, we depend on observational data with a high spatial resolution, fine enough to resolve the small mountain glaciers (kilometer scale), and a high temporal resolution, frequent enough to resolve seasonal changes.
CryoSat-2 is a promising tool for this task, given it provides interferometric synthetic aperture radar altimetry (SARIn) observations over glacier regions. While the traditional point of closest approach (POCA) method derives only a single high-precision surface elevation estimate from each waveform, i.e., the measured return power over time, the increasingly popular swath processing approach attempts to reconstruct all dominant echo origins that contributed to the signal. Calibration and validation studies, e.g., [3], [4], [5], tested the accuracy and precision of both strategies over the Devon ice cap, the Jakobshavn Glacier, and the Austfonna ice cap. Reported values for the swath elevation precision roughly range from 1 to 2 m [3], [5]. However, these results apply to the relatively smooth topographies of the studied test regions and may not be valid for mountain glaciers because of the complex interaction between the topography and the waveform.
For monitoring glaciers in the proximity of mountains, the POCA method leaves much of the glacier not surveyed making swath processing necessary. Therefore, swath processing has become the standard over ice sheet margins and ice caps because of the benefits in terms of spatial and temporal coverage [6], [7], [8], [9], [10], [11]. Separately, high temporal and spatial resolutions were already successfully obtained for several glacier regions. For example, Foresta et al. [7] and Tepes et al. [11] calculated long-term elevation change rates at a 500 m resolution and Morris et al. [8] and Jakob et al. [9] derived monthly elevation and mass change time series on regional scales. These studies show that enhancing the data quantity, even if it comes at the cost of precision, makes new insights possible. Simultaneously, obtaining a high temporal and high spatial resolution poses the next challenge in CryoSat-2 data usage. This would only be a question of data quantity if the elevation estimates This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ Fig. 1. Schematic of errors in CryoSat-2 swath elevation estimates. The white shape sketches a transect of a glacier in across-track direction, the dots (black, gray, or orange) indicate swath data, and the orange arc is equidistant to CryoSat-2. This range interval intersects the glacier surface at multiple places, such that the recorded signal is the sum of all reflections. The orange points correspond to this range interval. While the gray points overestimate the elevation, we will not focus on their cause in this study. Therefore, we have shaded them to distinguish them from the black and orange points.
were statistically distributed around the true surface and no systematic errors were present. Some previous studies [10], [12] acknowledge that a bias could be introduced by a change of the scattering horizon, i.e., the mean signal penetration depth, but explicitly report not to observe such a bias. Morris et al. [5] later show that for the glaciers on Svalbard this effect is observable and suggests that it could even be used to derive the thickness of the winter snowpack.
In this work, we investigate the impact of systematic errors on swath elevation estimates in mountainous regions that hamper advancements in increasing the spatiotemporal resolution. We focus on errors introduced by "range-ambiguous regions/footprints" as reported earlier by Gray et al. [6]. As we will show, these errors create an artificial signal with 369-day periodicity. In the presence of summits and valleys, errors on the order of 10 m occur. Fig. 1 illustrates a situation where this is the case.
Common data editing approaches that limit the impact of systematic errors exclude incoherent or weak signals to improve the data quality. However, this reduces the data quantity at the same time and affect higher level products such as surface elevation change rates, making the data selection a complicated trade-off. Improving our understanding of errors is essential to assess the significance of observations and knowledge of the origins of the errors could help to build effective algorithms to exclude them. Awareness of errors is, also, becoming more important with the growing use of CryoSat-2 swath data by non-radar altimetry experts. Not knowing about these errors, data users run the risk of misinterpretation.
In this article, we discuss and analyze the cause of one specific error source and its implications, alongside a study on how common data editing influences derived data products in the context of this kind of errors. The region shown in Fig. 2 serves us as a main testing and demonstration area, using a second area (see supplement Section S3.3) to ensure that our conclusions are not unique to the testing site. Our work contributes to the understanding of the challenges associated with using CryoSat-2 swath data in mountainous regions with a high spatiotemporal resolution. We will proceed with a Location and topography according to ArcticDEM of the test region in the spatial reference frame with ID 3413. The legend entry "Doppler cell" refers to Exp. A and is described in Section IV. (Left) Contains modified Sentinel-2 L1C imagery from August 23, 2017, provided by https://apps.sentinel-hub.com/eo-browser/EO Browser, Sinergise Ltd. summary of swath processing. We then discuss systematic errors focusing on range-ambiguous footprints. Next, we analyze conducted experiments to demonstrate the impact of the systematic errors. Finally, we discuss our findings.

II. SWATH PROCESSING CRYOSAT-2 DATA
This section briefly introduces the most important facts, concepts, and procedures related to the swath elevation estimates that are necessary to explain the origins of the systematic errors discussed in detail in Section III. Since we derive the elevation data ourselves from the L1b Baseline D data [13] provided by the European Space Agency (ESA), we describe our retrieval algorithm alongside. The supplement explains swath processing in further detail. We refer to ESA's documentation of CryoSat-2 [14] for details on the sensor, its measurements, and processing beyond the scope of the article and its supplement.
CryoSat-2 is a spaceborne radar altimeter [14], with an orbit inclined by 92 • , leaving only small data gaps at the poles, and with a 369-day repeat cycle and 30-day and 85-day subcycles. It sends Ku-band radar signals toward Earth's surface and records their echoes. We can derive the range from the sensor to the surface from the delay between sending and receiving the signals. For details about deriving the range we refer the interested reader to [14].
CryoSat-2 can operate in three modes: the Low-Resolution Mode, mainly active over the ice sheets' interior and the oceans, the SAR mode for sea ice and ocean currents, and the SARIn mode to monitor the ice sheets' peripheries and mountain glaciers. We use SARIn data and focus on this mode in the following. CryoSat-2 transmits coherent pulse bursts in the SARIn mode, with each burst consisting of 64 pulses. The pulse repetition frequency is about 18 kHz. Within a pulse, the signal frequency is modulated to obtain a chirp that can be decomposed into (including oversampling) 1024 measurements later in the process. After the echoes are recorded, they are split up into parts, according to the Doppler shift introduced by the relative movement of CryoSat-2 and the echo origins. Using this advantage, the post-processing can focus on specific echo origins in along-track direction. ESA resolves Doppler sharpened beams with a width of about 400 m. Each segment on the ground, called Doppler cell, is multi-looked: it is illuminated 64 times by bursts. The effective number of looks, however, is between 50 and 60 for the waveforms evaluated in the experiments section. In the L1b product, the along-track nadir locations are given and the multi-looks are averaged to a single waveform. Each waveform consists of power measurements for 1024 range bins. We refer to the 23.42 cm wide range intervals, corresponding to the range bins, as range slices.
CryoSat-2 carries two antennas, mounted in across-track direction, that enable interferometric measurements. From the records, ESA calculates the phase difference of the signal between the two antennas and its coherence for their L1b product. The coherence gives the fraction of the signal that is in-phase.
Conceptually, the echo's origin in the across-track direction can be calculated by comparing the times at which the echo arrives at CryoSat-2's left and right antennas. This time difference is, however, recovered from the phase difference of the signals, which is wrapped to the range −π to π radians. This introduces an ambiguity which results in a number of physically possible angles of arrival, each of which will lead to a different location. Two to three of those lie within the 3 dB (half-power) antenna footprint. We consider the 7 most central locations to be "physically possible" and calculate the corresponding elevations of those. This upper-limit footprint has a beamwidth of 7.6 • and measures about 95 km in diameter on the ground.
To decide which location is the most likely, we adopt the strategy reported by Garcia-Mondéjar et al. [15]. This involves treating consecutive waveform samples together as groups if their coherence exceeds a threshold. Additionally to their threshold criterion, we use restrictions on the number of grouped samples, treat samples with a deviating phase difference individually, separate sample groups at rapid phase difference changes, and apply an extra step to check the groups for ambiguities. We report details on the algorithm in the supplement. We use ArcticDEM v3.0 [16], a photogrammetrybased digital elevation model, to resolve the location ambiguity by finding the best agreement for each group or individual sample.
In the processing step from L1b data to elevation estimates, some data are commonly excluded by a coherence threshold and sometimes an additional power threshold. A low coherence indicates that the phase difference, and therefore the echo origin, is not reliable. The power threshold removes samples without any signal and those where the echo origin is far outside the antenna footprint. The power measurement is usually not given in Watts but in decibel-compared to a reference, for which we use 1 W. We do not subtract the noise floor or the incoherent part of the signal (see supplement Section S2.4 for an explicit formula). The thresholds influence the frequency and severity of the studied swath issues directly and play a major role in the experiments. We cling to common choices for thresholds, stated in Table I. We disregard elevation estimates with coherence values below 0.8 or with power values below ten times the median noise power of the satellite track if we do not explicitly state different values. Beyond that, we also disregard data that differ in their elevation more than 50 m from ArcticDEM or show a "suspicious" coherence: concretely, suspicious coherence values, that we exclude, are equal to 1 while the power is less than −155 dB. Because the coherence is physically constrained to be ≤(1/(1 + 1/SNR)) [19], with SNR being the signal-to-noise ratio, the coherence of weak signals cannot be high. To validate our retrieval algorithm, we compared a set of 58 000 elevation estimates to data from the https://cryotempoeolis.org/elevation over land ice from swath (EOLIS) product provided by ESA's CryoTEMPO project. For 79% of the EOLIS data, we find a corresponding point in the derived product deviating less than 10 m in across-track and less than 5 cm in the vertical. The elevation estimates for which we do not find a coinciding counterpart are differently located when resolving the location ambiguity. The supplement Section S2.6 gives further information on the comparison.

III. SYSTEMATIC ERRORS FROM RANGE-AMBIGUOUS SURFACES
CryoSat-2 swath data are subject to a number of systematic error sources. These include volume scattering, surface roughness, and the surface slope. The impact of volume scattering has recently been studied by Morris et al. [5]. They find a −1 to −1.5 m bias for swath data for the Austfonna ice cap. Wang et al. [20] analyzed the influence of slope, roughness, and volume scattering on the waveform.
When comparing CryoSat-2 swath elevation data to Arc-ticDEM, we observe a typical pattern for the majority of the waveforms. The data from the leading edge, i.e., the first part of a waveform before the first peak, tend to show an overestimation of the surface elevation, which gradually decreases and eventually changes to an underestimation in Fig. 3. Differences between swath elevation estimates and ArcticDEM. In all waveforms, a typical pattern with increasing sample number of a positive deviation from the DEM, followed by agreement between swath elevations and ArcticDEM, and finally negative deviation is present. Consecutive waveforms tend to show similar patterns. We chose the extent such, that the pattern is simply identifiable. The supplement Fig. S10 shows the satellite track extent with its surroundings. later samples. Fig. 3 shows a set of example waveforms. This pattern may occur repeatedly within a waveform. In this example, the pattern is easily visible, because the topography is similar in the consecutive footprints corresponding to about 15 km along-track. In the case of irregular topography, the pattern is still visible, but less apparent because it occurs at different sample number ranges. The systematic overestimation associated with the leading edge may have different causes, e.g., the slope, the surface roughness, and the penetration depth [20]. Here, we investigate the effect that contributes to the described underestimation and mainly influences the waveform samples after the first waveform peak and, thus, after the POCA which is usually assumed to be at or before the first peak.
Beyond the POCA, any range slice intersects at least two locations on the surface, such that reflections from those surfaces contribute to the same waveform sample [6], schematically shown in Fig. 1. If none of the reflections is stronger by multiple orders of magnitude, coherent reflections lead to measuring an incorrect phase difference leading to an incorrect geolocation [6]. For complex topographies, we find these incorrect geolocations to translate to elevation errors on the order of 10 m.
For topographic peaks, we observe signals which seem to originate below the actual surface. This is the case for range slices that cut through a peak and intersect with the surfaces at the left-and right-hand sides. The sum of the signals appears to originate at the common center of the scatterer distributions. The same is true for topographic depressions, except that we observe echoes seemingly above the surface. Fig. 1 highlights an example for the valley case-but it also shows returns "from inside" a summit on the right side of the sketch. The signals deviate several meters from ArcticDEM. The affected signals are marked by a drop of the coherence (between about 0.8 and 0.95, cmp. Fig. 4).
The further the surfaces are separated, the lower the coherence of their reflections is at the satellite. For CryoSat-2, the recorded signal becomes completely incoherent if their distance approaches ≈10 km [6], [21]. The decorrelation distance depends on the sensor type; it is proportional to the wavelength and the reciprocal antenna separation. A threshold . As a reference it shows the cross section of the glacier according to ArcticDEM (black line). The remaining are three stacked plots that show the waveform together with the coherence and the phase difference. Fig. 2 includes the location where CryoSat-2 recorded the signal labeled as "Doppler cell." on the coherence values can exclude a part of the ambiguous signals, but it comes at the cost of the quantity and the spatial coverage of the elevation estimates. It is difficult, and maybe impossible, to exclude all affected signals because of the nature of the problem: beyond the POCA, range slices intersect the surface in at least two places. In the best case, one of the echoes dominates the recorded signal (cmp. [6] with the illustration in their Fig. 1). The less powerful echoes cause a small drop in the coherence. If we considered only signals with the physically achievable maximum coherence, we would need to discard close to all measurements. However, the small coherence drop in turn only leads to a small deviation of the measurement location from the true surface. In most cases, this is not a problem and the measurements are still valuable. However, the impact of this error source on higher level products will be sensitive to coherence-based data editing.
Errors from the signal superposition of range-ambiguous surfaces depend on the position of the satellite because it defines the location, size, and beam incidence angle of those surfaces.

IV. EXPERIMENTAL SETUP
To study the effects of systematic errors and especially those of range-ambiguous footprints on swath elevation estimates, we conducted five experiments, named A-E. They will zoomed-out view from a detailed view, which shows the reproducibility of errors in contrast to the influence of CryoSat-2's position, to why those errors are problematic on a larger scale, by looking at interannual elevation change rates of a region of 165 km 2 . For the experiments, we derive elevation estimates from ESA's L1b product, as described in Section II, over a test region at the margins of the Greenland ice sheet. Fig. 2 shows the extent of the test region. This region is part of a glacier and does not include mountaintops or any parts not covered by ice (see supplement Fig. S11). We construct the polygon alongside topographic features, constrained by mountain ranges to the east and west, by the 2350 m elevation line in the north, and by the locally highest elevations in the south. The demonstration site is entirely located in the accumulation zone. According to the firn densification model IMAU-FDM [22], the region undergoes only minimal intra-annual variations in elevation due to surface processes.
We show the generality of our results by applying the same methodology to a second test region with different characteristics, the Glazov glacier. The Glazov glacier is located on Novaya Zemlya. It covers about 400 km 2 on elevations between 0 and 1050 m and is marine-terminating [23]. While we report some of the results in the following for comparison, please refer to the supplement Section S3.3 for figures and further details.

A. Case Study of Errors' Dependence on Sensor Position
First in Exp. A, we conduct a case study to demonstrate that the errors from the example used in Section III are reproducible and that the sensor position determines the errors. Here, we compare measurements at the same location as already shown in Fig. 4 with the following other recordings.
1) Two exact orbit repeats, temporally separated by 4 and 8 years. 2) Three approximate orbit repeats, temporally separated by multiples of 369 days. 3) One 85-day subcycle track, separated by 3.8 km in across-track direction. 4) The 30-day subcycle tracks of 2019. The supplement Table S1 identifies the recordings in detail and Fig. 2 shows the location of the Doppler cell. We compensate elevation estimates for differences in the along-track satellite position, which are smaller than 150 m in all cases, by adding the elevation difference according to ArcticDEM.

B. Errors' Effect on the Leading Edge and on Monthly Time Series
In the next experiments, we analyze how systematic errors show throughout the waveform and in the average surface elevation of the demonstration site. In Exp. B, we use all elevation estimates to compare the deviation from ArcticDEM to the relative sample number. We define this quantity as the waveform sample number after subtracting the sample number of the first peak of the smoothed waveform (see supplement Section S3.1).
In the next step, we produce and analyze a time series of monthly average elevation differences to ArcticDEM to investigate the temporal systematic behavior.

C. Perspective Versus Temporal Correlation
In Exp. C, we investigate the impact of the perspective error correlation on a statistical basis. We compute the correlation of elevation estimates between pairs of either ascending or descending tracks and bin those according to their across-track distance and their temporal separation. In detail, we average differences between swath elevation estimates and ArcticDEM over 500 by 500 m grid cells, separately for each track. For each pair of tracks that covers 100 or more common grid cells (≥25 km 2 ), we compute the correlation.

D. Time Series' Dependence on Thresholds
In Exp. D and E, we turn to the impact of systematic errors on products of longer aggregation periods. We compare data above different coherence and power thresholds, as commonly used for data editing in the community (see Table I). At least range-ambiguous footprints have a fingerprint in the coherence (see Section III), but other potential error sources with correlations to coherence or power contribute to the results.
In Exp. D, we calculate a 3-monthly time series of the average surface elevation over the test region, comparable to, e.g., [7], [9], [10] and roughly corresponding to the 85-day subcycle, for a variety of thresholds in the range used in previous studies (see Table I). In detail, we processed the CryoSat-2 tracks using a coherence threshold of 0.6. We, then, constrain the obtained elevation estimates further to those above certain thresholds for the signal power and the coherence. As power thresholds, we use the values −165, −160, −155, and −150 dB and as coherence thresholds 0.6, 0.7, 0.8, and 0.9. We calculate time series for all 16 combinations of power and coherence thresholds. For the time series, we calculate the average surface elevation differences to ArcticDEM for each time step. For this, we assume that within this 25 by 7 km region, surface elevation estimates at similar elevations are similarly affected by the glacier dynamics, the surface mass balance, the firn density profile, and the surface roughness. This assumption allows us to aggregate data over elevation bands. We choose 50 m elevation bands for this task because they feature a reasonable number of CryoSat-2 measurements per band as well as a sufficient resolution to account for elevation-dependent processes, e.g., precipitation, melt, and firn compaction. The elevation difference of the entire test region is calculated from the average of all bands, weighted by their area. The supplement Fig. S13 contains histograms of data counts per band and per 3-month for each pair of thresholds.

E. Thresholds' Impact on Interannual Change Rates
In Exp. E, we test to what extent the threshold dependence persists for interannual elevation change rates. We proceed as in the previous experiment, but instead of building a time series, we calculate the differences of weighted average elevation differences aggregated over each glaciological year, starting September 1.

A. Case Study of Errors' Dependence on Sensor Position
Inspecting the swath elevation estimates in detail suggests a link between the agreement of the elevation estimates from different tracks and the inter-track spacing. Fig. 5 compares elevation estimates different overpasses to the December 6, 2019, track which is also shown in Fig. 4. Fig. 5(a) compares only the best matching orbit repeats. These are separated from each other by less than 200 m and occur roughly every 4 years. They are similarly affected by the systematic errors: at the right of the 0 km-mark, the pattern of overestimation is typical for data on the leading edge, while between elevations 2715 and 2730 m the deviations are in line with misplaced across-track echo origin due to range-ambiguous surfaces. Fig. 5(b) shows 4 approximate orbit repeats. These are separated by 369 days and up to 2.3 km. These elevation estimates show how the data location changes with the perspective of the satellite. The number of points that strongly deviate from the reference surface around the −2 and 0 km-mark is different for each track, as is the magnitude of the deviation. Those changes can plausibly be explained by the changing contribution of range-ambiguous surfaces to the measurements. The 2017record is taken from 600 m further to the right (considering the figure). This increases the range to the left summit, so that it contributes only to later waveform samples. Around 0 km, this record returns lower elevation estimates than the 2019record. Similarly but amplified, this holds for the 2018-and 2016-track that are even further right. In these records, the coherence does not drop below the threshold resulting in a comparably continuous trail of elevation estimates bridging the data gap of the 2019-and 2017-records.
The changes are even more pronounced in Fig. 5(c). Here, an 85-day subcycle is shown with an inter-track spacing of 3.8 km. Above 2730 m, the observations differ from each other with less pronounced positive deviations from Arctic-DEM in the February 28, 2020, data. Below that, the February data are located right of the December 6, 2019, data. Contrary to the December overpass, the February overpass yields data between the across-track marks −1 and 0 km. Further left, the two overpasses return different elevation estimatesdisagreeing by up to 10 m in the vertical. As above, we believe the observed changes result from changing contributions of range-ambiguous surfaces. There are three places in this plot where the February track returns lower surface elevation than the December track (0.3 km across-track, elevations above 2740 m and below 2730 m; −2.2 km across-track, elevations above 2710 m). A lowering of the scattering surface for the February track is not expected. The disagreement of the highest data concerns data on the leading edges of the waveforms. At the other locations, range-ambiguous surfaces are a likely error source.   Fig. 5(a) and (b), we conclude that the errors and the data coverage vary with the position of the satellite. Changes in the volume scattering and surface roughness on their own cannot plausibly cause the observed "surface evolution". For example, the changes at −0.5 km between February 25 and July 18 would require a melt layer that formed on top of light snow. This, however, seems unlikely because of the large negative difference to ArcticDEM at −3.5 km in July. To make sense of the observations, we need to consider the impact of range-ambiguous surfaces.

B. Errors' Effect on the Leading Edge and on Monthly Time Series
To analyze the relation between the deviation of swath elevation estimates from ArcticDEM to the relative sample number, we aggregate all data per sample number and calculate the median deviation. Fig. 6 visualizes the results. We find that waveform samples on the leading edge systematically result in higher elevation estimates than expected from ArcticDEM with about 12 m in average for relative sample numbers −100 to −80. The deviation from ArcticDEM, then, decreases for increasing relative sample numbers. Fitting a trend to the median deviations of relative sample numbers −20 to 0 results in an average rate of (−0.203 ± 0.001) m per sample and an offset of (0.33 ± 0.01) m at the first-peak sample. The seemingly linear relationship between the relative sample number and the deviation from ArcticDEM indicates that the data are not independent but carry mostly the same information. The deviation takes its lowest values at sample numbers 30 to 50 of about −3.6 m in average, before increasing again. Fig. 7 shows the resulting time series for monthly average elevation changes over the demonstration site in Fig. 7 (top). The time series shows irregular fluctuations, with no seasonal signal visible. Fig. 7 (middle), where the average elevation differences are shown with regard to the season, suggests that any potential seasonal signal would be hidden by the variance of the monthly means. An autocorrelation analysis, shown in Fig. 7 (bottom), reveals significant correlations for lags of 1, 3, 6, 9, and 12 months with coefficients of −0.22, 0.35, 0.21, 0.27, and 0.42 and p-values of 0.025, 0.0003, 0.033, 0.005, and 2×10 −5 , respectively. Recalling the position dependency found in the first experiment and comparing the autocorrelation coefficients, we find that the autocorrelation coefficients reflect rather the perspective correlation than the temporal one. We observe significant positive correlations every 3 months (3, 6, 9, and 12). The distance between "parallel" tracks separated by 85 days (about 3 months) is substantially smaller than that of the 30-day subcycles and leads to a higher perspective correlation [cmp. Fig. 5(c) and (d)]. Similar to finding a negative 1-month lag correlation, the contrast between the negative (−0.16) 11-and positive (0.42) 12-month lag correlations stresses that an annual cycle is likely not responsible for the signal. For any realistic annual cycle, the 11-month correlation coefficient should be comparable to the 12-month coefficient and should not be negative. We interpret that the difference in the perspective of CryoSat-2 (which roughly repeats yearly but is different otherwise) causes the unexpected contrast in the correlation coefficients. This consideration is important because when interpreting time series one might be biased toward assuming a temporal correlation while the perspective correlation might be more important. The periodicity of the orbit links the two correlations to each other such that it is difficult to assess them separately.
For errors from volume scattering an annual cycle would be expected [5]. These errors have not been reported to depend on the satellite perspective in the context of swath processing. Contrary, it is plausible that range-ambiguous footprints cause the observed perspective correlation and we discuss a potential contribution of data on the leading edge in Section VI. While we cannot quantify the error budget per source, especially the negative 1-month lag correlation suggests that perspective-dependent errors dominate.
For Glazov glacier, we conducted the same analysis and show the results in supplement Fig. S5. We find similar correlation coefficients. The ones significantly different from 0 are at 3, 5, 6, 8, and 12-month lags with values of 0.38, −0.19, 0.22, −0.20, and 0.34 and p-values of 1 × 10 −5 , 0.031, 0.012, 0.020, and 8 × 10 −5 , respectively. Despite the insignificant 9-month lag coefficient of 0.12, the 3-monthly scheme pointed out above is even more apparent in the mid and bottom panels of supplement Fig. S5.

C. Perspective Versus Temporal Correlation
In Fig. 8, we compare track pairs that return surface elevations for the same locations. It shows the correlation means within bins of distances and track delay times and their 95% confidence intervals as error bars. We find that the correlation decreases monotonously with increasing track spacing and stays below 0.2 for separations bigger than 8 km. The temporal correlation is the highest for the orbit repeat time and the second-highest for the 85-day subcycle (or 85 days before the orbit repeat). The pattern repeats for longer delays, visible in the supplement Fig. S12. The observable disagreement between the correlation of ascending and descending tracks highlights that the results are unique to a specific setting.
In general, the track separations and the delays occur repeatedly as tuples that belong together. Even though we show the two associated correlation assessments next to each other, they are entangled. Without a control dataset that resolves seasonal dynamics, it remains impossible to quantify how much of the correlation is caused by the errors that show in the perspective evaluation and how much is indeed a temporal correlation. Considering the contrast between the low correlations for the 30-day subcycle delay and the higher correlation for the 85-day delay, we conclude that for the test region the perspective correlation is more pronounced than the temporal correlation. We find similar results when analyzing the results for the second region, shown in the supplement Fig. S6.

D. Time Series' Dependence on Thresholds
The results in Fig. 9 show that the average surface elevations for our test region temporarily diverge for different power and coherence data filtering thresholds. Most prominently, the time series differ about 1 m around 2012/2013. How far they deviate from each other depends on the choice of offsets. We choose the offsets such that they minimize the deviations between the time series. For the third quarter of 2020 (the prominent spike in down-and upward direction), there is reduced data availability and even a lack of data points with signal power ≥ − 150 dB for some elevation bands. The features that all time series share are: 1) the slight decrease between 2011 and 2015; 2) the elevation gain around mid-2015; 3) the stable period between 2016 and late 2019; and 4) the decreasing trend afterward. The time series do not agree on annual patterns such as a seasonal cycle.
Annual patterns are, indeed, difficult to identify, except for the time series using signals with coherence ≥0.6 and power ≥ − 150 dB. In this case, the annual maximum surface elevation is between July and September of each year and the minimum is in February for most years. Unfortunately, this is difficult to verify by bare eye because of the large number of time series shown; the supplement Fig. S14 shows all time series separately. The IMAU-FDM does not show the observed fluctuations, but it agrees on the general trend.
To test whether there is an obvious impact of systematic errors, we exploit their dependence on the satellite perspective and position. We divide the ascending and the descending tracks from each other and find that the ascending tracks carry most of the "seasonal" fingerprint while the descending tracks show the features (1-4) more prominently. This is shown by the blue and red lines in Fig. 9 (bottom). The angle between ascending and descending tracks measures about 14 • at this latitude. The supplement Fig. S15 shows that the data coverage is similar for both sorts of tracks and does not explain the difference.
We conclude that systematic errors likely affect the ascending and the descending tracks differently and generate a part of the observed signal on the scale of 1 m with a yearly periodicity. This limits the detection of seasonal signals from CryoSat-2 elevations for regions of comparable (or smaller) size as our test region. Because of CryoSat-2's 369-day orbit repeat interval, errors depending on CryoSat-2's position can yield a high correlation with seasonal climatic conditions, e.g., surface melt, and lead to misinterpretation, if not identified correctly. One needs to stay alert for artificial and only seemingly seasonal signals in time series for regions with just a few tracks contributing data to each time step.
We observe that the time series only diverge temporarily. This is in line with the earlier finding (Section V-A) that the measurements are closely repeated every four orbit repeats. It could imply that for a 4-years aggregation period, there is little data editing dependence and that high spatial resolutions are possible. However, the surface elevations in the test region have been stable over the period of measurements (cmp. IMAU-FDM), and this conclusion may be invalid for dynamic regions where the shape of the topography changes substantially.
We repeated the analysis for the Glazov glacier and come to the same conclusions: an overall decreasing surface elevation, no annual patterns, and disagreement between the observations of ascending and descending tracks. For this glacier, the time series for the different threshold pairs deviate further from each other, e.g., about 10 m in an extreme case between October and December 2012. We show the time series in the supplement Fig. S7.

E. Thresholds' Impact on Interannual Change Rates
For most interannual change rates, the results for the different threshold pairs spread over more than 0.5 m yr −1 and the average standard deviation is 0.2 m yr −1 . Fig. 10 (left) gives an impression of the spread of change rates per period and Fig. 10 (right) contains an example of how the different results are distributed over the space spanned by the two thresholds. For most periods, the change rates include positive and negative values. This makes an interpretation of the regional changes difficult. For the positive change rates in 14/15-15/16, we suspect indeed climatological reasons. The IMAU-FDM shows above-normal accumulation in the winter 14/15 that is not entirely melted in the following summer (Fig. 9).
Turning to how the change rates are distributed over the threshold pairs, shown in Fig. 10 (right), we notice that the values are not randomly distributed but that there instead seem to be correlations. We recognize that varying the power threshold consistently covers a larger range of change rates than varying the coherence threshold. The change rates for the power threshold −150 dB are notably different from the rest (also supplement Fig. S16). We do not find any alerting or unsuspected shortcomings or steep decreases in the data availability between threshold values. Fig. 11 shows at the hand of the first neighboring pair of positive and negative change rates in Fig. 10    in both cases. We do not observe major differences between these periods. The data density is higher for the less strict threshold, especially in the south. In both cases, all bands are sufficiently covered with measurements. changes or sparsely sampled elevation bands that could cause the threshold-dependent observations. The strictest threshold cuts the available data to 10% compared to the laxest one, but still passes our threshold of minimum data counts (supplement Fig. S17). The differences actually result from scattered and locally confined changes in the estimated surface elevation (supplement Section S3.4).
We repeated the analysis for the Glazov glacier. Here, the interannual change rates spread in 9 of (here) 11 cases over more than 1.5 m and their average standard deviation is 0.8 m yr −1 reflecting the larger deviations observed in Exp. D. We show the results in the supplement Fig. S16.

VI. DISCUSSION
First, we discuss the implications for the goal of producing higher level data at a sufficient spatiotemporal resolution to resolve the elevation change signal of mountain glaciers. The results reported in Sections V-B and V-C suggest that a large part of the fluctuations in the time series is the fingerprint of the sensor perspective leading to systematic errors. Because of the almost yearly (369-day) orbit repeat, the errors are entangled with seasonal changes. The measurements, thus, can only be of any value, if we manage to correct or at least exclude the errors, or if we aggregate data from many tracks, such that the errors average out. Until now, to our knowledge, correcting or excluding erroneous data is not possible at the required level, leaving us with the "aggregation strategy".
Such an aggregation-based strategy has been applied in previous studies. Aggregating data over whole regions, e.g., Hindu Kush, Jakob et al. [9] produced elevation change time series using a 3-month moving average. Tepes et al. [10] aggregated all data over multiple years to estimate the elevation changes of 500 by 500 m grid cells. The "optimal" spatiotemporal resolution, i.e., the area and the period over which to aggregate data, for a specific area of interest depends on the local topography and the change rates.
We turn to discuss a potential further source of errors that occur periodically next. Considering any Doppler cell, what is a POCA to one overpass may not be a POCA to another overpass. However, in the POCA-case, the surface elevation is biased depending on the relative waveform number (see Section V-B). Since the same surfaces will be POCAs for orbit repeats every 369 days, such systematic errors would contribute to a periodic signal. To limit the contribution to a track-dependent bias, we advocate avoiding to use swath data around the POCA in mountainous regions.
In Exp. C, we observed differences between the perspective correlation for ascending and descending tracks. Whether different tracking points for ascending and descending tracks can explain this observation partly may be an interesting question to answer for future research.
Finally, we discuss our findings in light of performance analysis and uncertainty treatment in earlier research. In line with our findings, Andersen et al. [24] report standard deviations from 6 to 17 m for CryoSat-2 swath elevations compared to air-and spaceborne laser altimetry, and airborne X-band radar altimetry data. Previous studies, which investigated surface elevation change rates, usually did not consider the uncertainties of CryoSat-2 swath data directly. Instead, they calculated the uncertainties of the change rate products from the spread of the linear regression residuals [7], [8], [9], [10], [11], [12]. This approach assumes uncorrelated and normally distributed residuals, which is invalid as shown in this study. Morris et al. [8] assumed different tracks to be uncorrelated. Jakob et al. [9] weight CryoSat-2 elevation estimates based on the signal's coherence and power which can act as implicitly accounting for the correlation of the data by reducing the effective sample size. On the one hand, the relations between the coherence, the power, and the measurement correlation are not explained, on the other hand, in the light of our findings, it seems to carry potential as a pragmatic choice to account for the likelihood of data to be affected by the here presented errors. The other studies, [7], [10], [11], [12], do not mention how they handle the implicit assumptions in the uncertainties.

VII. CONCLUSION
Using swath processing for CryoSat-2 data provides a wealth of elevation estimates. The number of data points suggests that high spatial and high temporal resolutions are possible, despite the large errors. A substantial part of these deviations results, however, from correlated errors and forms an obstacle to achieving such resolutions.
We analyzed the impact of separated range-ambiguous surfaces on CryoSat-2 swath elevation estimates in mountainous regions. Signals of which the origins are separated on the order of 1 km in the across-track direction interfere. Independent of the retrieval algorithm, this interference leads to an inaccurate geolocation. The measured coherence of affected signals is reduced and, in theory, a coherence threshold can be used to exclude affected data. However, excluding data does not improve the product quality, if the removed data still carries information. We showed that applying coherence and power thresholds influences the surface elevation calculated from many swath elevation estimates.
A share of the errors on the order of 10 m is systematic with respect to the sensor perspective, i.e., the relative position and orientation of CryoSat-2 toward the topography. Because of CryoSat-2's 369-day repeat orbit, these errors can lead to a signal with yearly periodicity in derived surface elevations. This needs to be considered if one studies correlations with seasonal climate patterns. They should in theory affect SARIn altimeter measurements in general, but their impact depends on the SAR focusing, the radar band, and the antenna distance. Future SARIn altimetry missions could consider different orbit repeat periods to alleviate the entanglement between errors and the annual cycle. However, for multipurpose missions such as CRISTAL [25] the benefits of yearly orbit repeats for ice sheets or oceans, i.e., comparability, may be more important.
Although we focused on errors related to topographic undulations in the across-track direction, we note that in the presence of along-track slopes on the order of 1% or steeper additional systematic errors may be introduced. In such cases, the extent of the along-track footprint could lead to a dispersion of elevation estimates on slopes. However, the returns from the footprint margins will be weaker and can easily be hidden by competing returns. Because there are no competing returns above the POCA, such effect could bias measurement aggregates toward overestimation. A proper analysis of the hypothesized phenomenon is beyond the scope of this study.
We conclude that over mountain glaciers, where because of the spatial variability and the fast response to climate change a high spatiotemporal resolution would be extremely valuable, the current processing strategies do not allow to derive high-resolution maps of seasonal elevation changes. A possible strategy might be to develop analytical models describing the error correlations in space and time using reference data and use these to build the full noise variance-covariance matrix needed in estimating the elevation change rates. Implementing such an approach to obtain estimates over large areas requires, however, dedicated algorithms and sufficient processing capacity to handle the vast amounts of data obtained in swath processing.