**2. Materials and Methods**

#### *2.1. SAR Offset Tracking*

The calculation of displacement fields using Synthetic Aperture Radar (SAR) data is possible using offset tracking methods. Offset tracking is normally used for analyzing rapidly moving objects, such as glaciers, because it overcomes the loss of signal, due to decorrelation, when using InSAR techniques over long time intervals [18–21]. As the glaciers in the Cordillera Blanca of Peru are relatively small, offset tracking is only feasible using high-resolution SAR images with a resolution of ~2–3 m. Using a unique, continuous time series of archived TerraSAR-X (TSX) radar imagery for the period 15 January 2017–6 April 2020 (3.2 years) we reconstructed ice surface velocities for the Palcaraju glacier. The dataset consists of 51 TSX radar images that were acquired on a descending orbit with a ground resolution of 3 m (Table 1). The radar images were acquired at 10.53 GMT and the distribution of re-visit intervals for the 50 image pairs was 11 days (2 pairs), 22 days (43 pairs), 33 days (1 pair), and 44 days (4 pairs).

**Table 1.** TSX SAR offset tracking displacement maps time intervals, perpendicular baseline, and resolution.



**Table 1.** *Cont*.

Utilizing a normalized cross-correlation of the amplitude component of the SAR images, offsets were measured with rectangular windows at a set of positions uniformly distributed over the scene. To obtain an accurate estimate of subpixel precision of the correlation peak, correlation function values were fitted using a biquadratic polynomial surface. The time interval of the image pairs can be adjusted according to the expected maximum displacement over the glaciers from 11 days to several months. Mismatches or errors were filtered by applying a threshold to the correlation coefficient [22], by (1) iteratively discarding spurious matches based on the angle and size of displacement vectors in the surrounding areas, (2) applying a low-pass filter on the resultant displacement fields, and (3) applying a 2–98% cumulative cutoff to remove potentially uncertain velocity values and outliers. Such values cannot be realistically validated and are likely to be artifacts of the radar data processing method. Hence they have no physical meaning.

Slant range and azimuth offset displacement fields were geocoded and transformed to 3D displacements along the terrain surface using the Advanced Land Observing Satellite (ALOS) World 3D (AW3D30), Digital Elevation Model (DEM) [23]. TSX images were processed in series with offset-tracking procedures [18–21] to 3D ice surface displacement maps. This involved combining slant-range and azimuth offsets by assuming that flow occurs parallel to the ice surface, as estimated from the DEM. Matching window sizes of 128 × 96 pixels (e.g., 202 m × 192 m) were applied with steps of 16 × 12 pixels (e.g., 25 m × 24 m). The displacement maps in m/a were geocoded to a posting of ca. 60 m. For an estimation of the uncertainties in the ice displacement maps, a precision of 1/10th of a pixel in the offset estimation can be assumed [24,25]. The displacement error of TSX data with pixel sizes in ground-range and azimuth direction of 0.9 m × 2.0 m respectively, and a time interval of 11 days, are thus on the order of 10 m/year. For longer time intervals, the noise level increases, hence, a similar displacement error of about 10 m/year can be assumed for other image pairs. The displacement values were extracted from individual pixels from the spatial displacement maps.

#### *2.2. Satellite Image Analysis*

Surface features and glacier extent were mapped in a geographic information software system [26] using high-resolution satellite Pléiades imagery (Table 2). In addition, we visually inspected Google Earth and ESRI World Atlas imagery (Table 2) captured between 1999 and 2020. Aspect and slope were calculated using the AW3D30 DEM, and, together with velocity and elevation data, exported with one data point per pixel for further statistical analyses.


**Table 2.** Satellite imagery analyzed in this study.

\* as published on https:\livingatlas.argis.com/wayback, accesed on 8 July 2021.

### *2.3. Statistical Analysis of Time Series*

The residuals from a linear fit of the glacier surface velocities over the 3.2 year observation period were used to conduct various correlation analyses. The analyses were based on about 1300 measurement points. Data time series with less than 50% completeness were removed from the statistical analysis. For the remaining data, where gaps were evident, a linear interpolation was applied to ensure consistent time intervals. The filtered and interpolated data allowed us to perform spatial and temporal correlation analysis between the residuals, precipitation, and SST. The Pearson correlation coefficient was determined for all correlation analyses.

As an independent test, the processed time series were also subjected to a time series analysis, without a prior assumption of correlations with climate or weather data. The entire data processing workflow was performed in Python using standard packages for transformation and optimization. The dominant wavelength was extracted for each time series using a Fourier transformation. From this dominant signal, a fit to the time series was performed with a Levenberg-Marquardt or dogleg algorithm, resulting in a single fitted trigonometric function with an estimate of best-fitting values for amplitude, phase, offset, and period, as well as the covariance matrix. Physically unreasonable outliers were removed based on the following criteria:


In total, 28 data points were excluded based on the above criteria.
