Supplement of Changes in ground deformation prior to and following a large urban landslide in La Paz , Bolivia , revealed by advanced InSAR

Abstract. We characterize and compare creep preceding and following the 2011 Pampahasi landslide (∼ 40 Mm3 ± 50 %) in the city of La Paz, Bolivia, using spaceborne RADAR interferometry (InSAR) that combines displacement records from both distributed and point scatterers. The failure remobilised deposits of an ancient landslide in weakly cemented, predominantly fine-grained sediments and affected ∼ 1.5 km2 of suburban development. During the 30 months preceding failure, about half of the toe area was creeping at 3–8 cm/a and localized parts of the scarp area showed displacements of up to 14 cm/a. Changes in deformation in the 10 months following the landslide are contrary to the common assumption that stress released during a discrete failure increases stability. During that period, most of the landslide toe and areas near the headscarp accelerated, respectively, to 4–14 and 14 cm/a. The extent of deformation increased to cover most, or probably all, of the 2011 landslide as well as adjacent parts of the slope and plateau above. The InSAR-measured displacement patterns – supplemented by field observations and by optical satellite images – indicate that kinematically complex, steady-state creep along pre-existing sliding surfaces temporarily accelerated in response to heavy rainfall, after which the slope quickly achieved a slightly faster and expanded steadily creeping state. This case study demonstrates that high-quality ground-surface motion fields derived using spaceborne InSAR can help to characterize creep mechanisms, quantify spatial and temporal patterns of slope activity, and identify isolated small-scale instabilities. Characterizing slope instability before, during, and after the 2011 Pampahasi landslide is particularly important for understanding landslide hazard in La Paz, half of which is underlain by similar, large paleolandslides.



Landslides in La Paz
The 2011 Pampahasi landslide is one of seven historical landslides exceeding 1 M m 3 in the La Paz area.Landslides possibly similar in size to the 2011 failure happened in 1582 and 1873 in the southwest part of the La Paz valley system, but little is known of these events aside fro m sparse written accounts.The first landslide affected an area of 2 km 2 or mo re and buried the villages of Cano ma and Ango Ango1 , which were likely located in Llojeta Valley (Fig. S1).Th is landslide claimed about 200 lives (Cabeza de Vaca, 1586; transcribed in Jiménez de la Espada, 1965, p. 342-351 andin Arispe, 2011).The second event involved extensive ground movement over an area o f about 8 km 2 (Markham, 1874) and more rapid, localized movement that destroyed an area named Temb laderani (Markham, 1874) located southwest of the city centre, causing 32 deaths (Crespo, 1902).Its location may correspond to the modern area o f the city ca lled Tembladeran i (Fig. S1), although Dobrovolny (1962) applies this name to a landslide deposit in Llo jeta Valley.The loss of life fro m these two events suggests velocities exceeding several metres per second or, at the very least, a lack o f pro mpt evacuation.Four smaller landslides happened in the twenty -first century in the Llo jeta and Allpaco ma valleys near the margins of the La Paz and Achocalla basins (Hermanns et al., 2012;Fig. 1B).They destroyed about 51 ho mes but did not claim any lives (Robert s, 2016) (Tab le S1).Most of these failures involved weakly lith ified, fine-grained sediments of the middle part of the La Paz Format ion, including large areas that had been previously mobilized by large paleolandslides.S1.See Figure 1B for context within the La Paz area.

Table S 2: Historic landslides exceeding 1 Mm 3 in the La Paz area
Note: Locations, where known, are shown in Figure S1. a No data denoted by '-'.b Mechanisms and rates after Hungr et al. (2014) are approximate and based on local topography, relations to older large landslide deposits, and field observations (Dobrovolny, 1962;Anzoleaga et al.,1977).Rate estimates apply to the main failure mechanism.c Total area affected by the landslide.Areas of pre-20th century landslides are based on written accounts and are approximate.Areas of 21st century landslides were accurately m easured using highresolution satellite imagery (see Roberts, 2016).d Order-of-magnitude volumes based on area and depth mentioned in written accounts or estimated in the field.
e Material type based on field observations (denoted by '*' in Source column) or on a comparison of the landslide location with the geologic map of Anzoleaga et al. (1977).Order of listed materials reflects their relative contribution to the overall landslide area, with the most extensive material types listed first.Material types involving less than 10% of the landslide area are not included.f So urces: 1, Cabeza de Vaca (1586); 2, Markham (1874);3, Crespo (1902);4, Hermanns et al. (2012);5, Roberts (2016);6, Constance (2005); 7, Mobarec et al. (2008);8, Quenta et al. (2007);9, Quenta et al. (2008);10, Aguilar (2013 2 Geomechanical properties of principle units Anzoleaga et al. (1977) prov ide detailed sedimentological and geotechnical characterizat ion of lithostratigraphic units throughout the La Paz area, including units underlying the Pampahasi slope (Table S2).Their data were derived from laboratory testing supplemented by field observations of structure and lithology.Most days during the week preceding the landslide had little or no precipitation (Fig. 4B).However, 25 February, the day before the 2011 Pampahasi landslide began, was one of the wet test (39.2mm) in the entire 72-year period of the record and was one of only 15 days with mo re than 35 mm of rain (Tab le S3).Cu mu lative seven -day rainfall peaked over the next few days (26 February to 1 March) due to continued, although lesser, precipita t ion; this period corresponds to the four-day period of the landslide.Thus, heavy rainfall on 25 February, in co mbination with precipitation during the preceding days, likely played a role in triggering the failure.
Despite the high variability of short-term precipitation evident fro m eyewitness accounts of rainfall intensity and fro m daily precipitation records, the longer term (monthly to annual) amount of precipitation before and after the 2011 landslide was similar (Fig. S2).Precipitation during the 10-month, post-failu re period (March to December, 2012) totalled 282.2 mm and ranged on a monthly basis from 0 to 127.6 mm.These conditions are typical of Marchto-December precip itation in each of the three years fully covered by the pre -failure InSAR stack (i.e. 2008, 2009, and 2010: 0 to 137 mm monthly; 235 to 366 mm 10-month total; 282 mm 10-month mean) (Table S4).   4 InSAR processing methodology

InSAR background
The ability of InSAR to measure sub-centimetre-scale ground displacements over large areas makes it a particu larly useful tool for characterizing and monitoring slow-moving landslides.In recent decades, InSAR has been increasingly used in landslide investigations due to improvements of RADAR systems and processing techniques that allow for better phase recovery and increased spatial resolution.
HDS-InSAR increases the spatial density of displacement histories relative to conventional InSAR techniques through two approaches (Eppler and Rabus, 2011;Rabus et al, 2012).First, it utilizes both coherent distributed scatterers (e.g.pavement, agricultural land, and natural slopes) and pers istent point scatterers (primarily built structures and fortuitously oriented rock exposures).Second, it applies adaptive filtering to preserve spatial resolution as much as possible while optimally suppressing noise fro m surface decorrelation.Adaptive filtering improves the spatial quality of the phase as well as RA DAR echo intensity.The HDS-InSAR technique also includes several enhanced steps to remove errors introduced by other phase components, particularly the error fro m atmospheric water vapor variat ions.The HDS-InSA R processing approach generates stronger statistics for stacks of approximately 15 or more scenes.

Datasets
We applied HDS-InSAR to Fine Beam mode RADA RSAT-2 scenes acquired fro m an ascending orbit (look direction slightly north of east) with an incidence angle of 36.3°fro m vert ical at the image center.The dataset comprises 44 scenes collected between September 2008 and December 2011 (Table S5).S6).

Processing summary
The processing chain that we used is an updated version of the prototype described by Rabus et al. (2012).It includes specific steps to remove differential phase contributions due to sensor squint (Doppler) variations, incomp lete knowledge of surface topography, baseline erro r, stratified at mosphere effects, temporal at mospheric heterogeneities, and thermal expansion, thereby isolating phase differences due to line -of-sight ground motion.All processing was conducted at MacDonald, Dettwiler and Associates (MDA) in Rich mond, Brit ish Colu mb ia, using programs developed by MDA on top of GAMMA RS modules and scripts for data ingestion and basic interferometric processing.

Pre-processing
After resamp ling all single -look co mp lex images to align with a master image acquired on 12 Ju ly 2010, we produced multi-looked images (2x3 looks), which p rovide the highest possible resolution for roughly square grounddimension pixels.We further produced scaled, lower resolution (8x12 looks) versions of these images suitab le for atmospheric modeling.
We stitched and cropped three-arc-second Shuttle RADA R Topographic M ission V3.0 (SRTM3; Farr, 2007), reference digital terrain model (DTM) t iles to fit the scene footprint.The SRTM 3 subset for the La Paz area is free of holes and other obvious abnormalities.We projected the DTM to SA R coordinates, converted elevation fro m a geoid to ellipsoid vertical datum, and finally calculated topographic phase for single -look and both multi-look complex phase data.

Differential InSAR
We applied differential InSA R (D-InSAR) to determine t ime-series phase statistics and to identify phase differences on scales greater than ~10 km, wh ich are presumed to be at mospheric phase contributions.For each unique acquisition pair, we generated a two -dimensional (2D) interferogram over the entire scene area for both the intermediate (2x3) and low-resolution (8x12), mu lti-look dimensions.The numbers of unique interferograms were 496 and 66, respectively, for the 32-scene pre-failure stack and the 12-scene post-failure stack.Topographic phase calculated fro m the SRTM DTM and flat-earth phase were removed, yielding init ial 2D interferograms corrected for terrain and Doppler effects.We also calculated spectral coherence for each unique image pair.
We then applied a coherence-maximizat ion algorith m to each interferogram in wrapped-phase space to identify and correct errors in each of the four co mponents of the interfero metric baseline (parallel baseline, perpendicular baseline, and the velocities of each), as well as parameters of a static at mosphere model.The algorith m imp roves initial estimates of orbital parameters for each of the generated 2D interferograms.
We applied a network inversion solution and smoothing function to the coarsest multi -look interferograms to characterize large spatial-scale phase differences, presumed to result fro m at mospheric heterogeneities.We masked out areas lacking coherence (e.g.water bodies) or with suspected ground motion (e.g.known landslides or areas of small spatial-scale phase fringes remaining after baseline optimization) to remove them in modeling the atmospheric phase.To reduce the likelihood of including unknown (a -priori to InSA R processing) ground deformation wh ile also allo wing adequate variation within the scene, we chose the scale of the smoothing function to be between the scale of deformation and half the scale of the scene (i.e., between ~3 km and 12.5 km).We removed the generated phase screens from baseline-optimized interferograms to create atmosphere-corrected interferograms.
Three main effects influence the remaining coherent phase: high -frequency topographic variability not adequately represented in the moderate-resolution DTM, cyclic movement related to thermal expansion, and permanent ground movement.Decorrelation also remains as an incoherence phase because the atmospheric screens do not contain, and thus do not remove, high-frequency noise.

Homogeneous Distributed Scatterer InSAR
We performed HDS-InSA R processing on a subset area (~ 3.4 x 4.1 km) near the center of the scene covering Pampahasi Plateau and the adjacent part of the Río Irpavi valley.As for D-InSAR processing, all steps were completed independently for the pre-failure and the post-failure stacks.Considering a spatial subset of single-look images helped facilitate t imely processing throughout the remaining processing chain and limited the co mputational resources required.
'Multi-looking'a form of spatial filtering within a weighted kernelimp roves estimates of interfero metric phase and image correlation (Parizzi and Brcic, 2011).HSD -InSA R emp loys a continuously weighted, spatially adaptive filter to group and average contiguous pixels that show similar time -series amplitude variability.The filter uses a similarity measure derived fro m a non-parametric statistical test (Anderson Darling) for thick scene stacks, such as the pre-failu re stack, or a one-parametric statistical test (Generalized Likelihood Ratio Test assuming Rayleigh distribution) for smaller stacks (fewer than ~15 scenes; Rabus et al., 2012), such as the post -failure stack.The resulting neighbourhoods display spatially ho mogeneous spectral behaviour over time.Due to their point nature, persistent scatterers yield single-pixel neighbourhoods.Coherent distributed scatters are characterized as spatially continuous neighbourhoods, limited in extent by the dimensions of the filtering kernel applied and represented by a point at the kernel's centre.We applied a 13 x 21 p ixel kernel to the stacks, which restricts spatial averaging to a roughly square ground area measuring ~105 m on a side.The adaptive filtering improves speckle and phase while preserving sharp boundaries and point targets.
We generated 2D interferograms for the sub-set area using refined baselines and the atmospheric screen developed during full-scene multi-look processing.These interferograms included each single-look network interferogram and each multi-looked (i.e., filtered using the previously defined weighted neighbourhoods) interferogram.The later are HDS differential interferograms that provide phase change for each HDS candidate.
For every interfero metric scene pair, we co mpared the differential p hase of each pixel to a reference phase of nearby pixels just outside its neighborhood candidate window (i.e., just beyond 13 x 21 pixels for these stack).We then calculated the overall differential coherence of the pixel throughout the scene stack, which provides a measure of the quality of coherence of each HDS neighbourhood by characterizing how much a pixel varies over time fro m adjacent pixels outside its neighbourhood.
Solving and correcting for residual height for each HDS led to significant improve ment of coherence prior to phase unwrapping and network inversion processing.As HDS -InSA R relies on prior unwrapping of the network, it contains a subsequent correction step that utilizes network redundancy to eliminate phase unwrapping errors.During this step, an init ial network inversion of the unwrapped phase identifies as outliers all p ixels that were incorrect ly unwrapped in some of the layers of the network.Iterative removal of the outliers by the inversion process then allo ws calculation of their integer 2 residuals caused by the phase errors in the corresponding network layers.
Phase unwrapping correction is completed by correcting the residuals for the affected pixels.
Search functions estimate the height error and linear deformation rate for each p ixel of each interferogram.For each factor considered, we calculated the peak-to-peak d ifferencethe difference between the highest and second highest maxima found during the searchfor each pixel.We then co mbined coherence maps of the differential phase components into an overall coherence map that shows the differential phase components with the lowest peak-topeak difference for each pixel location.Th is map shows the quality of the residual differential phase component at each pixel location.
We applied the aforementioned search function to at mosphere-corrected HDS interferograms, but using narrow search limits.This modeling estimates residual height corrections resulting fro m phase correlated with perpendicular baseline, and linear deformation rates resulting from temporally consistent phase.The former co mponent enables correction of the reference DTM to represent true terrain.In addit ion to height corrections and linear deformation rates, the modeled phase includes residual phase time series.
We separated spatially coherent neighbourhoods from incoherent ones using the temporal coherence of the differential phase as a measure of quality.In separate trials, we applied mu ltip le thresholds to the linear deformation component of the phase model; a higher threshold decreases the number of HDS candidates remaining.We determined the optimal coherence level by visually co mparing incremental threshold trails and adopted the coherence value that provided the best balance between eliminating spat ially variab le deformat ion trends among HDS and preserving high-quality-phase HDS.We ext racted the phase model and network interferograms for each HDS candidate above this coherence threshold.
We further flattened interfero metric phase through phase demod ulation.We co mputed the appro ximate deformation signal fro m the 2D InSAR solution of a minimu m network of h igh -quality interferograms.The remaining linear deformation rate was calcu lated after removing this deformation trend fro m the overall phase.Phase demodulation has two important benefits.First, it mitigates temporal wrapping due to large temporal baselines in stacks with many acquisitions over a long period.Second, it improves preservation of non -linear deformation for which the onedimensional linear deformation fit during phase modeling performs poorly.
We further flattened the demodulated phase using the phase model for residual height, linear displacement, and thermal dilation.The resulting corrections of each HDS for topographic height error relat ive to the reference DTM, and temperature dilation if selected, improve coherence.We then unwrapped the flattened phase with a two -step iterative variant of the minimu m-cost-function (Costantini, 1998) unwrapper, after which we added back the modeled phase components onto the HDS.We added back the deformation trend removed prior to spatial phase unwrapping onto the unwrapped HDS of each network interferogram to restore the complete deformation signal.
Using linear least squares regression, we imp roved the phase model by estimating residual height corrections and linear deformat ion coefficients remaining in the unwrapped data.Residual height errors determined by the refined phase model were removed from the unwrapped phase for each interferogram.
We inverted the network interferograms using singular value decomposition (SVD) to generate an optimu m time series of interferograms referenced to a common master date.Instead of a single spatial reference point or region, phase changes in HDS-InSAR are determined relative to a cu mulat ive reference area co mprising all areas lacking movement on scales <10 km, with broader scale motion having been removed by atmospheric filtering .We smoothed the SVD-inverted time series by temporal filtering and then calculated line-of-sight displacement.The temporal filter is a least-square, low-pass filter and is optimal for irregularly sampled data, such as result fro m temporal gaps due to missing scenes.The filter produces weighted -average smoothing over a five-scene period spanning two scenes on both sides of each time -series position.The s moothing period is thus a minimu m o f 120 days for RADARSAT-2 due to its 24-day return interval.

Data output generation
We georeferenced the temporally smoothed HDS-InSAR time series for each point location using the residual height errors estimated fro m the phase model inversion to improve the input reference DTM.The high density of data points (e.g.Fig. 7B, C) allowed us to interpolate a linear deformation map (e.g.Fig. 7E, F) fro m HDS to provide a base image to efficiently represent long-term, average line-of-sight deformat ion over the entire area.Small HDS clusters (fewer than five points) receive no weight during interpolation, resulting in suppression of much of the localised (smaller than several metres by several metres in the case of the RADA R resolution we processed) variations that might be noise resulting fro m phase unwrapping errors .Finally, we exported HDS point data and the linear deformation map to a GoogleEarth-compatible plug-in (Fig. 7).

Alignment of pre-failure and post-failure displacement trends
Due to the division of scenes into two temporally discrete stacks, the magnitude of motion between the last prefailure scene and first post-failure scene is unknown.For locations that did not experience appreciable movement during the 2011 failure event (i.e.those beyond the limits of the failure) , we assume that motion between the two stacks was limited and thus used results of full-stack processing (using all 44 scenes) to align pre-failure and postfailure d isplacement data (Fig. 9A; 'iv' and 'v ', Fig. 9B).To reduce the effect of temporal filtering, wh ich can negatively affect data at the ends of stacks due to a lack of temporal support, we matched displacement trends fro m the middle portion of the partial stacks to the same trends in the full stack.For areas within the failure area, ground displacements during the failure event were on the order of metres to h undreds of metres (Fig. 2B) and thus well beyond the measurement threshold of spaceborne InSAR.For such locations the precise amount of displacement between stacks is unknown (represented by line break in Fig. 9A for 'i', 'ii' and 'iii').

Estimating true displacement from line-of-sight displacement
One-dimensional displacements measured in the RADAR line of sight record only one component of the true ground displacement vectors.Consequently, the line-of-sight displacements underrepresent true displacement magnitudes unless these two directions are parallel, in wh ich case the true displacement equals the line -of-sight displacement.We estimated true displacement magnitudes by calculating the angle between the RADAR line of sight and the expected ground-movement direct ion to determine what proportion of the overall three -dimensional ground displacement was recorded by RADA R phase shifts.For the slopes descending east from Pampahasi Plateau, we assumed that the true deformat ion vector is parallel to the down-slope (fall-line) direction.The average downslope direction of these slopes was then used to calculate the approximate true displacement rates in Figure 8A.
The displacement direction of the plateau surface behind the headscarp is uncertain.The expe cted declination (horizontal co mponent) of the dilation-induced motion of the plateau surface is likely to be perpendicular to the local headscarp.However, its inclination is neither horizontal nor vert ical.True d isplacement rates cannot, therefore, be estimated and are not suggested in Figure 8B.

Comparison with InSAR characterization of the 2012 Via Piave landslide, southern Italy
To our knowledge, the only other orbital InSAR analysis that characterizes landslide act ivity both prior to and following a slope failure is for the 2012 Via Piave landslide in southern Italy (Confuorto et al., 2017).Co mparing the results of that study with those reported here highlights the unity of the HDS -InSA R technique as well as the suitability of environ mental conditions at La Paz for InSAR analysis.Furthermo re, the occurrence of the Pampahasi landslide between successive RADA R acquisitions allowed us to co mpare g round motion soon after the failure with that shortly before it.In contrast, substantial gaps (exceeding one year) in RADA R acquisit ions before and after the Via Piave landslide preclude the opportunity to reliab ly quantify slope-activity changes immed iately before and after the failure.

Figure S1 :
Figure S1: Locations of historic landslides exceeding 1 Mm 3 in volume in the La Paz area.Landslide numbers correspond to those in TableS1.See Figure1Bfor context within the La Paz area.

Figure S 2 :
Figure S 2: The precipitation record for the Laykacota meteorological station (see Fig. S 1 for station location) during the full period of InS AR monitoring.
the look d irection of the satellite (ground -range resolution) d iffers across the study area due to slope angle variations in the satellite look d irection and the systematic increase in the incidence angle from near range (south-southwest) to far range (north-northeast).Average ground-range resolutions for Pampahasi Plateau (horizontal in the satellite look direct ion) and the Pampahasi slope descending east from the plateau (11° average slope facing away fro m the satellite look d irection) are 8.8 m and 7.1 m, respectively.Azimuth resolution across the study area is 7.7 m (Table Confuorto et al. (2017) applied PS-InSAR and SBAS-InSAR to TerraSAR-X StripMap images covering the village of Papanice in the Crotone municipality.The scenes span two time periods separated by a three -year gap during which two landslides occurred.The pre-failure stack comprises 66 scenes acquired between April 2008 and June 2010, and was previously analysed byConfuorto et al. (2015).The two landslides happened on 23 February 2012, following two days of heavy rain.The post-failure stack co mprises 34 scenes acquired between October 2013 and October 2014.Although the RADA R scenes used byConfuorto et al. (2017) are of h igher no minal resolution (3 x 3 m) than those used in our Pampahasi study (5.2 x 7.7 m, TableS6), spatial characterizat ion of slope motion in the Italian case study was hindered by a low density of deformat ion records.That analysis focused on the smaller of the two 2012 Papanice landslides (Via Piave landslide; ~2 ha), seeming ly in part because of a lack of usable permanent scatterers on the larger (~8 ha) landslide.The s mall nu mber of reflectors identified o n the Via Piave landslide in the pre-failure (~15 points) and post failure (~35 points) stacks, and their concentration in the head region limit insight into the spatial extent of activity and hinder characterization of the failure mode.An apparent, although very slight, increase in displacement rates in the head region of the Via Piave landslide after 2012 may indicate a post-failure decrease in slope stability similar to that reported here for the 2011 Pampahasi landslide.Confuorto et al. (2017, p. 51) su mmarize their post-failure deformat ion rates as analogous to the prefailure d isplacement rates (30-40 mm/a LOS) reported byConfuorto et al. (2015).They report slight increases in RADA R LOS d isplacement rates in the northwest sector of the landslide (p art of the headscarp) fro m 36 mm/a well before the failure(Confuorto et al., 2017, p. 62)  to 40 mm/a well after it (p.64).The numbers of reflectors present on other parts of the landslide in both the pre-failure and post-failure InSA R scenes are small.Consequently, changes in activity in the central and lower parts of the landslide are unknown.
). g Sites visited in the field denoted by '*'.

Table S 3: Times and magnitudes of the wettest 15 days recorded at the Laykacota meteorological station between 1945 and 2016.
Note: T he location of the Laykacota meteorological station (16.5044° S, 68.1231° W) is shown in FigureS1.

Table S
4: Comparison of rainfall (mm) during each year covered by InS AR monitoring.Note: T he location of the Laykacota meteorological station (16.5044° S, 68.1231° W) is shown in FigureS1.