Long-Term Subsidence in Lava Fields at Piton de la Fournaise Volcano Measured by InSAR: New Insights for Interpretation of the Eastern Flank Motion

Long-term deformation often occurs in lava fields at volcanoes after flow emplacements. The investigation and interpretation of deformation in lava fields is one of the key factors for the assessment of volcanic hazards. As a typical Hawaiian volcano, Piton de la Fournaise volcano’s (La Reunion Island, France) main eruptive production is lava. Characteristics of the lava flows at Piton de la Fournaise, including the geometric parameters, location, and elevation, have been investigated by previous studies. However, no analysis focusing on the long-term post-emplacement deformation in its lava fields at a large spatial extent has yet been performed. One of the previous studies revealed an instability of the Eastern Flank since the March–April 2007 eruption related to post-emplacement lava subsidence. However, it was only a preliminary investigation. In this paper, an InSAR time series analysis is performed to characterize the long-term deformation in lava fields emplaced between 1998 and 2007 at Piton de la Fournaise, and to conduct an in-depth investigation over the influence of post-emplacement lava subsidence processes on the instability of the Eastern Flank. Results reveal an important regional difference in the subsidence behavior between the lava fields inside and outside of the Eastern Flank Area (EFA), which confirms that, in addition to the post-lava emplacement processes, other processes must have played a role in the observed subsidence in the EFA. The contribution of other processes is estimated to be up to ~78%. The spatial variation of the observed displacement in the EFA suggests that a set of active structures (like normal faults) could control a slip along a pre-existing structural discontinuity beneath the volcano flank. This study provides essential insights for the interpretation of the Eastern Flank motion of Piton de la Fournaise.


Introduction
Lava flows emitted during volcanic events can continue to deform for years after their emplacements, generating measureable ground displacement on and adjacent to lava flow fields [1][2][3][4]. superimposed on a shaded relief map with the main areas mentioned in the paper, modified from Servadio [9]: the Enclos Fouqué caldera, the Dolomieu and Bory caters, the Plaine des Osmondes, the Grandes Pentes, and the Grand Brûlé areas. Location of flank GNSS stations, GPNG, GPSG, HDLG, and GBNG, installed on these lava fields is indicated by black triangles. The inset shows the location of La Réunion Island in the Indian Ocean and that of the Piton des Neiges and Piton de la Fournaise volcanoes in La Réunion Island, with the location of the REUN IGS station used in this study for subtraction of the plate motion from GNSS data. Coordinates in Geographic Lat/Lon WGS84.
Interferometric Synthetic Aperture Radar (InSAR), a rapidly developed technique, has been widely and successfully applied to monitor the surface displacements at volcanoes worldwide, such as the Etna, Kilauea, and Piton de la Fournaise (e.g., [10][11][12][13][14][15]). However, the number of InSAR studies focusing on the displacement due to the post-lava emplacement processes is very limited with respect to those focusing on the displacement associated with magmatic processes [16]. Among those few InSAR studies detecting lava displacement after emplacement (e.g., [1][2][3][16][17][18][19][20]), the studies carried out at Piton de la Fournaise are scant [20][21][22] and none of them provide a detailed analysis and comparison on the long-term post-emplacement displacement for multiple lava fields.
Furthermore, in a previous study of Chen et al. [23], a widespread active seaward motion was revealed affecting the Eastern Flank of Piton de la Fournaise during the 2009-2014 period using InSAR, indicating an instability of the volcano flank after the March-April 2007 distal eruption. Along with this major investigation, a preliminary analysis on displacement in lava fields emplaced between March 1998 and April 2007 was carried out. Results suggested that the post-emplacement lava subsidence might have contributed up to~66% of the flank downward motion, providing useful information for understanding the flank instability. However, this result was just derived from a preliminary investigation. The displacement patterns related to the lava fields were roughly selected by visual inspection and the resultant values were simply calculated by mean measurements.
Therefore, in this study, we perform a detailed InSAR time series analysis to characterize the long-term displacement associated with lava fields emplaced between March 1998 and April 2007 at Piton de la Fournaise, with the objective to investigate in detail the influence of post-emplacement lava subsidence processes on the instability of the Eastern Flank. The spatial pattern and temporal behavior of the displacement of each lava flow are presented, and the relationship between displacement, lava thickness, and age are statistically analyzed. This study makes it possible to infer the location of two faults that may have played a role in the seaward flank motion, which provides a new perspective for the understanding and assessment of potential hazards associated with flank destabilization.

Deriving InSAR Long-Term Ground Displacement
To derive the long-term ground displacement in lava fields, we perform an InSAR time series analysis on a large amount of X band SAR dataset acquired by the German Space Agency TerraSAR-X/TanDEM-X (hereinafter named TSX/TDX) satellites constellation. The dataset includes 56 acquisitions from an ascending orbit, spanning from 7 March 2009 to 15 November 2014, and 34 acquisitions from a descending orbit spanning from 13 December 2008 to 6 October 2014 (Table 1). A small baseline (SB) method based on the Stanford Method for Persistent Scatterers (StaMPS; [24]) was adopted to implement the InSAR time series analysis. Interferograms with perpendicular baselines less than 500 m and temporal baselines less than 365 days were firstly processed using DORIS software. In order to ensure no isolated clusters remained in the networks, we then added interferograms with larger baselines but good coherence to obtain the final SB subsets for both ascending (containing 325 SB interferograms) and descending (containing 196 SB interferograms) orbits. An IGN (French National Geographic Institute) Lidar Digital Elevation Model (DEM) with 5 m resolution was used for coregistration and for simulation and subtraction of topographic contribution. The StaMPS selected the coherent pixels referred to as slowly-decorrelating filtered phase (SDFP) pixels based on a statistical analysis of amplitude difference dispersion and the phase stability [12,25], and adopted the 3D unwrapping algorithm [26,27] to unwrap the phases. In order to invert the InSAR time series corresponding to each SAR acquisition, we performed a pixel by pixel least squares adjustment independently on each SB subset, based on the following expression: where φ sb are the interferometric phases from SB subset, φ ts the time series to be inverted and G the design matrix with zeros, plus and minus ones [28]. The time series were temporally referenced to the first acquisition date and spatially referenced to non-displacement areas about 10 km away from volcanic craters. The non-negligible long-wavelength artifacts remained in the time series maps were properly modeled and subtracted using a correction method based on the principal component decomposition implemented in the Principal Component Analysis-based Inversion Method package (PCAIM; [29]). The correction method is demonstrated in detail by [23], and will not be repeated here in this paper. The InSAR ascending and descending measurements provide displacement along two different lines of sight. They can be combined to resolve a 2-dimensional (2D) displacement problem under an assumption that the third dimensional component is negligible. Here we generated the vertical (U-D) and horizontal (E-W) components of ground displacement using the formulation of Wright et al. [30]. This is reasonable because the real north-south (N-S) component of displacement is generally lower or not significantly larger than the U-D and/or E-W components over the study area, according to the daily measurements of 4 GNSS stations located on the Eastern Flank of Piton de la Fournaise ( Figure 2). Before inverting the 2D components of the ground displacement time series, we interpolated the ascending and descending data at missing epochs over their overlapped period. Figure 2 shows the U-D and E-W components of the cumulative displacement maps and the time series observed from InSAR data superimposed on the measurements recorded by four continuous GNSS stations. It is worth noting that the summit area is affected by a large transient displacement associated with eruptive activities (November 2009, December 2009, January 2010, October 2010 and December 2010, [23]) during the study period. The cumulative displacement in the summit area is evident in Figure 2b. Since this study focuses on mainly long-term displacement in lava fields on volcano flanks, the transient displacement in the summit area will not be discussed afterwards and the color scales were adjusted to the visualization of long-term displacement. The uncertainties of the displacements estimated by calculating the average standard deviations of the displacements in non-displacement areas are 3.9 cm and 3.1 cm for the E-W and U-D components respectively. A good agreement between the resultant InSAR time series and GNSS measurements has been reached. The RMS of the differences between both datasets are estimated to be~1.0 cm and~1.6 cm for the E-W and U-D components respectively, and the correlation coefficient is 0.96, revealing the reliability of InSAR processing. The map is spatially referenced to a non-displacement zone that is located~10 km away from volcanic craters. The dashed brown ellipse indicates the Eastern Flank Area (EFA) that is affected by both downward and eastward displacement during the study period. The red rectangle denotes southern part of the EFA where no lava fields emplaced between 1998 and 2007 but significant displacement occurred. The dashed blue curve in (b) indicates the eastern limit of the area affected by the October 2010 eruption, which is inferred according to the observed eastward displacement and the observation provided by Bato et al. [20]. F1 and F2 (black lines) represent the approximate location of two structures discussed in this study. Location of GNSS stations, GPNG, GPSG, HDLG, and GBNG, is indicated by black triangles. The cumulative displacements (U-D, E-W and N-S components) recorded at these four GNSS stations are illustrated below. The plate motion was previously removed using the data recorded at REUN IGS station located 15 km away from the volcano (Figure 1 inset) where only plate motion is considered to take place. The InSAR-derived time series of ground displacements (U-D and E-W components) at GNSS stations are superimposed on the GNSS measurements (red circles). The GNSS time series are shifted to reduce the root mean square (RMS) discrepancies between InSAR and GNSS observations to make different geodetic data comparable in time. Blue areas and dark blue vertical lines represent timing of eruptions and intrusions, respectively, which occurred in the study period (

Extracting Single Lava Fields between March 1998 and April 2007
At Piton de la Fournaise, lava flows are often superimposed on those emitted during earlier eruptions ( Figure 1). Such superimposition of two or multiple lava flows can induce more complex processes, such as compaction in the rough interface between two flows, subsequent thermal contraction of the whole composite lava, and the acceleration of creep in the substrate resulting from the heat transfer from the above to the underlying flow [2]. In this study, we therefore only consider the displacement restricted to single lava flow fields (Figure 3), attempting to avoid the considerations of interactions between lava flows. The selection of lava flow fields followed three principles:

Estimating Lava Thickness
To evaluate the thickness of lava flows emplaced between 1998 and 2007 at Piton de la Fournaise, we used two available Digital Elevation Models (DEMs). One is a 25 m resolution DEM produced by IGN from aerial photogrammetry in 1997 and the other is measured via airborne LiDAR surveying in 2008 with a 5 m resolution. It is worth noting that no eruptive events generating surface displacement occurred between the acquisition time of the 1997 DEM and March 1998, and between April 2007 and the acquisition time of the 2008 DEM. The two DEMs were first downsampled to the same grid posting. We then used the vertical difference between the two DEMs to estimate topographic height change, which reflects the thickness of total deposits. Given that lava flow is the main product of Hawaiian eruptions and composes much of the deposits, as at Piton de la Fournaise, and no independent data providing information about the discrepancies between the total deposit thickness and the lava thickness is available, we used the difference between the two DEMs to approximately represent the thickness of the lava flows. It is reasonable to some extent according to the previous study carried out by Chaussard [16], who used the height change to represent the lava thickness at the Paricutin volcano, even though ash and pyroclastic are expected to account for up to 30% of the total deposit thickness. Additionally, the effects of erosion and redistribution of the less consolidated materials were not considered in this study. We used the extracted lava flow fields ( Figure 3) as a mask to obtain the thickness of each lava field. The uncertainty of thickness of the lava flows was estimated to be 0.97 m by calculating the average standard deviation of the thickness in the areas without any emplacements of lava flows.

Results
Griffiths & Fink [31] and Griffiths [32] showed that after the emplacement of lava flows, a solid crust quickly forms and grows at their surface. The flow spreads only in depth from a certain thickness. The strength in the plane parallel to lava flow surface is far larger than that in the plane perpendicular to its surface. Thus, the displacement due to post-lava emplacement processes should be largely vertical and horizontally negligible (except for some local areas with steep slope).  (Figures 1 and 2). It suggests that the horizontal displacement induced by post-emplacement of the lava fields is negligible. Since an evident regional difference in displacement arises, we divide the lava fields into two groups (inside and outside the EFA) for further analysis.  (Figure 5d), which is about twice larger than the magnitude of subsidence occurring in the 10 years older lava flow fields (the Mar 1988a-b lava fields; Figure 5a,b). The Mar 1998a lava field (Figure 5a) exhibits an absolute subsidence-thickness slope of 0.0014, approximate to that of the Mar 1998b lava field that emplaced in the same year (0.0013; Figure 5b) and a bit smaller than that of the Feb 2000 lava field that emplaced two years later (0.0016; Figure 5c).
To investigate the temporal behavior of subsidence, we selected one point of 10 m thickness in each lava field outside the EFA (see Figure 4 for location of points) as representatives and plot the displacement time series in Figure 6a. The four time series display a similar trend, revealing the four lava fields subside following a similar temporal behavior between 2009 and 2014-continuous subsidence without being clearly interrupted by volcano eruptive phases and a subsidence rate decaying with time. For the lava field that emplaced in Feb 2000 (green triangles in Figure 6a), the magnitude of subsidence rate (at the 10 m thick point) decays from~0.04 m/year at the beginning of study period (nine years after it emplaced) to~0.01 m/year about five years later, a reduction of three quarters.   It is worth noting that at a certain point of time, the subsidence-age relation in Figure 6 cannot signify the subsidence-age relation of lava fields in general since the plot is based on time series observed at one point for each lava field.

Subsidence in Lava Fields inside the EFA
The same analysis was performed on the lava fields inside the EFA (result figures: Figures 6b, 7 and 8). Similar to the lava fields located outside the EFA, the inside lava fields show a trend of downward displacement all over the area during the study period (Figures 7 and 8). However, the subsidence pattern and behavior in the EFA display huge variability contrasts with the uniformity shown outside the EFA. First, the magnitude of subsidence is generally larger in lava fields inside the EFA than those outside. The majority of pixels falls into the cumulative subsidence interval 0.1-0.2 m, regardless of lava thickness and age (Figures 7 and 8).  (Figures 7 and 8). Second, the subsidence-thickness relationship behaves differently among the five lava fields. The June 2000 and Oct 2000 lava fields will be two remarkable exceptions since their subsidence magnitudes decrease with the increase of lava thickness (Figure 8a,b). They both display a linear trend and the decreasing rate of subsidence with lava thickness (the subsidence-thickness slope) of the June 2000 lava field (0.0042; Figure 8a), however, is significantly larger than that of the Oct 2000 lava field (0.0007; Figure 8b). Third, the subsidence-thickness slope and the age do not show any clear correlation. Indeed, the absolute subsidence-thickness slope of the Nov 2002 lava field (0.0052; Figure 8d) is larger than that of a one year older lava field (June 2001, 0.0034; Figure 8c). Nevertheless, the youngest lava field (Aug 2004) has the lowest slope among these three (Figure 8c-e). Furthermore, the June 2001 and Nov 2002 lava fields subside (m/m of thickness) about twice and three times more significantly than the lava that emplaced just one and two years ago in Feb 2000 outside the EFA (Figures 5c and 8c,d). More remarkably, the absolute subsidence-thickness slope of the June 2001 and Nov 2002 lava fields (0.0034 and 0.0052; Figure 8c,d) are even greater than the six years younger lava field that emplaced outside the EFA (Figure 5d). Lastly, relatively large norms of residuals for the June 2001 and Nov 2002 lava fields indicate that linear fits seem unable to satisfactorily explain the relationship between their subsidence and thickness.   Despite all the differences, the time series of lava fields inside the EFA show a similar displacement behavior in time (continuous subsidence with rates decreasing in time) to those lava fields outside (Figure 6a,b). For a lava field that emplaced in June 2000 (blue squares in Figure 6b) is almost the same time as the Feb 2000 lava emplacement, however, the magnitude of the subsidence rate is more important (up to~0.07 m/year at the beginning of study period and~0.02 m/year at the end). As already mentioned, the subsidence at a certain time does not exhibit clear correlation with the lava age.
The complexities of displacement in the lava fields located inside the EFA are evident when compared with the lava subsidence outside the EFA, indicating that the post-emplacement subsidence of lava fields may not be the single component contributing to the downward displacement observed in the EFA.

Only Post-Lava Emplacement Processes Account for the Subsidence in the EFA?
The vertical displacement (2009-2014) affecting the lava fields that emplaced between 1998 and 2007 at Piton de la Fournaise has been derived and analyzed in this study. Results show that the cumulative subsidence in lava fields outside the EFA increase with the increase of lava thickness, presenting a positive linear correlation (Figure 7). Similar correlations have been previously observed from lava fields at other volcanoes, such as the Etna, Okmok, Santiaguito, and Paricutin volcanoes [2,16,18,19]. For every meter of lava, the magnitude of subsidence increases as the lava age decrease. For certain lava fields, the subsidence rate decays in time (Figures 5 and 6a). These observed subsidence behaviors are expected to be the results of post-lava emplacement processes (i.e., the ongoing thermal contraction and compaction induced by the cooling and solidification of lava fields and/or the substrate [2]). However, the subsidence behaviors in lava fields inside the EFA differ significantly from those outside the EFA and in previous studies (Figures 6b, 7 and 8). The correlations between subsidence and thickness vary from one to another. Subsidence behaviors in the June and Oct 2000 lava fields even display negative correlations with thicknesses. The subsidence per meter of thickness does not correlate with the lava age ( Figure 8). This irregularity is further illustrated in Figure 9, in which an exponential curve is fitted using values observed in lava fields outside the EFA to roughly explain the expected correlation between the subsidence (per meter of thickness) induced by post-lava emplacement processes and the lava age. The markers corresponding to lava fields inside the EFA deviate broadly (

Contribution of the Post-Lava Emplacement Processes to the Subsidence in the EFA
The portion of subsidence (per meter of thickness) due to post-lava emplacement processes in lava fields inside the EFA is expected to drop on the fitted curve (as shown by the transparent blue squares with approximated subsidence-thickness slopes in brackets in Figure 9). The discrepancies between the observed (blue squares) and the approximated (transparent blue squares) values reveal that the effects produced by other potential processes on lava fields inside the EFA vary spatially from north to south. The significant discrepancies of the June 2000, Oct 2000, June 2001, and Nov 2002 lava fields suggest that the potential processes might have accounted for a larger portion of the observed signals rather than post-lava emplacement processes. The two lava fields of 2000 located in the south, especially, display inversed subsidence-thickness slopes (subsidence decreasing with an increase of thickness) with respect to others. This suggests the potential processes could have played a dominant role in the pattern of the observed displacement affecting the 2000 lava fields. In particular, the areas covered by relatively thinner lava show a higher sensitivity in the accumulation of displacement in response to such potential processes. The Aug 2004 lava field located in the north, however, presents less discrepancy (close to the fitted curve; Figure 9) and reveals a relatively slighter effect from the potential processes. Based on this exponential law, we estimated that the contribution of post-emplacement lava flow subsidence to the subsidence observed inside the EFA could range from~22% (in the June 2000 lava field) to~79% (in the Aug 2004 lava field). In other words, the other processes could contribute up to~78% of the subsidence. In comparison with the rough estimate of 34% contributed by other processes derived by Chen et al. [23], the estimate of~78% from this study should draw more attention, since this large portion of flank motion triggered by unknown processes could be hazardous.

New Insights for Causes of Displacement on the Eastern Flank
As previously discussed in the work of Peltier et al. [7], Froger et al. [13] and Chen et al. [23], several possible origins, such as pre-existing structural discontinuity, summit stress associated with magmatism, and activated fault movement, have been proposed to explain the Eastern Flank motion of Piton de la Fournaise. However, the precise origin remains undetermined. According to the spatial variation of the displacement of lava fields revealed in this study and the investigation of previous studies, it is suggested that some potential structures (e.g., normal faults) could exist. Such potential structures (dip and strike) could make an important contribution to influencing the spatial distribution of strain accommodation. A structural feature (F1, as indicated in Figures 2 and 3) has been previously identified by early researches [13,33]. The significant difference in observed displacements between the Aug 2004 lava field and the other four in the EFA (Figure 9) also suggests an inclined structure/fault could underlie between the Aug 2004 and the June 2001 lava fields. Another structure located along the deformed boundary (F2 in Figures 2 and 3) was suggested by Froger et al. [14] that it could be related to a normal fault activated during the March-April 2007 eruption. It is also supported by the observations of this study by taking into account the red rectangle displacement area indicated by F2 in Figures 2 and 3. The two structures make the area covered by the 2000-2002 lava fields a mobile hanging wall, which indicates that the internal rocks (as well as their overlying lava fields) between the two structures would significantly slide along a plane beneath the Eastern Flank, contrasted to the external rocks (Aug 2004 lava field). This interpretation could explain the significant discrepancies of the June 2000, Oct 2000, June 2001, and Nov 2002 lava fields and the small discrepancy of the Aug 2004 lava fields in Figure 9. However, the existence of the interface beneath the Eastern Flank is still a subject of debate, which is mainly due to the rareness of clear shallow seismicity below the Eastern Flank and that the earthquakes below the mobile flank are not aligned along a plane [7,34,35]. Thus, according to available data and the results retrieved in this study, the deep mechanism driving the flank motion at Piton de la Fournaise is always a valuable issue to investigate.

Conclusions
We measured the long-term displacement in recent lava fields at Piton de la Fournaise volcano using TSX/TDX InSAR data between 2009 and 2014. The spatial pattern and temporal behavior of the displacement of each lava flow, and the relationship between displacement, lava thickness and age are investigated by a detailed statistical analysis. The subsidence behavior exhibits a significant difference between the lava fields outside and inside of the EFA (Eastern Flank Area). The former shows a clear linear relationship with lava thickness and it decays with lava age, which is attributed to the post-lava emplacement processes. The latter however displays an irregularity in the subsidence and the relationship between subsidence and thickness and age, which confirms that, in addition to the post-lava emplacement processes, other processes must have played a role in the observed subsidence on the Eastern Flank. Qualitative analysis show that the contribution of other processes to the flank downward motion could reach up to 78%. A detailed analysis of the displacement behavior in lava fields inside the EFA favors the interpretation that two pre-existing normal faults could have been activated during the March-April 2007 eruption and have resulted in a slip along a pre-existing structural discontinuity beneath the Eastern Flank, though the existence of the discontinuity is still an issue of debate. This paper presents for the first time a spatial and temporal analysis of ground displacement in multiple lava fields at Piton de la Fournaise, by characterizing the correlations between lava subsidence, lava thickness, and age. Although more studies are required to determine the dynamic process driving the instability of the Eastern Flank, the results from this study add new knowledge and essential insight for understanding and monitoring potentially hazardous flank instabilities.