Monitoring an effusive eruption at Piton de la Fournaise using radar and thermal infrared remote sensing data: insights into the October 2010 eruption and its lava flows

Abstract Accurate and fast delivery of information about recent lava flows is important for near-real-time monitoring of eruptions. Here, we have characterized the October 2010 lava flow at Piton de la Fournaise using various InSAR datasets. We first produced a map of the area covered by the lava flow (i.e. Arealava=0.71–0.75 km2) using the coherence of two syn-eruptive interferograms. Then we analysed two post-eruptive InSAR datasets (i.e. monostatic and bistatic data). The monostatic database provided us simultaneously with the displacement rates, lava thickness, volume and volume flux. We found that the lava flow was subsiding and moving eastward at maximum rates of 13±0.3 and 4±0.2 cm a−1, respectively. Also, it had a mean thickness of Zmean=5.85 m, VolDRE=1.77±0.75×106 m3 (1σ) and MOR=1.25±0.53 m3 s−1. The bistatic database provided us only with the thickness and volume information (i.e. Zmean=6.00 m, VolDRE=1.83±0.65×106 m3 and MOR=1.29±0.46 m3 s−1). Finally, we used a thermal remote sensing technique to verify the InSAR-derived measurements. Results show that the monostatic and bistatic datasets were both well within the range for the DRE volume obtained from MODIS data (2.44–4.40×106 m3). Supplementary material: Tables A1 and A2 give satellite images used in this study. Table A3 gives the parameters used for the calculation of the effusion rates. The figures give the data processing of the post-eruptive radar images. These are available at https://doi.org/10.6084/m9.figshare.c.2213563

One of the main hazards posed by volcanoes that erupt effusively is lava flow propagation in inhabited areas. Although most lava flows move slowly enough to allow people to be evacuated, they can still destroy property and can impact hugely on the economy. Performing near-real-time monitoring of an advancing lava flow involves a high level of risk and is not always easy, especially if the eruption occurs in a remote area or in bad weather conditions. Many volcano observatories dedicate much effort to predicting the most likely future courses of lava flows to aid civil protection groups in their decision making. As a consequence, several algorithms for simulating lava flow paths have been developed in the last decade. An accurate map of recently emplaced lava flows is obviously a key input in the a posteriori validation of these algorithms. Unsurprisingly, the reliability of the simulations depends heavily on the quality of the digital elevation model (DEM) used as a basis for the numerical propagation of the lava flow. In the case of frequently erupting effusive volcanoes, maintaining an up-to-date and accurate DEM can be a real challenge, since each new lava flow changes the topography. However, if an accurate thickness map of the lava flow can be produced soon after its emplacement, the DEM can be easily updated by simply adding this thickness map. Also, the thickness map can be used to estimate the volume and effusion rate, which provide relevant insights into the short-and long-term behaviour of the volcano, the mechanism of effusive eruptions and the dynamics of the magmatic system feeding the eruption (e.g. Wadge et al. 1975Wadge et al. , 2006Wadge 1981;Wadge & Guest 1981;Dvorak & Dzurisin 1993;Denlinger 1997;Harris et al. 2000Harris et al. , 2005Harris et al. , 2011. Another issue apart from the planimetric and volumetric characterization of lava flows is the ground surface displacement associated with a lava flow. These displacements affect not only the lava flow itself (i.e. thermomechanical compaction) but also the surrounding substratum (Murray 1988;Briole et al. 1997;Stevens et al. 2001). While they potentially provide relevant information about lava flow thickness and cooling processes or about substratum rheology, they are sometimes considered as unwanted signals if the main source of the displacements is related to internal processes such as pressure variation in a deep magmatic reservoir or magma transport within the volcano plumbing system.
In order to characterize the area, thickness, volume and volume flux of ongoing or recently emplaced lava flows, most volcano observatories use geodetic and remote sensing ground-based and airborne approaches, where techniques based on GPS, hand-held thermal imagery, photogrammetry and LiDAR are available (e.g. Harris et al. 1997Harris et al. , 2005Stevens et al. 1997;Mazzarini et al. 2005;Coltelli et al. 2007;Neri et al. 2008;Favalli et al. 2010). However, such techniques often face problems related to: (a) the cost and duration of the mission; (b) accessibility of the target area and cloud cover; (c) safety during data acquisition; and most importantly, (d) the accuracy of the acquired measurements. In this respect, the classic approach consisting of walking around or flying over the newly emplaced flow with a GPS to map the flow area and then calculating a volume assuming a constant thickness across the entire lava flow (generally the edge thickness) provides a very rough estimation, with uncertainties of up to c. 200% (e.g. Stevens et al. 1997;Harris & Neri 2002). Nevertheless, if an emplaced volume is obtained, the value can be divided by the duration of emplacement to obtain a time-averaged discharge rate (TADR), which becomes a mean output rate (MOR) when the total erupted volume is divided by the entire eruption duration from start to finish (Harris et al. 2007).
Following the launch of the first European radar satellites, ERS 1 and 2, in the 1990s, interferometric synthetic aperture radar (InSAR) became a widely used technique for deriving ground surface displacements, no longer requiring the deployment and maintenance of ground-based instruments as is the case for in-situ geodetic methods. In this work, we show that InSAR can be used not only to image the syn-eruptive ground surface displacement related to the dyke intrusion that fed the October 2010 eruption at Piton de la Fournaise, but also to obtain an early map of the ongoing lava flow during the eruption and to estimate, with a high level of accuracy, its surface area, thickness and volume. We also characterize the displacement behaviour of the lava flow in the months following its emplacement.
In InSAR processing, a DEM is normally used to correct for topography in the interferograms. Any topographic phase shifts, corresponding to changes in topography between the time of the DEM production and the interferograms, are considered as artefacts. Here, we exploit this idea to derive the morphology and volume of the October 2010 lava flow using two types of high spatial resolution X-band datasets: (a) monostatic Cosmo-SkyMed; and (b) bistatic TanDEM-X. These two databases are quite different, so we also compared them in terms of what information each can best deliver.
For further validation of our technique, we crossvalidated the InSAR-derived volume and volume flux with those obtained from satellite thermal data. Finally, we discuss the applicability of this study to near-real time monitoring of effusive eruptions.

Piton de la Fournaise and the volcano observatory
Piton de la Fournaise is one of the most active basaltic shield volcanoes in the world. It forms the southern half of La Réunion Island, an oceanic island in the south of the Mascarene Basin (Indian Ocean), 800 km to the east of Madagascar (Fig. 1). It is the most recent manifestation of the hotspot that created the Cretaceous Deccan Traps (e.g. Courtillot et al. 1986;Duncan & Pyle 1988). Its morphology is characterized by a series of large depressions bounded by a steep horseshoe-shaped scar opening to the east (e.g. Bachèlery & Mairine 1990;Lénat & Labazuy 1990;Lénat et al. 2001;Bachèlery & Michon 2010). Most of its historical eruptions have occurred within the Enclos Fouqué -Grand Brûlé structure, along the two zones of preferential fissures oriented at N108 E and N1708 E (e.g. Bachèlery 1981;Lénat et al. 2011). The Central Cone is found at the junction of these two zones and has an east -west elongated shape, with two summit pit-craters -the Bory crater to the west and the Dolomieu crater to the east (Bachèlery 1981;Lénat et al. 2000).
The Enclos Fouqué -Grand Brûlé structure is an unpopulated area bounded by steep cliffs to the north, south and west. These act as barriers, trapping lava flows within the structure. However, there are times where eruptions also happen outside the Enclos Fouqué -Grand Brûlé structure (Peltier et al. 2009b). These can threaten populated zones, as was the case during the 1977 eruption when lava flows reached the village of Piton Sainte Rose on the eastern shore of the island and devastated the village.

InSAR volcano monitoring of Piton de la Fournaise
Piton de la Fournaise is one of the few volcanoes in the world that is monitored on a regular basis from space by the use of InSAR data. Thirty of the 39 eruptions that occurred between 1998 and 2010 were imaged using data provided by the Canadian RADARSAT-1 satellite (Sigmundsson et al. 1999;Fukushima et al. 2005Fukushima et al. , 2010, the European ASAR-ENVISAT satellite (Froger et al. 2004(Froger et al. , 2015aTinard 2007), the Japanese ALOS-PALSAR satellite (Augier 2011), the German TerraSAR-X satellite (Froger et al. 2012) and the Italian Cosmo-SkyMed satellite (Bato 2012). InSAR monitoring of Piton de la Fournaise became routine in 2003, initially through the Observatoire de Physique du Globe at Clermont-Ferrand, France. In 2010, this monitoring was integrated into the National Service for Volcanological Observations of the French National Institute for Earth Science and Astronomy under the name OI 2 (Observatoire InSAR de l'Ocean Indien). Currently, the OI 2 database (https://wwwobs.univ-bpclermont.fr/casoar) includes around 900 radar images of Piton de la Fournaise acquired by various radar satellites from February 2003.
The resulting InSAR database represents an exceptional record of the ground surface displacements that have occurred at Piton de la Fournaise since 2003. It has revealed not only the displacements that are directly induced by magma transfer within the volcano (Sigmundsson et al. 1999;Froger et al. 2004;Fukushima et al. 2005Fukushima et al. , 2010, but also the possible role of the hydrothermal system as a source of the deformation as well as the existence of a large active detachment affecting the eastern flank of Piton de la Fournaise (Froger et al. 2015a, b).

The October 2010 eruption
The October 2010 eruption began on 14 October around 7:30 a.m. (local time) in the vicinity of the Château-Fort crater, at the southeastern base of the Central Cone, about 1.5 km SE of the Dolomieu crater (Fig. 1). The beginning of the eruption was marked by lava fountains that created four distinct eruptive cones. By the third day, only one cone remained active. It continued to feed a lava flow that slowly moved southeastwards until the eruption ceased at 2:00 a.m. on 31 October, after 16.4 days of effusive activity (Roult et al. 2012). Petrological and geochemical analyses of samples showed that the lava flows from this eruption were aphyric basalts with less than 7% MgO (Roult et al. 2012). This is the case for many eruptions at Piton de la Fournaise, but more specifically for the eruptions after 2007 (Michon et al. 2009;Staudacher et al. 2009;Peltier et al. 2009a, b). Soon after the end of lava flow emplacement, a GPS campaign was carried out to delineate the flow field margin. From these data, the surface of the lava flow was estimated at 0.68 km 2 and then, assuming a mean flow thickness of c. 4 m, its volume was estimated at 2.7 × 10 6 m 3 (Roult et al. 2012). These measurements also allowed a MOR of about 1.9 m 3 s 21 to be estimated (Roult et al. 2012).

The InSAR technique
The classic InSAR technique calculates the phase difference, Dw, between two SAR images acquired for the same area at different times (Massonnet & Feigl 1998). The result of the phase difference is the creation of a new image, an interferogram, where the phase depends on the changes between the two acquisitions both in the two-way propagation time of the radar wave between the ground and the radar antenna and in the interaction between the radar wave and the multiple backscatterer sources present on the ground. The phase changes owing to the radar wave propagation time can result from three different sources: (1) changes in the radar antenna position that induce phase changes commonly separated into a flat earth term (Dw orb ) and a topographic term (Dw GST ),; (2) changes in the ground surface position, that is, ground surface displacement (Dw GSD ); and (3) changes in the radar wave velocity related to changes in the properties of the atmosphere between the two SAR image acquisitions (Dw atm ). The changes in ground surface backscattering properties are generally considered as a source of noise in the interferometric signal, along with other sources of noise (1) so that each interferogram is described as the sum of the five terms mentioned above: Provided that the geometric and physical properties of the ground surface backscatterers do not change between the two acquisitions, the last term of equation (1) can be ignored. Then, depending on the objective of the study, the term of interest to be modelled is selected, and most of the other terms are removed in order to isolate and characterize it as accurately as possible. Most often in InSAR studies the selected term is the ground surface displacement.
The first term in equation (1) can be modelled and removed using precise orbital state vectors, which are provided by space agencies. The second term can be removed using a DEM. Since the DEM and the actual topography of the study area can differ from each other, residual topographic phases on the interferograms can remain, the magnitude of which depends on the horizontal orbital separation of the satellite positions at the time of image acquisition. When projected in the slant range, this horizontal separation is known as the perpendicular baseline, B perp . It has an inverse relation to the altitude of ambiguity that is a measure of the change in elevation required to generate a 2p difference (equivalent to one cycle) in the phase after the interferogram has been corrected for the curvature of the Earth. A low altitude of ambiguity means that there is a high sensitivity in the interferogram to topographic change.
If the geometric or physical properties of the ground surface backscatterers change significantly between the two radar acquisitions, the last term of equation (1) dominates the interferometric signal, which then becomes very noisy, making it hard to characterize the other terms properly. Many factors can contribute to a change in the ground surface backscatterers. These include growth of vegetation, surface erosion, snowfall, strong ground surface displacement at sub-pixel scale and, in our case, emplacement of a volcanic deposit. Coherence is an inverse function of the local variance of the interferometric phases and is typically used as a means of estimating the temporal stability of surface backscatterers, which can be exploited to assess the quality of an interferogram (Massonnet & Feigl 1998).
Generally, during InSAR monitoring of an effusive crisis, two signatures are expected to be observed: † The interferograms that span the early stage of the eruption may have recorded the ground surface displacement related to the intrusion of the dyke that fed the eruption. These displacements can be modelled to infer the dyke geometry and location (Sigmundsson et al. 1999;Massonnet & Sigmundsson 2000;Amelung & Day 2002;Froger et al. 2004;Fukushima et al. 2005Fukushima et al. , 2010. † The interferograms that span entirely, or partly, the period of the lava flow emplacement are expected to be incoherent in the area covered by the lava flow because the backscattering properties of the ground after the emplacement will change radically (Zebker et al. 1996Dietterich et al. 2012). Interferometric coherence maps can be calculated at this stage and used to carry out an early mapping of the emplaced lava flows. The interferograms that span the period after the end of lava flow emplacement are expected to become increasingly coherent as the lava flow surfaces become progressively stable, so that we expect two further, essentially different, signals: † Topographic residues that correspond to the changes between the original DEM and the new topography of the volcano as modified by the emplacement of the lava flow field (Rowland et al. 1999(Rowland et al. , 2003Stevens et al. 1999). † Displacement signals that are related to the thermal compaction of the lava flow and the visco-elastic flexion of the substratum owing to the loading of the lava flow itself (Murray 1988;Briole et al. 1997;Stevens et al. 2001;Lu et al. 2005).
The thermal remote sensing method Pieri & Baloga (1986) proposed the following steady-state relationship to relate active lava flow area to volume flux: where Q is lava volume flux, 1 is the emissivity, s is the Stefan-Boltzmann constant, T e is flow surface effective radiation temperature, r is flow density, C p is the lava specific heat capacity, T o is the eruption temperature, T f is the temperature at which forward motion ceases, A is the flow area and L is the maximum possible distance from the vent reached by the lava fed at volume flux Q. This relationship was later applied by Crisp & Baloga (1990) to planetary lava flows and was first adapted by Harris et al. (1997) to satellite data to convert spectral radiances to TADR. Wright et al. (2001), Harris et al. (2007), Harris & Baloga (2009) and Harris (2013) all point out that, when applied to satellite thermal data, the conversion reduces to a linear relation between spectral radiance, lava flow area and TADR, when the flow is assumed to have reached steady state (Garel et al. 2012).

DEM
In 2008, the French Geographic Institute (IGN) carried out two airborne LiDAR surveys over La Réunion Island. They used this data to produce a 1 m resolution DEM along a 1 km-wide coastal strip, and a 5 m resolution DEM of the entire island. Through comparison with 93 precise groundreference points, taken from the geodetic network of the Observatoire du Piton de la Fournaise (i.e. the permanent GNSS station and benchmarks included in the non-permanent network), the vertical precision of the IGN LiDAR DEM is estimated to be 1.6 m (mean deviation between DEM and GNSS elevations). With respect to the previous La Réunion IGN DEM, produced from aerial stereophotographs taken in 1997, the LiDAR DEM provides exceptionally high accuracy, with a spatial resolution improved by a factor of 5. Additionally, the LiDAR DEM takes into account most of the changes that occurred between 1998 and 2008 owing to lava flow emplacement, cinder cone construction and caldera collapse.

InSAR data
We used three different types of X-band (l ≈ 3.1 cm, f ¼ 9.6 GHz) radar data for the InSAR characterization of the October 2010 eruption. These were: (a) TerraSAR-X monostatic data; (b) TerraSAR-X -TanDEM-X bistatic data; and (c) COSMO-SkyMed data.
All of the interferograms used in this study were produced using the DIAPASON software (CNES 1996) that applies the two-pass method described by Massonnet & Feigl (1998). The contribution of the orbital trajectories in the interferograms was modelled and removed using the orbit state vectors provided by DLR (German space agency) and ASI (Italian space agency). We used the IGN 5 m LiDAR DEM to model and remove the topographic phases from the interferograms and to provide a precise geographic frame (i.e. UTM-WGS84). As the DEM dates from 2008, any topographic changes that occurred after 2008 will produce residual topographic fringes in the interferograms.
TerraSAR-X (TSX) monostatic data. TerraSAR-X (TSX) is a German Earth-observation satellite launched on June 2007. Its primary payload is an X-band radar sensor operating at different look angles (20 -608), resolutions and polarizations. The TSX active radar allows data acquisition with a return period of 11 days, both day and night, regardless of the weather conditions. In this study, we used two ascending and two descending TSX monostatic Stripmap images (resolution c. 3 m) that spanned from 1 September to 26 October 2010. From these we calculated two syn-eruptive interferograms.
TerraSAR-X -TanDEM-X bistatic data. On 21 June 2010, TanDEM-X (TDX), a twin satellite of TSX, was launched by DLR. While radar data can still be acquired by TDX in a classic monostatic mode similar to that of the TSX (i.e. one single satellite successively emits and receives a radar pulse), the real innovative characteristic of the TanDEM-X mission (TDM) is its bistatic acquisition mode where the two satellites are flying in closely controlled formation, separated by distances of between 250 and 500 m. In the bistatic configuration one of the satellites emits the radar pulse and the radar echo backscattered from the ground is received and recorded simultaneously by the two satellites from slightly shifted positions, resulting in two radar images. There are two main consequences of this particular mode of acquisition. Firstly, since the two satellites receive the radar pulse simultaneously, they do not record ground surface displacement signals. Secondly, there is no difference in the atmospheric phase delay between the two recorded radar images, since the backscattered radar pulse received by TSX and TDX is passing through the same atmosphere. Furthermore, it is also worth noting that, in this mode of acquisition, the coherence of TDM interferograms is expected to be very high. As a consequence, the differential phase between the two simultaneous bistatic images is only due to the difference in range between the ground surface and the antenna of the satellites, hence the bistatic interferograms can be used to determine accurately the ground surface topography. In fact, the main objective of the TDM mission is to produce high-resolution DEMs globally. If a DEM is used in the bistatic interferogram processing, in the same way as for monostatic interferogram processing, the residual phases will record only the changes in topography that occurred between the time of the DEM production and the bistatic image acquisition. This approach obviously has great potential for mapping newly emplaced volcanic products.
We used 41 TSX-TDX pairs of Stripmap bistatic images that were acquired from both ascending and descending orbital passes between 6 September 2011 and 12 November 2012 (i.e. about one year after the end of the October 2010 eruption). From the 41 image pairs we produced 41 bistatic interferograms. Note that in this paper, we refer to the satellites as TSX or TDX, and to the TSX-TDX bistatic mission as TDM. Hence, TDX data refers to radar images acquired by TDX alone and TDM data are the pairs of TSX and TDX radar images acquired bistatically.
COSMO-SkyMed data. COSMO-SkyMed, Constellation of small Satellites for the Mediterranean basin Observation, is an Italian Earth-observation constellation of four X-band radar satellites. Each of the satellites has a return period of 16 days. By using different satellites in the constellation, the routine return period can be reduced to 8 days. In this study, we used 54 Stripmap images acquired during ascending and descending orbital passes between 15 February 2011 and 29 December 2013. From these images we computed 412 interferograms.
Thermal remote sensing data MODIS data. The Moderate Resolution Imaging Spectroradiometer (MODIS) is a NASA sensor dedicated to observing the global dynamics of the Earth's surface and the lower atmosphere. It is flown on two platforms, the Terra platform and the Aqua platform. For our purposes, the 1 km data acquired in the mid-infrared (bands 21 and 22, both 3.929-3.989 mm) and thermal infrared (band 32, 11.770-12.270 mm) are the most useful. We analysed 44 MODIS images acquired during the eruption: 24 were from MODIS/Terra and 20 from MODIS/Aqua (up to four images per day during the October 2010 eruption).

Syn-eruptive displacement maps
The TSX syn-eruptive interferograms gave a very clear image of the October 2010 eruption (Fig. 2a,  b). The ascending interferogram (Fig. 2a) spans a 33-day period between 19 September and 22 October 2010. It exhibits two main fringe patterns. The westernmost pattern, elongated in a NNW-SSE direction, joins Dolomieu crater to the location of the October eruptive fissures. It corresponds to a decrease in the satellite-ground distance between the two acquisitions by about 13.5 cm in the line-of-sight (LOS) direction. The easternmost pattern is elongated in an NNE -SSW direction, implying an increased distance in the LOS of about 7.5 cm. Two minor fringe patterns are also visible just to the west of the Château Fort crater indicating 3 and 21.5 cm displacement in the LOS direction for the upper and lower set of fringes, respectively.
The descending interferogram (Fig. 2b) spans a 55 day period between 1 September and 26 October 2010. It shows two distinct fringe patterns. The main pattern indicates a decrease in the satellite -ground distance of about 42 cm in the LOS direction. The second fringe pattern to the SW of the Château Fort crater probably has several fringes, but only two are discernible in the increasing LOS view, revealing an increase of at least 3 cm in the LOS direction. Figure 2c and d illustrates the east-west (EW) and vertical component maps of the syn-eruptive displacements as derived from the two interferograms following the method of Wright et al. (2004). Here, we assumed that the northsouth (NS) displacements vanished when projected along the LOS. The horizontal component (Fig.  2c) indicates an opening of up to 40 cm which is mostly eastward, and the vertical component (Fig. 2d) shows a general inflation of up to 30 cm. This presence of both opening and uplift is a typical displacement signature of a dyke forcing itself towards the Earth's surface. Assuming an elastic host rock with a Poisson's ratio of 0.25, and a ratio of 3:4 between the uplift volume and the injection volume (i.e. typical of a dyke, Delaney & McTigue 1994), the integrated ground surface uplift yields a minimum volume of injected magma of 0.8 × 10 6 m 3 .

Syn-eruptive lava flow maps
On Figure 3 we show the coherence maps for the two TSX syn-eruptive interferograms. While the flanks of Piton de la Fournaise are incoherent owing to the presence of dense tropical vegetation, the upper part of the Enclos Fouqué -Grand Brûlé structure exhibits relatively good coherence. In fact, since Piton de la Fournaise is a very active volcano, the ground surface of the Enclos Fouqué -Grand Brûlé structure, where most of the historic eruptions occurred, is frequently renewed by the emplacement of new lava flows, hence there is not enough time for vegetation to develop. In contrast, the interior of the summit crater is characterized by a very low coherence. This is mainly due to recurrent rock falls from the unstable walls of the crater and associated debris deposits since the collapse of 2007. The northeastern region also exhibits incoherent pixels that are related to the dense vegetation present in this area.
The October 2010 lava flow is characterized by a strong contrast in coherence with its environment. The radical change in the backscattering properties of the ground after the emplacement of the flow induces a total loss of coherence for all the area covered by the lava flow. The boundaries of the incoherent area can be traced in order to derive the extent and shape of the lava flow field, providing an estimate of the area it covered. From the ascending coherence map (Fig. 3a) we obtained an area of 0.71 km 2 , compared with 0.75 km 2 for the descending map (Fig. 3b). The variation in the estimated values is due to the difference in time covered by the two coherence maps. While the ascending coherence map only covered the first 9 days of the eruption, the descending one covered the entire eruption, meaning that a further 0.04 km 2 of new ground was covered after the ninth day. These values are remarkably close to the 0.68 km 2 deduced from the GPS measurements of the lava flow contour. The slight difference could arise from two main reasons: † The fact that we may have incorporated some incoherent pixels within the calculated area that do not actually belong to the lava flow (e.g. pixels close to the cinder cones where pyroclastic fall deposits accumulated during the eruption induce loss of coherence). † The fact that the 5 m resolution of the coherence maps is less accurate than ground delineation by GPS in terms of lava flow-margin-delineation. We used the descending coherence map to produce a preliminary lava flow mask. This mask was used in the lava flow characterization procedure presented above.

Post-eruptive lava flow displacement maps, thickness maps and volume estimates
We present two InSAR volumetric characterization cases: one with a large database of monostatic CSK data and another with a small number of bistatic TDM data. Because each CSK monostatic interferogram has potentially recorded signals that are related to topographic residuals, ground surface displacements, atmospheric phase delays and other sources of noise, we applied a statistical approach to the dataset so as to reduce the contributions of any unwanted signals present in each interferogram. We intentionally used as complete a set as possible of monostatic interferograms in order to provide a large range of sensitivities to both the spatial and temporal baselines. Such an approach is obviously not compatible with near-real-time characterization of lava flows since it requires a large interferometric dataset, which takes time to establish.
As the TDM bistatic data are theoretically not affected by unwanted signals, they potentially allow for characterization of the lava flow soon after its emplacement. Therefore, in contrast to the first approach, and despite the fact that the first TDM bistatic data were acquired almost one year after the end of the October 2010 eruption, we intended to use as few interferograms as possible to characterize the lava flow, in order to test the near-real-time capability of this type of data.
The third case presented here uses the MODIS data for comparison with the volumes obtained from the InSAR-based approaches.
Case 1: using a large dataset of repeat-pass CSK data Our approach is not fundamentally different from that proposed by Ebmeier et al. (2012), who characterized the thickness and subsidence of the lava flows at Santiaguito using CSK InSAR data. However, in our work, we were able to use a larger interferogram dataset with three different acquisition geometries. This allowed us to retrieve the up-down (UD) and EW components of the displacement associated with the lava flow, whereas Ebmeier et al. (2012) were only able to retrieve the LOS displacement.
The processing begins by extracting a square region of interest, about 3.4 km wide, centred on the October 2010 lava flow field. We found two end members for the interferometric signals in our CSK dataset. Firstly, in interferograms with a low altitude of ambiguity and spanning a short duration, the area of the October 2010 lava flow is marked primarily by topographic signals. Secondly, in interferograms with a high value of altitude of ambiguity and covering a longer duration, the area of the October 2010 lava flow is marked mainly by ground surface displacement fringes. Most of the CSK interferograms are between these two types of end members. As a consequence we developed an approach to discriminate one from the other.
Given that interferograms record a phase signal that has a 2p modulus of an unknown absolute phase, we converted the ambiguous phase into an absolute phase, through phase unwrapping. Here we adopted the phase-unwrapping algorithm, SNA-PHU, developed by Chen & Zebker (2000. However, because the phase cannot always be unwrapped correctly, it can result in a poorly estimated value, mainly owing to the spatial aliasing of the phase. We encountered this difficulty in unwrapping several CSK interferograms because of the high spatial displacement rates that were present in the thickest part of the lava flow. After unwrapping, to reduce the possible contributions of large wavelength interferometric signals, such as atmospheric artefacts, residual orbital ramp and/or large wavelength ground surface displacements, we de-trended each interferogram using a six-degree regression polynomial surface (Bato 2012). This surface is assumed to account for the large wavelength phase variations as a function of elevation, latitude and longitude. For each interferogram, we obtained the parameters of the polynomial by calculating the regression on the whole region of interest, except for the area covered by the lava flow. The use of a high-degree polynomial was based both on the empirical observation and on the misfit value between the original phase and the simulated n-polynomial surface (i.e. n is from 1 to 10). Bato (2012) showed that, for the CSK-dataset, the misfit is at its minimum at the sixth-degree polynomial. Moreover, because badly unwrapped data are mostly observed at low coherence values, we also masked areas with a coherence of lower than 6.25%, an empirically estimated threshold, prior to de-trending.
The final step in the processing chain is phase inversion. This allows us to finally obtain the characteristic parameters of the lava flow (thickness plus horizontal and vertical displacement rates). We assume that, after de-trending, the phase of each interferogram can be reduced to: where Dw GST is the phase related to changes in the ground surface topography, Dw GSD is the phase related to ground surface displacement and 1 is the residue. In full, this becomes: in which n and i are the interferogram and pixel indexes, respectively. The first term of equation (4), on the right hand side, is the phase due to topography, in which h a n,i is the altitude of ambiguity and z i is the difference between the actual topographic elevation and the elevation of the DEM. In the case of the October 2010 lava flow, the difference reflects the thickness of the lava flow. The second term of equation (4) represents the phase owing to ground displacement between the two acquisitions. Here, Dt n reflects the duration of the interferograms in days, l is the wavelength of the satellite (for CSK, l ¼ 3.1 cm), X i is the displacement rate in the EW direction, Y i is the displacement rate in the UD direction, and R EWn,i and R UDn,i are the radar-looks, this being the angle from the horizontal plane in which the radar is pointing in the vertical and horizontal directions, respectively. The contribution of the NS direction was ignored in this study. Finally, the error resulting from disregarding this contribution can be described by the residual value, 1 n,i , which also includes other noise values related to failures in phase unwrapping or uncorrected atmospheric component.
By definition, the phase owing to topography and the inverse of the altitude of ambiguity are linearly correlated, as shown in equation (4). Displacements related to lava flow compaction and loading on its own substratum have been observed to decrease exponentially with time (Peck 1978;Stevens et al. 2001). However, in our case, empirical observations of the actual phase as a function of time and interferogram duration show that the phases within the area of the lava flow behave linearly, hence the decomposed ground surface displacement phase in equation (4). Here, we attribute the linearity to the short period of time covered by the whole dataset. The period of observation is probably too short to recognize the exponential behaviour. That is, the start of the dataset is 108 days after the eruption ended, and the whole dataset covers a total of 1048 days.
Since we are interested in deriving the characteristic behaviour of the October 2010 lava flow, we isolated the parameters on a pixel-by-pixel basis using a least squares solution following: where P i is the column vector of the parameter set that contains: (a) the thickness of the lava flow; (b) the EW displacement rate; (c) the vertical displacement rate; and (d) the residual noise, at the ith pixel. Parameter G is an n ×4 coefficient matrix, defined with respect to equation (4), V is the n × n full variance -covariance matrix of the unwrapped and de-trended interferograms and F i is an n×1 column vector that represents the unwrapped and de-trended phases at pixel i. V is calculated for the whole region of interest after masking the lava flow area. It is used to weight the inversion. The result of the inversion provided maps of each of the aforementioned parameters. In addition, we produced maps of parameter uncertainty, as estimated from the square root of the diagonal elements of the posterior covariance matrix (Strang 1986). The InSAR-derived thickness map (Fig. 4a) matches both the extent and shape of the October 2010 lava flow derived from the aerial survey conducted by OVPF shortly after the end of the eruption. It reveals maximum thickness values of 15.45 + 2.43 (1s) m in areas close to the vents. The mean lava flow thickness calculated using the preliminary lava flow mask is 5.80 m and the mean error (mean of the thickness standard deviation calculated for each pixel) is 0.71 m (1s). We then refined our estimation of the lava flow area using the following procedure. By dilating the preliminary lava flow mask with a one-pixel-width margin, while keeping all pixels within the new contour with a thickness greater than 20.71 m, we obtained a maximum lava flow area of 0.71 km 2 . Then, by retaining all pixels with a thickness greater than 0.71 m, we obtained a minimum lava flow area of 0.67 km 2 . We observed that the a posteriori mean error on the lava flow thickness using either the maximum or the minimum area is not significantly different from that obtained using the preliminary mask.
By summing the thickness values within the lava flow and multiplying by the spatial resolution of the map (in our case 5 × 5 m), we estimated the bulk volume to be 3.94 + 0.49 × 10 6 m 3 (1s). The uncertainty in the bulk volume was derived taking both the standard errors on the thickness and the uncertainty on the lava flow contour into consideration. If corrected for a vesicularity of 55.46 + 13.53% (1s), an estimation of the mean vesicularity for the recent lava flows on Piton de la Fournaise (Gurioli pers. comm., 2015), then the dense rock equivalent (DRE) volume is 1.75 + 0.75 × 10 6 m 3 (1s).
The EW displacement map (Fig. 4b) shows that the lava flow has a horizontal displacement maximum of 4 + 0.2 cm a 21 (1s), which is towards the centre of the lava flow, with a stronger eastward component. We explain this horizontal displacement as being due to the subsidence of the lava flow after its emplacement, and possibly also to gravity-driven compaction controlled by the local slope of the volcano oriented at N1308. The vertical displacement map (Fig. 4c) illustrates that the lava flow is subsiding at a rate of 13 + 0.3 cm a 21 (1s) over the time spanned by our CSK dataset (i.e. on a c. 3 years period of time, where the first interferogram used is 108 days after the eruption ended). We attribute this to three possible mechanisms. The first is the porous and viscoelastic relaxation and compaction of the surrounding medium in response to the gravity load, owing to the emplacement of the lava flow. This phenomenon is particularly clear near the border of the lava flow where there is narrowing relative to the actual area affected by the subsidence map. The second is thermal contraction of the lava flow coupled with thermal expansion of the surrounding medium owing to the cooling of the lava. Note that these two phenomena have competing effects. The third is porous and viscoelastic compaction of the lava flow itself owing to its own loading, plus settling of clasts in the 'a'a crusts. The effects of these mechanisms have been proposed for the post-emplacement behaviour of the lava flows of Okmok (Lu et al. 2005) and Etna (Murray 1988;Briole et al. 1997;Stevens et al. 2001).
In addition to the thickness and displacement maps, we also derived a map of the coefficient of determination, r 2 (Fig. 4d). Essentially, this map quantifies how well the actual data agree with the assumed inversion model (i.e. equations 4 & 5). Not surprisingly, the area covered by the lava flow exhibits a high value of r 2 , whereas the area beyond the lava flow is characterized by low r 2 values. These results support, a posteriori, our initial assumptions (linearity of displacements as a function of time, and linearity of Dw GST as an inverse function of altitude of ambiguity). In addition, the r 2 map shows intermediate values indicating partial agreement with the inversion model in a zone a few pixels wide surrounding the lava flow. In this boundary area, there is no correlation between the phase and the altitude of ambiguity, but, on the other hand, because the interferograms also recorded visco-elastic flexion of the substratum related to the loading of the nearby lava flow, then we can observe a correlation here between the phase and the duration of the interferogram.
The 1-map (Fig. 4e) shows the residual phase depicting erroneous and/or noise values that have not been fully accounted for by the assumptions of the inversion equation. The residuals are low across the whole flow field area, apart from in the centre, where they reach up to 3.5 radians. A possible explanation for the high residuals in this part of the flow can be seen in the plot of the vertical displacement rate against flow thickness (Fig. 5). Between 0 and 10 m the plot shows a clear positive correlation between lava flow subsidence rate and thickness. Such a positive correlation is to be expected because the subsidence measured on the lava flow is, at the first order, the sum of the three mechanisms mentioned initially. However, there is a significant deviation from this overall positive trend for pixels where the thickness exceeds 10 m. Note that these are also the same pixels where we observe strong residuals. This suggests that the subsidence may have been underestimated for the thickest part of the lava flow, most likely owing to the fact that strong displacements may have induced phase aliasing in some interferograms, making them difficult to unwrap properly.
As a consequence, we developed an iterative approach for improving the unwrapping of the interferograms. The basis of this approach was to reduce fringe aliasing by removing a simulated fringe pattern, assumed to account for the larger wavelength component of the signal, from the wrapped interferograms. This simulated fringe pattern was obtained from our first estimate of lava thickness and ground surface displacement, calculated using equation (4) and a low-pass filter. The residual interferograms were then re-unwrapped and the phases of the simulated fringe pattern added back into the result of the unwrapping. We ran this operation twice, hence significantly reducing the residual values in the 1-map and improving the r 2 map in several locations, particularly within the area of the flow emplacement (1 1st iter. max = 3.2 rad ⇒ 1 3rd iter. max = 3.0 rad, r 2 1st iter. mean = 0.79 ⇒ r 2 3rd iter. mean = 0.8). After improving the unwrapping step in this way, we obtained a final maximum thickness of 16.62 + 2.49 m (1s), a final average thickness of 5.85 m, and a mean error of 0.69 m (1s). The final volume estimate from the thickness map is 3.97 + 0.48 × 10 6 m 3 (1s), giving a DRE volume of 1.77 + 0.75 × 10 6 m 3 (1s) assuming a 55.46 + 13.53% vesicularity (1s). Although the iterative approach indisputably contributes to reducing the residual values in the 1-map, some unwrapping errors remain in the interferometric dataset, so the thickest part of the lava flow is still affected by high values of 1 (i.e. 1 3rd iter. max = 3.0 rad). As a consequence, a systematic bias on the thickness exists for these badly unwrapped areas.

Case 2: using TDM data
Since all the pairs were acquired after the October 2010 eruption, it was not possible to differentiate a post-eruption from a pre-eruption TDM interferogram, as used in the method of Poland (2014) to characterize emplaced lava flows at Kīlauea Volcano, Hawai'i, between 2011 and 2013. Instead, we applied here an approach similar to that for the CSK data processing, where the IGN 5 m LiDAR DEM is used as the pre-eruption topographic reference. However, because the images of bistatic pairs are acquired simultaneously in TDM, the atmospheric and ground surface displacement signals are negligible and the coherence of the radar data is generally close to 100%. The phase reflects only topographic residuals, so that equation (4) reduces to: The 0.5 ratio between the first term on the righthand side of equation (6) and the first term on the right-hand side of equation (4) accounts for the fact that the difference in the round-trip range of the radar wave between the master and slave in bistatic acquisitions is half the difference in the round-trip range in the case of monostatic acquisitions.  5. Plot of the vertical displacement rate as a function of the CSK-thickness map. Each point corresponds to a pixel within the lava flow. Assuming that the compaction of the lava flow at a given depth is a linear function of the load above, the total compaction is a quadratic function of the thickness; hence, we used a degree-2 polynomial fit. This curve fit presents a positive correlation between the subsidence rate and the thickness of the flow at 0 -10 m. However, beyond 10 m the subsidence rate deviates to a constant behaviour. The loss of positive correlation at the thickest part of the flow may be due to unwrapping problems encountered during the processing of the CSK interferograms, resulting in an underestimation of displacement and/or an over estimation of thickness, as these thickness values coincide with the pixels where strong residuals occur.
We derived the thickness on a pixel-by-pixel basis using the same least squares approach as used for the CSK data, but using in place of the variance-covariance matrix a weighted matrix with diagonal terms inversely proportional to the altitude of ambiguity of the interferograms and off-diagonal terms equal to zero. In this way, we gave a higher weight to the interferograms that are more sensitive to topography.
From the resulting thickness map (Fig. 6), we obtained a maximum thickness of 22.20 + 0.93 m (1s) at the location of the scoria cone built during the October 2010 eruption. This does not correspond to the location of the thickest part of the CSK thickness map. Actually the 1-map produced from the CSK dataset exhibits high residual values at the location of the October scoria cone, suggesting possible unwrapping errors there that probably result in an underestimation of the thickness at this location. The average thickness obtained from the TDM thickness map is 6.00 m and the mean error (mean of the thickness standard error calculated for each pixel) is 0.31 m (1s). The total bulk volume and the DRE volume were calculated to be 4.10 + 0.21 × 10 6 m 3 (1s) and 1.83 + 0.65 × 10 6 m 3 (1s), respectively. The thickness and the volume uncertainties were calculated using the same approach that we used in the CSK data processing.
Note that it is always possible to obtain thickness and volume estimates using one pair of bistatic images. However, here we have simultaneously used more interferograms to reduce the Gaussian noise present in the radar data. In Figure 7 we illustrate the relationship of the bulk volume uncertainty to the number of interferograms used. Even with two bistatic interferograms, one can already have a statistically satisfying estimate of the thickness and volume of the October 2010 lava flow when compared with volume estimations derived from ground measurements.

Case 3: TADR and volume estimations with MODIS data
Although MODIS operates in a total of 36 bands (i.e. 10 VIS (visible) bands, six NIR (near infrared) bands, four SWIR (short-wave infrared) bands, six MIR (mid-infrared) bands and 10 TIR (thermal infrared) bands), here we only used four bands: band 6 (1.628-1.652 mm); band 21 (3.929 -3.989 mm); band 22 (3.929-3.989 mm); and band 32 (11.770-12.270 mm). We adopted the normalized thermal index (NTI) of Wright et al. (2002) to automatically detect the hot spots, or anomalous radiances, owing to the presence of the active lava flow: If the NTI is greater than a fixed threshold, then the pixel is flagged as 'hot.' If a saturated radiance (i.e. pixel value of 21.000) is observed in band 22, the corresponding band 21 radiance value is utilized as a substitute. Band 21 has the same wavelength as band 22, and fortunately saturates at a higher temperature. Wright et al. (2001) noted that lava flow emplacement in a 1 km MODIS pixel resolution can increase radiance by up to c. 200% at 4 mm but by only 1% at 11 mm. However, by day, although the 4 mm signal will be contaminated by solar-reflected radiance, the 11 mm will not. This will increase the apparent brightness temperature difference between the two wavebands even further (Harris 2013). We used the MODIS band 6 reflection to convert to a 4 mm reflection-equivalent, and subtracted this from the band 22 and 21 radiances to correct the effect of solar reflection for MODIS images captured during the day. Although Wright et al. (2002) flagged pixels as hot spots if NTI values were greater than 20.6 (by night) or 20.8 (by day), these were conservative values set to allow complete, global detection. For local precision, case-specific values can be selected and set (Kervyn et al. 2008). For La Réunion, based on empirical observation, we set a threshold of 20.6. Thus, If NTI . −0.6, then pixel is ' hot' A phenomena peculiar to MODIS is the 'bow-tie effect'. This results from the image being acquired in sequential scans that are 10 pixels wide. Although at nadir there is no overlap between adjacent scans, this increases with scan angle up to an overlap of 100% at a scan angle of 558. This means that, at a scan angle of 458, pixel 1 of scan 1 would also appear in pixel 8 of scan 2. Over an active lava this means that, when the pixels are reprojected to their correct position, the hot spot is duplicated, appearing as a mirrored hot spot separated by two pixels from the real hot spot (Harris 2013). This results in double counting. Thus, to eliminate the spurious pixel clustering, we selected the cluster of anomalous pixels that contained the highest radiance value.
To estimate the hot area in a pixel we need to apply a single-band mixture model to estimate the portion of the pixel ( p) occupied by the hot source, whereby (Harris et al. 1997;Harris 2013) Here, M int (l) is the atmospherically and emissivity-corrected pixel-integrated radiance at wavelength l, and M hot (l) and M ambient (l) are the radiance emitted at wavelength l by the hot spot, in our case lava, and ambient background surrounding the hot spot, respectively. This required two assumptions: M hot (l) and M ambient (l). For M hot (l) we used the radiance equivalent of the minimum and maximum effective lava temperatures (2008C Fig. 7. The TDM-derived mean standard deviation of the bulk volume as a function of the number of interferograms used during the calculation of the thickness parameter. The plot indicates that even with two interferograms the mean accuracy is still satisfactory (c. 0.71 × 10 6 m 3 , or around 17% of the total estimated volume; this is c. 20 times less than when the area is estimated with a GPS and the thickness is inferred from spot measurements). Notice too that the range of the values is quite narrow. This further supports the reliability of the volume measurements derived using TDM data. and 8008C) for poorly crusted lava used by Harris & Neri (2002). This gives two end-member solutions within which the reality is assumed to lie (Harris 2013). For the background, Wright & Flynn (2003) suggested using the surface temperature from the band 32 pixel for the cold background value, under the assumption there was no contribution from the hot source at this wavelength. However, saturated pixels in band 22 can be powerful enough to raise the radiance value above ambient in band 32. To avoid this problem, for those pixels that were saturated in band 22, a four-pixel-neighbourhood averaging was performed on band 32 pixels. The band 22 or 21 radiance was then used for M int (l) to solve equation (8). The pixel portion occupied by the hot source was then multiplied by pixel area to obtain the hot area in each pixel. These were then summed for all hot pixels to obtain two values, one for each of the hot temperature conditions, for the total active lava area (A 200 and A 800 ). Now, assuming a steady state, we can potentially use equation (8) to convert the areas to equivalent TADRs leading to the derivation of two new relations for the Piton de la Fournaise, these being: TADR min and TADR max being the minimum and maximum bounds on our estimated TADR range, and 3.72 × 10 26 and 5.56 × 10 25 m s 21 being the linear relation between TADR and area for the 200 and 8008C temperature assumptions, respectively.
In Figure 8, TADR is presented as a function of Julian day. TADR starts at a very low value, which may be a result of a steady state not yet having been reached (Garel et al. 2012), but then increases rapidly to peak towards the third day. This waxing behaviour has, in the past, been explained by a pressurized magma body forcing itself to the surface by propagating a crack upwards (Wadge 1981;Rowland et al. 2003). The intrusion of a dyke during the October 2010 eruption is consistent with observations made by the OVPF during the eruption. Following the chronology of events based on the OVPF internal reports there was an increase in the depth and frequency of tremors before the eruption, with dyke propagation to feed a pressurized release of lava being the logical result. However, by 17 October volcanic tremors had decreased by 30%, and only one cone remained active; this implies waning activity. We thus expect that TADR should have dropped by this time (i.e. Julian day ¼ 290), which is indeed the case in the TADR plot. Following Harris et al. (1997), the single point peak on day 299 is likely to be an artefact, and thus should not be trusted (Harris 2013).
The total erupted volume can be calculated by integrating the two TADR curves, using a simple trapezium rule (Harris 2013). This gives us a MODIS-derived range for the total DRE volume of between 2.44 and 4.40 × 10 6 m 3 . This overlaps with the InSAR-derived values.

Conclusion and further discussions
We have presented three different techniques that use InSAR and thermal remote sensing data to derive lava flow area, volume and subsidence information during and following the October 2010 effusive eruption of Piton de la Fournaise.

During the eruption
The syn-eruptive TSX SAR data, when coupled with an updated and high spatial resolution DEM, can produce interferograms that allow deformation induced by an eruption to be detected. Modelling the observed displacements from the interferograms can be very useful to provide constraints on the source geometry and, if combined with the mass budget of the October 2010 eruption, it can then take into account both the effusive and intrusive part of the emitted lava. Here, we measured a maximum EW displacement of 40 cm and an uplift of up to 30 cm. The presence of both opening and inflation suggests that the magma body forced itself towards the Earth's surface as a dyke which then fissure-fed the lava flow. This can be further supported by the waxing and waning signature found in the MODIS-derived TADR plot. The integrated ground surface uplift yields a minimum volume of magma injection of 0.8 × 10 6 m 3 .
By exploiting the coherence of two syn-eruptive interferograms, we were also able to map the emplacement of the lava flows regardless of weather conditions and time of acquisition (i.e. day or night). Provided that one can obtain SAR data in nearreal time, such an exploitation of interferometric coherence can be utilized to determine flow area, shape and extent during an eruption. Moreover, as compared with airborne surveys, this method is cheaper and overpasses are guaranteed at known and regular times.

After the eruption
We used two types of InSAR datasets to characterize, a posteriori, the October 2010 lava flows: a large database of CSK monostatic interferograms and a limited number of TDM bistatic interferograms. Each of the data sets has its own strengths and weaknesses in terms of information that it can provide about the lava flow. We compared the volumes of the lava flow obtained from both InSAR datasets with an approach where the volume was derived from a thermal remote sensing method based on volume flux estimates. We also compared our results with several sets of published data. All results are of the same order of magnitude, but the volume determined using the InSAR data has an uncertainty level at least one order of magnitude lower than the other approaches.
From the CSK dataset, we were able to measure displacements associated with subsidence of the October 2010 lava flow 8 months after its emplacement. After attempting to correct phase unwrapping problems, we found that the subsidence of the lava flow occurred at a maximum rate of 13 + 0.3 cm a 21 . There was also a small horizontal movement which resulted from a combination of displacement towards the centre of the lava flow and a slight eastward sliding, at a rate of 4 + 0.2 cm a 21 , owing to the local slope (i.e. N1308). We attribute the subsidence to the following causes: † viscoelastic relaxation of the surrounding medium owing to gravity loading and compaction; † thermal contraction of the lava flow on cooling; † porous and viscoelastic compaction of the lava flow itself owing to self-loading.
The thickest part of the flow was found near the vents, as depicted by the thickness maps obtained from the two InSAR databases. The final mean thicknesses are about 6 m for the CSK data and the TDM data. Thickness (flow front thickness) can be used along with the underlying slope to derive the yield strength and viscosity of the lava flow. It is also very helpful when it comes to lava flow simulation models. Another important application for the thickness map is the regular updating of DEMs, which is necessary for lava flow simulations. This is especially useful for volcanoes that erupt frequently as their morphology is known to change rapidly, so that DEMs need to be changed regularly, even daily, during ongoing effusion. Additionally, using the thickness data, we estimated the volume of the erupted lava flow with accuracy. In Table 1 we summarize the different volume estimates that have been derived from this and other studies, indicating a DRE volume of between 1.02 and 4.40 × 10 6 m 3 erupted over a duration of 16.4 days, to give a MOR of 0.7-3.1 m 3 s 21 . The volumes calculated from CSK and TDM datasets fit well within the range of the thermal remote sensing technique. In terms of volume uncertainty, the TDM and CSK datasets are also very close, but it is worth mentioning that the CSK approach requires many interferograms that cover a large range of altitude of ambiguity and long duration to allow simultaneous characterization of the thickness and displacement rates of the lava flow in a statistically satisfactory way. Obviously the production of such a large interferometric dataset requires several months, which is incompatible with near-real-time monitoring. On the other hand, the TDM technique requires a minimum of two interferograms to produce more accurate thickness map and a volume estimate in comparison to ground measurement-derived estimates.
The procedure that we have implemented in order to reduce the unwrapping errors in the CSK dataset significantly enhances results, as shown by the improved 1-map and r 2 -map towards the third iteration. This is another advantage of using bistatic TDM interferograms rather than the classical monostatic interferograms if the signal of interest is ground surface topography, for which the unwrapping step is achieved without any particular difficulties.