Geodetic Study of the 2006–2010 Ground Deformation in La Palma (Canary Islands): Observational Results

: La Palma is one of the youngest of the Canary Islands, and historically the most active. The recent activity and unrest in the archipelago, the moderate seismicity observed in 2017 and 2018 and the possibility of catastrophic landslides related to the Cumbre Vieja volcano have made it strongly advisable to ensure a realistic knowledge of the background surface deformation on the island. This will then allow any anomalous deformation related to potential volcanic unrest on the island to be detected by monitoring the surface deformation. We describe here the observation results obtained during the 2006–2010 period using geodetic techniques such as Global Navigation Satellite System (GNSS), Advanced Di ﬀ erential Synthetic Aperture Radar Interferometry (A-DInSAR) and microgravimetry. These results show that, although there are no signiﬁcant associated variations in gravity, there is a clear surface deformation that is spatially and temporally variable. Our results are discussed from the point of view of the unrest and its implications for the deﬁnition of an operational geodetic monitoring system for the island. J.E., J.F., J.F.P., A.G.C. and M.P.; writing—original draft preparation, J.E., J.F., J.F.P. and A.G.C.; writing—review and editing, J.E., J.F., J.F.P., A.G.C., M.P., G.R.-V. and E.A.; project administration, J.F.; and acquisition, J.F. and A.G.C.


Introduction
The detection of the accumulation and ascent of magma from depth can have a practical application providing advance warning of future eruptive activity [1]. A basic tool for this is to use geodetic deformation and gravity data [2,3]. Their use, together with appropriate inversion techniques, allow exploring the geometry, volume and location of magmatic reservoirs along the volcano plumbing systems, thus helping to understand the mechanisms and characteristics of unrest and eruption [4]. Recent decades have seen an explosion in the quality and quantity of volcano geodetic data [3], in particular using satellite radar interferometry (InSAR). This technique gives the option to have high precision and high spatial resolution, spanning decades, measurements of deformation for a large option to have high precision and high spatial resolution, spanning decades, measurements of deformation for a large number of active volcanic areas and allows carrying out a time study of magma movements along their plumbing systems. Many volcanic active regions are characterized by complex patterns of ground deformation resulting from multiple natural (e.g., inflation, deflation, dike intrusion, active faulting, flank instability and landslides) and anthropogenic sources [5][6][7][8][9]. Therefore, geodetic data are usually applied to study the main natural and anthropogenic sources of deformation as well as their associated hazards on active volcanic areas. Considering this, we applied geodetic observation techniques to study La Palma, Canary Islands, during 2006-2010.
The Canary Islands are located in the northwestern part of the Nubian (African) tectonic plate, relatively far from major plate boundaries and close to the thickened Western African craton continental lithosphere. They comprise a group of seven major islands forming a roughly east-west trending archipelago, and dated rocks show an age progression from Fuerteventura and Lanzarote in the east to La Palma and El Hierro Islands in the west (Figure 1). In general, older islands offer clearer evidence of longer erosional periods.  landslide and to study the most appropriate geodetic monitoring approach for the island. This paper therefore continues the work started by Fernández et al. [11] on La Palma. We update the works (Table 1) prior to 2010 and use an A-DInSAR processing technique including an error estimation for the line of sight (LOS) mean velocity and deformation time series to determine the background deformation status. More specifically, we use the satellite C-band SAR data from the ENVISAT image archive of the European Space Agency (ESA) acquired in ascending and descending orbits, supplemented with GNSS and gravity measurements collected during the same time interval. This study is the first part of a more ambitious work that also includes the interpretation by state-of-the-art inversion techniques of the deformation measured with A-DInSAR techniques integrated with the 3D crustal density structure of the island.

Reference Time Period Studied Technique Used Main Results
Moss et al. [5] 1994-1997 Observation of a geodetic network on the western flank of Cumbre Vieja using EDM and rapid-static GPS.
A coherent pattern of displacement vectors but with no statically significant values [4,22] Massonet and Sigmundsson [32] 1992-1995 Classical InSAR without topographic correction using only 1 pair of ERS radar images.
No deformation of the order of tens of cm. Smaller deformations could be masked by topographic fringes.
No deformation at the level of tens of cm in LOS. Smaller deformations could be masked by topographic fringes.
Subsidence on the Teneguía volcano, where the last eruption (1971) took place. No deformation on the northern part of the island.
Prieto et al. [28] 2006-2008 GNSS surveys with a geodetic network comprising 27 stations Subsidence and horizontal displacement detected in the area of the Teneguía volcano.
Arjona et al. [34] 1992-2008 A-DInSAR using CPT and ERS and ENVISAT radar images Negative LOS displacement in the Teneguía volcano area and on the western flank of the Cumbre Vieja volcano.
González et al. [23] 1992-2000; 2003-2008 DInSAR using stacking with ERS and ENVISAT radar images Two clear negative LOS displacements in the Teneguía volcano area and on the western slopes of the Cumbre Vieja volcano.

GNSS
To complement the previous geodetic GNSS works carried out on La Palma Island prior to the 2000s [2], a dense spatial sample of the island's ground deformation field was obtained from the La Palma episodic GNSS Monitoring Network established in 2006 [28,35]. It includes the national high-precision REGENTE geodetic network of seven sites, originally designed as a geodetic reference frame [11,36], and a set of new geodetic markers for monitoring volcano and geodynamic activity associated with the geology and structure of the island (Figure 1).
The geodetic network was designed to cover the areas of greatest susceptibility to potential significant surface deformation [37]. These are the rift-type structures and polygenic volcanic edifices such as the north-south rift of Cumbre Vieja; the areas of tensional stress caused by the continued injection of dikes on both sides of these rift structures; and calderas and valleys caused by debris avalanches such as Taburiente and Cumbre Nueva, as well as possible emerging calderas such as the hypothetical caldera resulting from the failure of the eruption of the San Juan volcano in the Llano del Banco area (Figure 1), after the 1949 eruption [31].
In addition to all these areas, the GNSS monitoring network covers the entire island with 26 GPS stations [28], as shown in Figure 1. There are two different types of geodetic network markers. The first Remote Sens. 2020, 12, 2566 5 of 23 type corresponds to the stations belonging to the former IGN network which consist of 1.20-m concrete pillars on top of a 1.00 m square concrete foundation anchored to the rock (Figure 2a). The pillars are equipped with a centering device that allows us to set the GNSS antenna with repeatability of over 1 mm. The second type of monumentation of the new monitoring stations consists of a benchmark steel bolt anchored to the rock (Figure 2b), as used routinely in other geodetic studies [38].
Remote Sens. 2020, 11, x FOR PEER REVIEW 6 of 24 introduced in the processing. Station coordinates and Earth orientation parameters were estimated with daily double differenced GNSS phase observations. The observations were weighted according to the elevation angle, with a cut-off angle of 10°. Absolute IGS antenna phase center and atmospheric zenith delay models [44] were used, and Global Mapping Functions [45] were adopted for the neutral atmosphere. The results of this processing step are daily estimates of loosely constrained station coordinates and other parameters, along with the associated variance-covariance matrices. In a subsequent step, the loosely constrained daily solutions were used as quasi observations in a Kalman filter (GLOBK) to estimate a consistent set of daily coordinates (i.e., time series) for all sites. All the time series were inspected to eliminate anomalous solutions (based on the maximum uncertainty) and to correct possible jumps related to antenna changes. The time series for each station are shown in Figures S1-S26. As a final step, the loosely constrained solutions and their full covariance matrices were combined using the GLORG module of GLOBK to estimate a consistent set of positions and velocities in the ITRF2008 reference frame [46], by minimizing the horizontal velocity of the 15 continuously operating global tracking stations. To adequately show the crustal deformation pattern over the study area, we aligned our estimated GNSS velocities to a fixed Nubian reference frame (see [47] for details). Annual mean velocities and their associated uncertainties are shown in Figure 3 and summarized in Table S1.
It should be noted that the deformation detected at the GNSS stations is, as a general rule, nonlinear. However, when the GNSS time series of these deformations are studied and compared ( Figure 4 and Figures S1-S26), it can be seen that they depart very little from the linear trend. This is because the deformations are not large, and the interval studied is limited to five years. For this reason, the annual mean velocities estimated from the time series represent well the general movements that have taken place during this period.

A-DInSAR
Advanced Differential Synthetic Aperture Radar Interferometry (A-DInSAR) is a microwave remote sensing technique that allows us to investigate surface deformation phenomena with a centimeter to millimeter precision and with a large spatial coverage capacity [48,49]. In particular, the A-DInSAR technique uses the phase difference, often referred to as an interferogram, between two suitable SAR images to temporally separate observations in a study area, and it provides a measurement of the ground deformation projection along the radar line of sight (LOS) [50]. For our A-DInSAR analysis, we used Synthetic Aperture Radar (SAR) data from the ESA's ENVISAT satellite to retrieve ground deformation maps in the time period covered by the four GNSS campaigns. In all four campaigns, at least two different sessions were carried out per station to minimize systematic local and/or user errors. For the bolt anchored stations, leveled fixed-height poles were used as antenna monument [38] to ensure a height repeatability of 1 mm with a random horizontal error of less than 2 mm [28]. The observations covered periods of 5-8 h, depending on the baseline lengths involved in each session.
The first GNSS data results were obtained by Prieto et al. [28]. Here, we report a new and comprehensive processing of all the acquired GNSS data by using the GAMIT/GLOBK 10.4 software [39,40] and by taking into account precise ephemerides from the IGS [41] and Earth orientation parameters from the International Earth Rotation Service [42]. To improve the overall configuration of the network and link the regional measurements to an external global reference frame, data from over 15 continuously operating global tracking stations, mainly permanent IGS networks [43], were introduced in the processing. Station coordinates and Earth orientation parameters were estimated with daily double differenced GNSS phase observations. The observations were weighted according to the elevation angle, with a cut-off angle of 10 • . Absolute IGS antenna phase center and atmospheric zenith delay models [44] were used, and Global Mapping Functions [45] were adopted for the neutral atmosphere. The results of this processing step are daily estimates of loosely constrained station coordinates and other parameters, along with the associated variance-covariance matrices.
In a subsequent step, the loosely constrained daily solutions were used as quasi observations in a Kalman filter (GLOBK) to estimate a consistent set of daily coordinates (i.e., time series) for all sites. All the time series were inspected to eliminate anomalous solutions (based on the maximum uncertainty) and to correct possible jumps related to antenna changes. The time series for each station are shown in Figures S1-S26. As a final step, the loosely constrained solutions and their full covariance matrices were combined using the GLORG module of GLOBK to estimate a consistent set of positions and velocities in the ITRF2008 reference frame [46], by minimizing the horizontal velocity of the 15 continuously operating global tracking stations. To adequately show the crustal deformation pattern over the study area, we aligned our estimated GNSS velocities to a fixed Nubian reference frame (see [47] for details). Annual mean velocities and their associated uncertainties are shown in Figure 3 and summarized in Table S1.  Envisat satellite described a Sun-synchronous low-Earth orbit. This type of orbiting satellites has an ascending pass that occurs when the satellite moves northward toward the equator and a descending pass, on the other side of the planet, where the satellite moves southward. Since A-DInSAR technique retrieves the displacement in one dimension (in the LOS direction), studying both geometries allow to retrieve more information on the deformation.
The descending orbit of the same SAR dataset was previously used in other studies (see Table  1). González  In this study, the A-DInSAR analysis was performed using the Subsidence software, based on the Coherent Pixels Technique (CPT) [51]. This software allows us to improve on previous results in several ways: by studying two geometries (ascending and descending) and including a deformation time-series for each coherent pixel within our study area and an error estimation of the velocity results [52]. The CPT algorithm and Subsidence software have been successfully used to perform multiple studies of different scenarios [8,[53][54][55]. We also increased the time period considered, from 2004 to 2010, and the number of radar images. The use of a more advanced technique and a larger number of images covering a longer time period, led to more precise results.

SAR Dataset
The dataset used in this study comprises 68 ENVISAT Single Look Complex (SLC) images: 42 obtained in descending orbit and 26 in ascending orbit. These images were acquired and provided by ESA and the dataset includes all the available images for this sensor, acquisition mode (IS2) and our region of interest. All the images are from tracks 359 and 431 for descending and ascending orbits, respectively. The list of all the SLC images used is given in Table S2. ERS-2, a C-Band satellite launched in 2002 and operational until 2011, products were also provided by ESA for the same tracks and frames. They were originally considered for inclusion in It should be noted that the deformation detected at the GNSS stations is, as a general rule, nonlinear. However, when the GNSS time series of these deformations are studied and compared ( Figure 4 and Figures S1-S26), it can be seen that they depart very little from the linear trend. This is because the deformations are not large, and the interval studied is limited to five years. For this reason, the annual mean velocities estimated from the time series represent well the general movements that have taken place during this period.
Remote Sens. 2020, 11, x FOR PEER REVIEW 9 of 24 due to their approximately 1-km correlation distance. However, the atmospheric contribution of a given pixel can be considered as a white process in time, as atmospheric conditions can be considered a random variable for each acquisition date. For instance, troposphere characteristics vary from one day to another. However, nonlinear movements have a narrower correlation window in space and a low-pass behavior in time. In view of all these considerations, atmospheric artifacts can be estimated with a filtering process in both the spatial and time domains [57,59,60]. A low-pass temporal filtering was applied to nonlinear estimation results to further reduce residual noise that was not removed in the second step of the processing (estimation of the nonlinear deformation component).

A-DInSAR
Advanced Differential Synthetic Aperture Radar Interferometry (A-DInSAR) is a microwave remote sensing technique that allows us to investigate surface deformation phenomena with a centimeter to millimeter precision and with a large spatial coverage capacity [48,49]. In particular, the A-DInSAR technique uses the phase difference, often referred to as an interferogram, between two suitable SAR images to temporally separate observations in a study area, and it provides a measurement of the ground deformation projection along the radar line of sight (LOS) [50]. For our A-DInSAR analysis, we used Synthetic Aperture Radar (SAR) data from the ESA's ENVISAT satellite to retrieve ground deformation maps in the time period covered by the four GNSS campaigns.
Envisat satellite described a Sun-synchronous low-Earth orbit. This type of orbiting satellites has an ascending pass that occurs when the satellite moves northward toward the equator and a descending pass, on the other side of the planet, where the satellite moves southward. Since A-DInSAR technique retrieves the displacement in one dimension (in the LOS direction), studying both geometries allow to retrieve more information on the deformation.
The descending orbit of the same SAR dataset was previously used in other studies (see Table 1). González et al. [23]  In this study, the A-DInSAR analysis was performed using the Subsidence software, based on the Coherent Pixels Technique (CPT) [51]. This software allows us to improve on previous results in several ways: by studying two geometries (ascending and descending) and including a deformation time-series for each coherent pixel within our study area and an error estimation of the velocity results [52]. The CPT algorithm and Subsidence software have been successfully used to perform multiple studies of different scenarios [8,[53][54][55]. We also increased the time period considered, from 2004 to 2010, and the number of radar images. The use of a more advanced technique and a larger number of images covering a longer time period, led to more precise results.

SAR Dataset
The dataset used in this study comprises 68 ENVISAT Single Look Complex (SLC) images: 42 obtained in descending orbit and 26 in ascending orbit. These images were acquired and provided by ESA and the dataset includes all the available images for this sensor, acquisition mode (IS2) and our region of interest. All the images are from tracks 359 and 431 for descending and ascending orbits, respectively. The list of all the SLC images used is given in Table S2. ERS-2, a C-Band satellite launched in 2002 and operational until 2011, products were also provided by ESA for the same tracks and frames. They were originally considered for inclusion in this study, but subsequently discarded due to high Doppler centroid frequency differences between the images for this period. This is due to the fact that a single gyroscope operation mode was used from February 2000 owing to the failure of five of the six onboard gyroscopes [56].
An area covering the entire island of about 33 km × 47 km was cropped from the original 100 km × 100 km SLC images. All SLC images were co-registered to make the dataset appropriate for subsequent A-DInSAR processing. To minimize perpendicular and temporal baselines, images 20070711 and 20080204, for ascending and descending orbits, respectively, were selected as reference images for co-registration [57].
Interferograms were generated using all possible pairs within 350 m perpendicular baseline and 400 days temporal baseline, discarding pairs with large perpendicular and temporal baselines (over 150 m and 150 days simultaneously). This criterion was used for both geometries. A Delaunay triangulation was also considered, although a larger set of interferograms was found to offer better results. This selection mode generated a total of 150 interferograms for descending orbit and 54 for ascending (see Tables S3 and S4 for the complete list). The number of descending interferograms was larger because of the higher number of SLC images, which reduces the temporal baseline between Remote Sens. 2020, 12, 2566 8 of 23 interferometric pairs. Since we were unable to use images acquired before 2006, our study needed to be reduced to cover the time span from March 2006 to April 2009 for ascending orbits and up to October 2010 for descending orbits. Figure S27 represents the spatial and temporal baselines obtained for the ascending and descending orbits.
An external high-resolution Digital Elevation Model (DEM) was used to remove the topographic phase from the interferograms. This DEM was generated with data from the Spanish public agency Instituto Geográfico Nacional (IGN), which provides these data for scientific and civilian use. The data consisted of orthometric heights for the entire La Palma Island and a geoid model (EGM08-REDNAP). A 25 m × 25 m ellipsoidal height DEM was generated with these data and used to remove the phase component due to the topography.

A-DInSAR Processing
We used the Coherent Pixels Technique algorithm (CPT) to obtain the surface displacements [51]. This is a coherence-based method that works with distributed scatterers at low-resolution over the multi-looked interferograms, similar to the widely-used Small Baselines Subset (SBAS) [50,58]. The geological characteristics of the surface and the available dataset make this kind of analysis more suitable than a full-resolution approach such as the Point Scatters (PS) method [59]. A mean coherence map was processed to establish a coherence-based pixel selection, using a multilook window of 15 lines in azimuth and three samples in range. This multilook results in low-resolution pixels obtained from the average of 45 pixels from the original interferogram, corresponding to 60 m × 60 m of ground spatial resolution.
A coherence criterion was chosen to select the pixels with sufficient phase quality to obtain the surface deformation. A medium coherence threshold of 0.35 was applied, corresponding to a phase standard deviation of 18 • [49]. This value provides good spatial coverage and sufficient phase quality to obtain a convergent solution. A Delaunay triangulation was used between pixels, and a limit of 800 m was established to reduce the atmospheric artifacts in the linear processing. To estimate linear velocity, CPT requires velocity and DEM error seeds, points with known velocity and altitude for the entire time period studied [51]. The GNSS stations JEDE and TIRI were used as velocity seeds, as their locations (on the eastern and western flanks of the Cumbre Vieja volcano) and the knowledge of their almost zero displacement rates from the GNSS measurements indicated that they would be good choices (see Figure 4 for the JEDE station and Figure S26 for the TIRI station). An additional seed was placed in Los Canarios village in the southern part of the island for descending processing, which was required to obtain LOS displacements in the Teneguía area due to a zone of very low coherence between this area and the TIRI seed. Points with low topography were selected for DEM error seeds (see Figure 5b).
The atmospheric contribution to the phase must be calculated to obtain the nonlinear velocity. Atmospheric perturbations are considered as a low spatial frequency signal in each interferogram due to their approximately 1-km correlation distance. However, the atmospheric contribution of a given pixel can be considered as a white process in time, as atmospheric conditions can be considered a random variable for each acquisition date. For instance, troposphere characteristics vary from one day to another. However, nonlinear movements have a narrower correlation window in space and a low-pass behavior in time. In view of all these considerations, atmospheric artifacts can be estimated with a filtering process in both the spatial and time domains [57,59,60].
A low-pass temporal filtering was applied to nonlinear estimation results to further reduce residual noise that was not removed in the second step of the processing (estimation of the nonlinear deformation component).  Velocity and DEM error seeds are marked as a rhombuses and triangles, respectively. Note that these images are in radar coordinates, and there is a mirroring effect along the horizontal axis due to the acquisition geometry in descending orbit.

A-DInSAR Results
La Palma Island is characterized by its high topography and areas of dense vegetation, which are the main reasons that A-DInSAR analysis is not optimal to study certain areas of the island. We processed the entire island in our study, but, due to the low number of pixels selected for the northern part of the island, we limited the results to the southern part, south of Caldera de Taburiente and the Cumbre Nueva volcano. Figure 5a shows the mean coherence of the multilook used (15 lines in azimuth and three samples in range direction); the northern part of the island is characterized by very low coherence values.
In our analysis, 19,082 and 10,533 pixels were selected in ascending and descending orbits, respectively. Some pixels shown in Figure 5b do not appear in Figure 6, as they were not assigned a linear velocity value. To obtain the linear velocity, each pixel must fulfil certain coherence criteria and be connected to a seed, and since no seeds are used in the northern part of the island the velocity of these points is not calculated.
The linear results for the region show a range of velocities between +1 and −2 cm/year in LOS direction. Positive values indicate pixels that move toward the satellite and negative values correspond to pixels that move away from the satellite. Note that some previous works discussed here [23] use the opposite sign convention.
The next step in our study consists of atmospheric filtering, calculating nonlinear displacements, and estimating statistical errors [52]. The final result is the time series of LOS displacements with an error estimation for each selected pixel in the region of interest. We selected three pixels for each geometry to show the results (see Figures 6 and 7). Pixels 1 and 4 present the maximum displacement rates on the flank of the Teneguía volcano. Pixels 2 and 5 are located in a stable zone during the period considered, on the western flank of the Cumbre Vieja volcano, with close to zero displacement. Pixels 3 and 6 are located in the northern part of the study area, in the southern part of the Cumbre Nueva volcano in the Aridane Valley, and also show significant LOS displacements. error estimation for each selected pixel in the region of interest. We selected three pixels for each geometry to show the results (see Figures 6 and 7). Pixels 1 and 4 present the maximum displacement rates on the flank of the Teneguía volcano. Pixels 2 and 5 are located in a stable zone during the period considered, on the western flank of the Cumbre Vieja volcano, with close to zero displacement. Pixels 3 and 6 are located in the northern part of the study area, in the southern part of the Cumbre Nueva volcano in the Aridane Valley, and also show significant LOS displacements.  indicate pixels that move toward the satellite and negative values correspond to pixels that move away from the satellite.
These results show that the uncertainties are very homogeneous over the study area, in the 0.25-0.35 cm interval range. This is a good result considering the number of images, the density of acquisitions for the time period and the perpendicular baselines between images.  Figure 6). Subfigures 2 and 5 show two stable points (pixels 2 and 5 on Figure 6). Subfigures 3 and 6 show the displacement of two pixels (3 and 6 in Figure 6) in the southern part of Cumbre Nueva with significant LOS displacement. Vertical axis corresponds with the displacement measured in cm. The left-hand column shows ascending orbit results and the righthand column descending orbit results. See Figure 6 for pixel locations.

Micro-Gravimetry Observations
To detect any possible signals related to additional masses or mass redistributions under the surface of La Palma, and taking into account the results reported in the literature [5,6,28,30,31], a micro-gravity network consisting of 12 stations was established in the SW part of the island and  Figure 6). Subfigures 2 and 5 show two stable points (pixels 2 and 5 on Figure 6). Subfigures 3 and 6 show the displacement of two pixels (3 and 6 in Figure 6) in the southern part of Cumbre Nueva with significant LOS displacement. Vertical axis corresponds with the displacement measured in cm. The left-hand column shows ascending orbit results and the right-hand column descending orbit results. See Figure 6 for pixel locations.
These results show that the uncertainties are very homogeneous over the study area, in the 0.25-0.35 cm interval range. This is a good result considering the number of images, the density of acquisitions for the time period and the perpendicular baselines between images.

Micro-Gravimetry Observations
To detect any possible signals related to additional masses or mass redistributions under the surface of La Palma, and taking into account the results reported in the literature [5,6,28,30,31], a micro-gravity network consisting of 12 stations was established in the SW part of the island and measured for the first time in June 2008 [35]. The geographic locations of the stations are shown in Figure 8, while their altitudes are represented in Figure 9. Coordinates are given in Table S5.  Four observation campaigns were carried out in June 2008, August 2009, July 2011 and June 2014. The network was measured ten times during each survey, over a five-day observation period. One full measurement was taken in the morning and another in the afternoon, and the total observation time was around 13 h. The stations were represented by paint marks on rocks or concrete (see Figure 10). The observation methodology was the usual one for microgravimetry surveys [12,61], and each observation was corrected considering calibration factors, tidal variations, drift and possible jumps in the recording. A precise GNSS survey was carried out in parallel with each gravimetric observation to obtain high-quality altitudes for the microgravimetry stations. Two GNSS survey methods were applied (e.g., rapid-static and real time kinematic), depending on the station conditions. Due to the different precisions of the two methods and the geoid model used (EGM08-REDNAP), a final adjustment was made with a different weighting strategy based on each variance factor, following the same procedure as [41]. This is a very controlled technique that has been proven in high-precision studies when calculating recorded observations and coordinates with different  Four observation campaigns were carried out in June 2008, August 2009, July 2011 and June 2014. The network was measured ten times during each survey, over a five-day observation period. One full measurement was taken in the morning and another in the afternoon, and the total observation time was around 13 h. The stations were represented by paint marks on rocks or concrete (see Figure 10). The observation methodology was the usual one for microgravimetry surveys [12,61], and each observation was corrected considering calibration factors, tidal variations, drift and possible jumps in the recording. A precise GNSS survey was carried out in parallel with each gravimetric observation to obtain high-quality altitudes for the microgravimetry stations. Two GNSS survey methods were applied (e.g., rapid-static and real time kinematic), depending on the station conditions. Due to the different precisions of the two methods and the geoid model used (EGM08-REDNAP), a final adjustment was made with a different weighting strategy based on each variance factor, following the same procedure as [41]. This is a very controlled technique that has been proven in high-precision studies when calculating recorded observations and coordinates with different Four observation campaigns were carried out in June 2008, August 2009, July 2011 and June 2014. The network was measured ten times during each survey, over a five-day observation period. One full measurement was taken in the morning and another in the afternoon, and the total observation time was around 13 h. The stations were represented by paint marks on rocks or concrete (see Figure 10). The observation methodology was the usual one for microgravimetry surveys [12,61], and each observation was corrected considering calibration factors, tidal variations, drift and possible jumps in the recording. A precise GNSS survey was carried out in parallel with each gravimetric observation to obtain high-quality altitudes for the microgravimetry stations. Two GNSS survey methods were applied (e.g., rapid-static and real time kinematic), depending on the station conditions. Due to the different precisions of the two methods and the geoid model used (EGM08-REDNAP), a final adjustment was made with a different weighting strategy based on each variance factor, following the same procedure as [41]. This is a very controlled technique that has been proven in high-precision studies when calculating recorded observations and coordinates with different levels of precision [62]. The adjusted relative gravity values for each survey are listed in Table S6. Based on the relative gravity values obtained, we opted to make a single global adjustment which allowed us to estimate simultaneously the corrections to the calibration factors for the different campaigns, the gravity variation rates (µGal/year) and the corrected gravity values at each station. All this was done without assuming any station to be stable, with a zero-variation rate (considering we have no available information that properly identifies a station as being stable throughout the study time period covered by the study) or any gravimeter as being perfectly calibrated.
Let gij be the relative gravity values obtained (after reducing the observation data) for each station i (i = 1, …, n, where n is the total number of stations, n = 12) and for each campaign j (j = 1, …, m, where m is the number of campaigns, m = 4). These values can be somewhat heterogeneous, and they are affected, among other things, by minor inaccuracies in the calibration factors used for the gravimeters. If we denote the homogenized gravity values as Gij, we can write: Let be the periods (in years) when the campaigns were carried out. Let be a mean reference date (for example 2010) and let be the homogenized gravity value corresponding to that central moment j = 0. In these circumstances, we can introduce the annual variation rate in each station as given by: The value is expressed as a function of the average value for the m campaigns, in the form: Combining Equations (1)-(3) we can write: where the unknowns are (i = 1,…, n), (i = 1,…, n) and (j = 1,…, m). The resulting system of equations is clearly ambiguous or range-deficient (an ill-posed problem) as no scale is fixed. We solve this problem using the Tikhonov damping factor technique, a method of regularization of ill-posed problems [63]. It is also known as ridge regression and is particularly useful to mitigate the problem of multicollinearity in linear regression, which commonly occurs in models with a large number of parameters. The method generally provides improved efficiency in parameter estimation problems in exchange for a tolerable amount of bias.
In essence, for the system of observation Equations (4), which we now write in vector form as − = , with being the design matrix, denoting the vector of unknowns, being the Based on the relative gravity values obtained, we opted to make a single global adjustment which allowed us to estimate simultaneously the corrections to the calibration factors for the different campaigns, the gravity variation rates (µGal/year) and the corrected gravity values at each station. All this was done without assuming any station to be stable, with a zero-variation rate (considering we have no available information that properly identifies a station as being stable throughout the study time period covered by the study) or any gravimeter as being perfectly calibrated.
Let g ij be the relative gravity values obtained (after reducing the observation data) for each station i (i = 1, . . . , n, where n is the total number of stations, n = 12) and for each campaign j (j = 1, . . . , m, where m is the number of campaigns, m = 4). These values can be somewhat heterogeneous, and they are affected, among other things, by minor inaccuracies e j in the calibration factors used for the gravimeters. If we denote the homogenized gravity values as G ij , we can write: Let t i be the periods (in years) when the campaigns were carried out. Let t o be a mean reference date (for example 2010) and let G i0 be the homogenized gravity value corresponding to that central moment j = 0. In these circumstances, we can introduce the annual variation rate v i . in each station as given by: The value G i0 . is expressed as a function of the average value for the m campaigns, in the form: Combining Equations (1)-(3) we can write: where the unknowns are v i (i = 1, . . . , n), δ i (i = 1, . . . , n) and e j (j = 1, . . . , m).
The resulting system of equations is clearly ambiguous or range-deficient (an ill-posed problem) as no scale is fixed. We solve this problem using the Tikhonov damping factor technique, a method of regularization of ill-posed problems [63]. It is also known as ridge regression and is particularly useful to mitigate the problem of multicollinearity in linear regression, which commonly occurs in models with a large number of parameters. The method generally provides improved efficiency in parameter estimation problems in exchange for a tolerable amount of bias.
In essence, for the system of observation Equations (4), which we now write in vector form as Ax − b = r, with A being the design matrix, x denoting the vector of unknowns, b being the vector of independent terms and r denoting the vector of residuals. Solution x is computed as: where I is the identity matrix (corresponding to the number of unknowns) and α > 0 is the Tikhonov factor. In this case, the usual minimization of Ax − b 2 is replaced by the minimization of Ax − b 2 − α x 2 , which, in addition to the size of the residuals, takes into account the size of the unknowns. The g ij data are the values contained in Table S6 for each of the n = 12 stations and each of the m = 4 dates. Applying the previously described adjustment process (and considering, after tests, an α = 1 value), we arrived at the values shown in Table S4 and Table 2. The root mean square of the adjustment is 5.9 µGal, which is quite acceptable given the physical conditions of the network (very rugged relief and connected by means of winding forest tracks), but perhaps insufficient for studying such subtle variations.  Figure 11 shows the estimated annual rate and standard error for each of the 12 gravimetric stations in the network. It can be seen that the adjusted annual rate values and the standard error values are generally small. Two stations have values exceeding the standard error: Pista 20 (6.5 ± 4.5 µGal/year) and Cancajos (4.6 ± 2.0 µGal/year), although neither clearly exceeds twice the value of the standard error. This points to two conclusions: (1) the quality of the gravimetric network is acceptable (estimated standard errors of the order of ± 5 µGal/year for the annual rate values and mean root square residual of 5.9 µGal), although perhaps insufficient due to the small magnitude obtained for the variations measured; and (2) none of the adjusted annual rate values is significant. No significant variations in gravity are therefore detected in the gravity network, at least at the level of precision obtained. standard error. This points to two conclusions: (1) the quality of the gravimetric network is acceptable (estimated standard errors of the order of ± 5 µGal/year for the annual rate values and mean root square residual of 5.9 µGal), although perhaps insufficient due to the small magnitude obtained for the variations measured; and (2) none of the adjusted annual rate values is significant. No significant variations in gravity are therefore detected in the gravity network, at least at the level of precision obtained.

GNSS
As is the case for many of the Canary Islands, it is not evident where future eruptions will take place in La Palma because their likelihood (with varying probability) is not limited to a specific area but affects the whole island [11,21,64,65]. Another phenomenon that must be monitored on the island is the possibility of displacements that could be associated with a landslide motion in the area of the Cumbre Vieja volcano [6,66,67]. It was therefore necessary to install a geodetic monitoring system that extensively covered the island [11,21,68] and was the reason for defining a GNSS network, as previously described.
Two important aspects must be considered when assessing the applicability of this network for monitoring displacements on the island. The first is related to the comparison of the observed/expected magnitude of the displacements (when these values are available) in terms of the precision of the observation technique (between 0.5 and 1 cm in the different displacement components), while a second aspect concerns the possibilities for interpreting these displacements.
Some conclusions can be drawn from a study of the LOS displacements detected using A-DInSAR, which are of the order of 0.5-1.5 cm/year in absolute values, and their locations with respect to the GNSS stations. The first is that we are already at the detection limit of the survey mode GNSS observation used in the network [28,37], and the second is that the separation between the GNSS stations, which varies from approximately 1 to 10 km, combined with the magnitude of the displacement itself, implies that the temporal GNSS observation of the network defined for the interpretation/inversion of the displacements is of very limited applicability during the initial stages of unrest. This technique would therefore be more suitable when the displacement is on the order of tens of centimeters and a longer decay wavelength.
These limitations, coupled with the cost of GNSS observations [69], pointed to A-DInSAR as the primary technique for monitoring deformation on the Island of La Palma in an operational way. The main applicability of the GNSS network at the present is therefore the validation and scaling of the displacements, comparing the GNSS and A-DInSAR results. It would be very useful if some of the network stations were with continuous recording to complement the campaign mode. Obviously, the station locations must be selected taking into account volcanological, structural, geological and geodynamic information as well as the results obtained from the radar observations. Unfortunately, this last aspect was not considered for the definition of the GNSS network in the study, as it was not available at the time of configuring the network.

GNSS-InSAR Comparison
It is interesting to note that, in the areas where the A-DInSAR results were obtained using ascending and descending radar images, the A-DInSAR results sometimes have a different sign for the LOS velocity in each geometric configuration. The signs of both LOS displacement velocities agree for the western flank of the Cumbre Vieja volcano, but they differ for the Teneguía volcano and Aridane Valley areas (see Figure 5). This is a very rough estimate because the LOS depends heavily on the look angle.
In a qualitative way, this clearly indicates that, in the first zone-the western flank of the Cumbre Vieja volcano-it is basically the vertical displacement component that contributes to the LOS and comprises mainly subsidence. This is confirmed by the GNSS results for the same area in Figure 3. This is not the case for the Teneguía area, where there is a contribution of both vertical and horizontal displacement components (E-W) in the LOS. This is again confirmed in Figure 3 for the GNSS stations located in this area. In the Aridane Valley area, there is only one station unable to represent the entire deforming area (see Figure 6).
To compare the GNSS and A-DInSAR information in a quantitative way, the three-dimensional variation of coordinates from the GNSS values for each station must be projected into LOS direction. We can do this using the equation [70]: where v x , v y , v z are the three components of the GNSS velocity, θ is the incidence angle and α is the azimuth angle of the satellite track. In our sign criterion, positive values of d LOS indicate that the pixel moves toward the satellite and negative values indicate that the pixel moves away from the satellite. Figure 12 shows both values for several GNSS stations from the network. It was basically observed from the error bars that the agreement is quite good, validating the conclusions described previously in the discussion of the GNSS results. Because SAR is sensitive in the perpendicular direction to its azimuth and the satellite moves along an almost polar orbit, the LOS displacement can be assumed to be produced by vertical and E-W motion, and the N-S motion is neglected. To obtain the vertical and E-W components of the displacement, the ascending and descending LOS motions must be combined in areas where both are present [8]. However, if we consider the magnitude of the displacements observed using A-DInSAR and GNSS (Figures 3 and 6), it is clear that the computation of this decomposition will not produce good (i.e., sufficiently precise) results due to methodological aspects, as well as because the magnitude of the E-W motion is at sub-centimeter level in many of the coherent pixels, or very close to zero, i.e., of the same order or even lower than the A-DInSAR uncertainty [8]. We therefore did not calculate it.
A comparison of Figures 3 and 6 shows that the GNSS and A-DInSAR observations are also complementary in terms of coverage. There are GNSS pixels in some areas where there are no coherent SAR pixels, although the number of stations is very limited. Considering the magnitude of the displacements observed, as previously mentioned, it would also be more convenient to use continuous GNSS observations in place of the survey observation mode in some very well selected areas to increase the precision.
In view of all the aspects discussed in this section, our results show that, as in other environments [8], the ad-hoc establishment of survey mode GNSS networks improves the spatiotemporal monitoring of the 3D displacement of the island subjected to volcano hazard and complements the A-DInSAR observations. As is well known [3], this is a valuable monitoring technique with a variety of applications at different stages of volcanic activity/unrest, ranging from complementing A-DInSAR (facilitating its validation and scaling) to data inversion, helping to obtain the characteristics of the deformation sources.
Because SAR is sensitive in the perpendicular direction to its azimuth and the satellite moves along an almost polar orbit, the LOS displacement can be assumed to be produced by vertical and E-W motion, and the N-S motion is neglected. To obtain the vertical and E-W components of the displacement, the ascending and descending LOS motions must be combined in areas where both are present [8]. However, if we consider the magnitude of the displacements observed using A-DInSAR and GNSS (Figures 3 and 6), it is clear that the computation of this decomposition will not produce good (i.e., sufficiently precise) results due to methodological aspects, as well as because the magnitude of the E-W motion is at sub-centimeter level in many of the coherent pixels, or very close to zero, i.e. of the same order or even lower than the A-DInSAR uncertainty [8]. We therefore did not calculate it.

Comparison of the DInSAR Results with Previous Ones
We now compare our A-DInSAR results with the previous results obtained from images from the ESA's ENVISAT satellite [23,34]. However, it should be noted that we use a different processing methodology, consider a different time period and process a different number of images. Each of these aspects can produce differences in the results.
González et al. [23] estimated only the mean velocity using stacking, without time series or error estimation and corrected the atmospheric contribution by applying elevation-phase dependence. Their study focuses only on descending geometry and a total of 19 ENVISAT SLC images acquired between 2004 and 2008, which they compared with the results of the first analysis of only three GNSS surveys carried out between 1994 and 2007. It should also be noted that, in the case of their estimates for the linear rate from the 2004-2008 period, the images were noisier, mainly due to the smaller dataset, thus they recommended they should be considered with caution.
Arjona et al. [34] used a previous version of the Subsidence software without error estimation, processing both ascending and descending geometries. They included 15 images acquired between 2004 and 2007 for the ascending orbit and 18 images in the period 2006-2008 for the descending one.
The present work used the most recent version of the Subsidence software, which offers the latest developments in atmospheric filtering, error estimation, interferogram formation and velocity calculation. We used the complete available ENVISAT catalog for La Palma and processed a greater number of SLC images than in previous works, with 25 (Table S2). Updated GNSS results obtained from five surveys carried out between 2006 and 2011 were used for comparison, in parallel to the radar observation. Table 3 shows a comparison of the results obtained in these three works.  Table 3 shows that, with a few exceptions, our findings agree with those obtained in previous works. Even considering that the results presented by Arjona et al. [34] were obtained using a smaller number of images and an older version of the Subsidence software, they concur with our result, except for the Aridane Valley area, which they described as stable while our processing shows displacements in both geometries. This difference may be due to the location of their seed. This area could be a good choice for locating a velocity seed since it is mostly flat with high coherence (large urban areas) and with no deformation signal detected by GNSS observation throughout the study period. Since the location of the seed is not described in their work [28], this is the most likely cause of the discrepancy between the results, coupled with the aforementioned differences between the two processing methods.
The discrepancies in the rest of the cases appear to be more closely related to the processing methodology or software version, the number of images and the time period. In this last aspect, if (in addition to Table 3) we consider the results of the ERS-1/2 images [23,30], it is clear that between 1992 and 2011 the magnitude and location of the LOS displacement rates (mean linear velocities) have changed across the study area, and they may be associated to some unrest in the island, as proposed by Torres-González et al. [20].

Microgravimetry Results
If there has been some unrest on the island, there may be new masses (magma) at shallow depths below the surface. This can be checked by using the gravity variation values obtained from the micro-gravity network observations. The possible variations in gravity could derive from: (i) free-air effects derived from vertical deformation produced by volcanic unrest, among other possible sources; or (ii) the contribution of additional masses, magma intrusions and/or density redistribution below the surface produced by internal and surface deformation.
The displacements detected using A-DInSAR in the gravity observation area (see Figures 3,6 and 8) show that the deformation in the gravimetric stations is of the order of few mm/year, except perhaps in the Cancajos-Hotel, Gasolinera and Teneguia 10-20 stations, where LOS displacement values of the order of 0.8 cm/year may be reached. Assuming that these are all elevation rate values, the gravity variation corresponding to this maximum displacement would be lower than 3 µGal, which is clearly below the measurement precision.
Despite the difficulties posed by the observations with the La Palma microgravimetry network (given the rugged nature of the relief, and the distances and slopes to be traveled from point to point along forest dirt tracks), the adjusted quality parameters suggest that the observational work was satisfactory, with a standard deviation of around 5 µGal. At this level of precision, no significant geographically extensive variations in gravity were detected. It is conceivable that the variations in gravity associated with density redistribution must be very small (the deformations are minor) and the expected changes in mass, if any, are therefore not large enough to be detectable. As we did not measure any significant gravity changes, non-significant additional masses (magma) appear to be present at the shallow crustal level in the study area and time period.

Conclusions
A review and follow-up of previous InSAR studies is justified to determine the background surface displacement field for La Palma Island considering the existing hazards and risks [65,71]. We therefore studied this phenomenon for the period 2006-2010, using several last-generation geodetic observational techniques simultaneously (GNSS, A-DInSAR and microgravimetry).
The observational results show a broad deformation but with low magnitude and rates of 0.5-2 cm/year, without any significant gravity changes. This implies that there is no additional mass (magma) in the shallow crust below the Cumbre Vieja Volcano.
These displacements are within the detection limit of the survey mode of the GNSS method used in the network. The magnitude of the displacement and the separation between the GNSS stations (1-10 km approximately) indicates that GNSS has limited application in the interpretation/inversion of the displacements during quiescent periods and early stages of unrest, but is more useful at more evolved stages with shallow intrusions, when displacements of the order of tens of centimeters with a longer wavelength decay may appear.
In view of these limitations and the cost of GNSS observations, A-DInSAR is the primary method recommended for an operational volcanic geodetic monitoring system on La Palma Island, using advanced techniques as seen in this work. The main applicability of the GNSS network observations during quiescent periods or in the early stages of unrest is currently the validation and scaling of the A-DInSAR displacements. It would also be useful for some of the stations in the network to carry out continuous observations to complement the survey mode. The observational results described here suggest that the station locations must be selected based not only on volcanological, structural, geological and geodynamical information, but also considering the results of the radar observations. The GNSS stations should be installed in the most effective locations for validating and complementing the A-DInSAR results, potentially leading to significant reductions in the cost of geodetic monitoring. The GNSS network used in this study should therefore be redesigned in some areas of the island, and according to the current distribution of continuous GNSS stations.
We obtained new A-DInSAR results for La Palma, including time series of LOS displacements and error estimation, using images obtained in ascending and descending orbits. It should be noted that in areas where we have deformation results, the LOS obtained using ascending and descending radar images sometimes have different signs, as in the case of the Teneguía volcano and Aridane Valley areas. However, the sign of both LOS displacement velocities is the same for the eastern flank of the Cumbre Vieja volcano. This provides qualitative evidence that it is essentially the vertical displacement component that contributes to the LOS on the eastern flank of the Cumbre Vieja volcano, and it is mainly subsidence, as confirmed by the GNSS results for the same area. However, for the other two areas-the Teneguía volcano and Aridane Valley-both vertical and horizontal displacement components (E-W) contribute to the LOS. This is confirmed by the GNSS results for the Teneguía volcano area, although in the case of the Aridane Valley there is only one station, and the results cannot be checked against the GNSS results. In any case, there is good agreement between the GNSS and A-DInSAR results at a regional scale.
The GNSS and A-DInSAR observations can also be seen to be complementary in terms of coverage. There are GNSS stations in areas where there are no SAR pixels, although the number of stations is very limited, and it is worth studying whether the observation cost in terms of the results obtained would compensate for the increase in the number of GNSS stations. In summary, our results again show that the ad-hoc configuration of survey mode GNSS networks combined with satellite radar observation can improve the spatiotemporal monitoring of the 3D displacement field of areas subjected to volcano hazard. Both observations techniques have different and complementary applications throughout the various stages of volcanic quiescence and unrest.
These findings, in conjunction with the results of previous works (see Table 1), clearly show that LOS displacement rates (mean velocities) and displacement time series for coherent pixels have changed over time from 1992 to 2011, and spatially along the surface of the study area. This may be associated to some unrest on the island, as recently proposed by Torres-González et al. [20]; this same conclusion is suggested by considering our results together with the degassing results obtained in La Palma during a similar time period. Padrón et al. [72] reported the results of three soil helium surveys undertaken at the Cumbre Vieja volcano (2002, 2003 and 2004), which support the hypothesis of contemporary degassing at Cumbre Vieja from magmas stored at different depths under the island.
An increase in the solid 3 He/ 4 He ratio was also observed toward the southern part of the Cumbre Vieja volcano [72]. It was interpreted as the result of residual degassing of volatiles from magma bodies stored at lithospheric levels beneath the southern part of the volcano, which were responsible for the last volcanic eruption of Cumbre Vieja in Teneguía in 1971. This result agrees with those obtained for the Teneguía volcano area both here and in previous works [23,30].
As mentioned above, this is the first part of a more ambitious work intended to also include the interpretation of deformation by means of state-of-the-art inversion, integrated with the three-dimensional crustal density structure of the island. This additional study is the subject of a future manuscript.