The Use of Massive Deformation Datasets for the Analysis of Spatial and Temporal Evolution of Mauna Loa Volcano (Hawai’i)

Abstract

In this work, we exploited large DInSAR and GPS datasets to create a 4D image of the magma transfer processes at Mauna Loa Volcano (Island of Hawai’i) from 2005 to 2015. The datasets consist of 23 continuous GPS time series and 307 SAR images acquired from ascending and descending orbits by ENVISAT (ENV) and COSMO-SkyMed (CSK) satellites. Our results highlight how the joint use of SAR data acquired from different orbits (thus with different look angles and wavelengths), together with deformation data from GPS networks and geological information can significantly improve the constraints on the geometry and location of the sources responsible for the observed deformation. The analysis of these datasets has been performed by using an innovative method that allows building a complex source configuration. The results suggest that the deformation pattern observed from 2005 to 2015 has been controlled by three deformation sources: the ascent of magma along a conduit, the opening of a dike and the slip along the basal decollement. This confirms that the intrusion of the magma within a tabular system (rift dikes) may trigger the sliding of the SE portion of the volcanic edifice along the basal decollement. This case study confirms that it is now possible to exploit large geodetic datasets to improve our knowledge of volcano dynamics. The same approach could also be easily applied in other geodynamical contexts such as geothermal reservoirs and regions with complex tectonics.

1. Introduction

Earth Sciences have experienced in the last decades a significant growth of in situ and remote sensing measurements useful for understanding the ongoing natural and anthropogenic phenomena. Major geodetic and seismic networks UNAVCO (http://www.unavco.org) and IRIS, (https://www.iris.edu) data archives, as well as large remote sensing data archives (Landsat satellite, http://landsat.gsfc.nasa.gov; ERS, ENVISAT and Sentinel satellites) represent only an example of the large volume of data presently available to the Earth Sciences community.

Recently, several initiatives devoted to ease the dissemination, processing and exploitation of remote sensing data have become operative (Supersites, http://supersites.earthobservations.org; Geohazards Exploitation Platform https://geohazards-tep.eo.esa.int; Google Earth Engine, https://earthengine.google.com; Amazon AWS, https://aws.amazon.com/it/public-data-sets/landsat/). The availability of Earth Sciences large dataset has produced an imbalance between the data production and their effective exploitation in research activities. This is particularly crucial in the Earth Sciences, where the partial use of the available information for a specific area implies a partial understanding of the observed phenomena since the available data monitor different effects of the same phenomena. For instance, during a volcanic eruption, it is possible to collect several kinds of data [1] (e.g., gas emissions, ground displacements, seismic tremor, etc.), all of them describing a particular aspect of the same phenomenon.

The joint analysis of heterogeneous data is the desirable solution for understanding geological phenomena in the broad sense. Nevertheless, the joint exploitation of heterogeneous data to effectively extract value added information is still an open task [2,3,4,5,6,7] and, therefore, the development of algorithms and procedures that synthetize this large variety of data is needed.

The goal of this work was the analysis and inverse modeling of Volcano surface deformation through the joint exploitation of geological information, in situ dataset (GPS) and large Differential Synthetic Aperture Radar Interferometry (DInSAR) measurements.

In this regard, one of the most interesting areas for the Earth Science community is the Mauna Loa Shield Volcano, located on Island of Hawai’i (USA) (Figure 1). It is a very active and deforming vlcano for which a large variety of studies and data is available. Indeed, the structural geological and geophysical information of Mauna Loa Volcano, available in literature, to constraint geometry and position of the source configurations can be retrieved, among others, from: [8,9,10,11,12,13,14,15].

Moreover, Mauna Loa is monitored by a permanent GPS network (http://hvo.wr.usgs.gov/maunaloa/current/gps.html) and, since the early 2000s, a large SAR archive has been collected from several satellites, such as the C-band ENVISAT of the European Space Agency (ESA), RADARSAT-1/2 of Canadian Space Agency (CSA), the L-band ALOS-1/2 of Japan, until the X-band German TerraSAR-X and Italian COSMO-SkyMed (CSK) constellations. More recently, the C-band Sentinel-1 sensors of the European Copernicus program have also started to systematically acquire over the area with a repeat pass of 12 days. In particular, we exploited the ENVISAT and CSK satellite data, which represent one of the most comprehensive and available SAR archives over the Volcano. To efficiently process the SAR datasets, we used the Differential SAR Interferometry (DInSAR) method of Parallel Small BAseline Subset (P-SBAS) [16], which, benefiting from High Computing Performance facilities, allowed us to generate mean deformation velocity maps and time series by the automatic processing of large volume of SAR data, thus making the generation of DInSAR result extremely fast [17].

Finally, to model the time-varying source volumes at Mauna Loa Volcano, we performed linear inverse modeling by exploiting the large amount of GPS and DInSAR measurements. We proposed a complex model having a specific spatial configuration and constrained by the available geological information. In particular, within the volcanic system, we simultaneously modeled three spatially and temporally variable elements: the tabular dike system, the central pipe conduit and the basal decollement. This represents an innovative method compared to previous works such as [18,19,20], where simple active sources were considered to simulate the observed surface deformation.

2. Geological Background

Mauna Loa is the largest active Volcano on Earth, having erupted 33 times since its first well-documented historical eruption in 1843.

Its morphology is characterized by the Moku’aweoweo summit caldera and by two opposing elongated ridges: the Northeast Rift Zone (NERZ) and the Southwest Rift Zone (SWRZ) (see Figure 1, south and east branches coincide with SWRZ and NERZ, respectively) [21]. The rift zones have widened because of repeated dike intrusions and gravity driven deformation [22]. This “spreading” has been accommodated by displacements along a low-angle thrust fault at the base of the edifice, the so-called “decollement”, which coincides with the interface between the oceanic crust and the volcanic edifice [23,24]. Dike intrusions into the rift zones are thought to compress the adjacent flanks of Mauna Loa [11,25], causing intermittent flank movement and large earthquakes at the decollement [26]. Large earthquakes, in turn, may open magma pathways and facilitate dike intrusions and eruptions [12,27].

According to the geophysical model of Mauna Loa’s magma reservoir proposed by [15], most of the magma is generated in the mantle, possibly at the depth of the subcrustal seismicity, pond in the deep and intermediate section of the rift zone with only a small percentage collecting into the shallow magma chamber. Intrusion of the magma changes the stress field within the volcanic edifice encouraging dike propagation into the shallow SWRZ and faulting along subhorizontal faults in most parts of the volcanic edifice, including the decollement. Mauna Loa exhibits a cyclic behavior where the deformation due to earthquakes and intrusions helps new intrusions in the rift zone and fault slip under the flanks. [28] showed how the dike propagation at Mauna Loa occurs both vertically and laterally and it is controlled by numerous factors: magma pressure, the preexisting stress field, magma viscosity and host rock properties [29,30,31]. When dike driving pressure (defined as the difference between the magma pressure and the least compressive stress in the crust) is sufficiently high, the dike will reach the surface and erupt; however, geologic evidence suggests that most dikes arrest before they reach the surface, probably because mechanical and stress barriers inhibit propagation [32,33,34,35,36,37]. Earlier models based on Mauna Loa’s two recent eruptions made a similar prediction, placing the chamber at 3 km deep at the approximately same geographic position.

The Global Positioning System (GPS) measurements obtained by permanent stations of the Hawai’i monitoring network, processed by the Nevada Geodetic Laboratory of the Hawai’i monitoring network (http://geodesy.unr.edu/NGLStationPages/GlobalStationList) (used in this paper) show southeastward movement of the southeast flank of about 4 cm/year, comparable to displacements measured on adjacent Kilauea Volcano’s south flank over the same interval. The upper west flank appears to be stable, producing a strong asymmetry of motion on the summit. Gradients of motion on the southeast flank result in about one micro-strain/year of compression and shear across the Kaoiki seismic zone, an area of persistent seismicity that has produced large historic earthquakes. The flank motions observed could be caused by the combined effects of slip along the basal Kaoiki decollement and inflation of a deep source [21,38].

3. Data Analysis

We investigated the source of the ground deformation observed at Mauna Loa by exploiting both GPS (above mentioned) and multi-platform/multi-orbit SAR sensors. In particular, we analyzed the deformation time series of 23 GPS (Figure 2) permanent stations of the Hawai’i monitoring network, operating since 2000. All the GPS measurements were referenced to the MAUI GPS station, located about 155 km away from Mauna Loa Volcano.

In addition, we considered seven SAR data stacks, acquired by the C-Band (wavelength of ~5.6 cm) ENVISAT satellite in 2003–2010 (Figure 3 and

Table 1. ENVISAT and COSMO-SkyMed SAR datasets processed with the G-POD environment and stand-alone processing, respectively.

Sensor	Mode	Track	Image Swath	Number of SLC	Number of Interferograms	Time Interval *
ENVISAT	Descending	472	I1	36	106	31122005–11092010
ENVISAT	Ascending	136	I3	12	34	27012005–23092010
ENVISAT	Descending	343	I7	18	50	17112005–07102010
ENVISAT	Descending	200	I2	42	123	08122003–27092010
ENVISAT	Ascending	408	I4	36	105	27122005–12102010
ENVISAT	Descending	114	I6	12	48	25012005–17082010
ENVISAT	Ascending	365	I2	27	76	08022003–04092010
CSK	Ascending	------	HI_06	51	154	13032012–07092014
CSK	Descending	------	HI_04	67	200	03032012–05092014

* dates are in the format: ddmmyyyy.

) along ascending and descending tracks, with different look angles (from I1 to I7, see Table 1) for a total of 189 SAR images. Finally, we employed two data stacks collected by the X-Band (wavelength of ~3.1 cm) COSMO-SkyMed constellation during the 2012–2015 time span, for a total of 118 images (Figure 3 and Table 1).

The SAR datasets were processed through the P-SBAS technique [16], which allowed us to retrieve the Line of Sight (LOS) projection of the surface deformation and analyze its temporal evolution by generating displacement time series.

The ENVISAT datasets were processed by using the on-line P-SBAS web tool available within the ESA’s Grid Processing On Demand (G-POD) environment, in the framework of the ESA’s Geohazards Exploitation Platform (GEP) [17,39], while the CSK datasets were processed employing the same P-SBAS processing chain on in-house computing facilities, because of the unavailability of CSK data on the G-POD platform.

It is worth noting that the G-POD approach not only allows quick and remote processing of a large quantity of ENVISAT data in an automated and unsupervised way, but also offers the possibility to reproduce the P-SBAS results employed in this work, since the platform is openly accessible to the scientific community [17].

The ENVISAT datasets consist of SAR stacks ranging from 12 to 42 acquisitions (Table 1), with interferograms sequences spanning the range from 33 to 123. The satellite look angles among the exploited ENVISAT tracks varied from 15° (I1) to 38° (I7), allowing us to image the ongoing displacement from different viewing angles and, therefore, to better characterize the contributions of the different ground deformation components (Figure 4). The P-SBAS results from the ENVISAT data were spatially averaged (multilooked) to obtain a pixel size on the ground of about 80 m by 80 m.

The CSK datasets are composed by 51 and 67 Stripmap images which have been acquired by ascending (look angle of ≈ 32°) and descending (look angle of 29°) orbits, respectively. We computed 154 interferograms from the ascending orbits and 200 from the descending orbits. These were subsequently inverted by applying the P-SBAS approach to produce average velocity maps and deformation time series for both orbits. A multilooking factor of 10 in both azimuth and range directions was applied, thus reducing the final resolution to about 30 m by 30 m.

For the purposes of the inverse modeling, all the DInSAR-SBAS dataset have been spatially dowsampled using a regular grid with a 1000 m spacing.

Figure 4 shows the ENVISAT (Figure 4A–G) and COSMO-SkyMed (Figure 4H–I) mean velocity maps acquired from both orbits. The C-band sensors reveal a large range of LOS mean velocities. Such a variation is mainly due to the usage of many datasets with various look angles, which implies a different projection of the deformation over the sensor LOS. In other words, the swath with a narrower look angle has a larger sensitivity to the vertical component of the deformation with respect to the larger look angle swath that, instead, is more sensitive to the corresponding horizontal component [39,40,41,42,43,44]. The retrieved velocity measurements move from a maximum of about 4 cm/year (corresponding to a sensor-target range decrease) at the Volcano summit to a minimum of 1 cm/year. Similar velocity values were also observed along both ENVISAT descending (Figure 4D–G) and COSMO-SkyMed ascending and descending orbits (Figure 4H–I). It is worth noting that the reference points for all the mean velocity maps, located in a stable area of the Big Island, are different for each dataset, but are all very close the MLPR GPS station.

Because of the lack of both GPS an DInSAR measurements during 2003–2005, we considered in our modeling only the data from 2005 onward.

4. Inverse Modeling

4.1. Methods

The method was based on the joint inversion of GPS and DInSAR measurements to image a complex and time-varying deformation source.

We considered each observation d of the dataset, being a relative measurement of an absolute deformation d^A [45]

$$d\left(x,y,t\right)=d^A\left(x,y,t\right)-d^A\left(x,y,t_0\right)-\left[d^A\left(x_R,y_R,t\right)-d^A\left(x_R,y_R,t_0\right)\right]$$

(1)

where d is the data at time t and at the point (x,y) on the surface. Conversely, d^A represents the absolute ground deformation in a given reference system. In the second term of Equation (1), we subtract the ground deformation measured at a reference time t₀; specifically, in our analysis, the reference time differs for each dataset: either 1 January 2005, if the dataset were already available before this date, or the time of the first time series date (Figure 3D). The terms within the square parenthesis express the how time-varying deformation, starting by reference time t₀.

Now, we express generic data as a simple linear model:

$$d_k\left(x,y,t\right)={\displaystyle \sum }_p{a}^p\left(t\right)G_k^p\left(x,y\right)$$

(2)

where d_k is the data belonging to the k-th dataset. The sum is intended over the p elementary sources responsible for the time-varying ground deformation. The term a^p is a time-varying amplitude coefficient for each source (e.g., the volume change for a point source or the slip for an Okada dislocation model), while G_k^p is the Green’s function for the spatial ground deformation pattern of the p-th source for the k-th dataset.

We express both GPS and DInSAR Green’s functions G_k^p in a similar way as:

$$G_k^p\left(x,y\right)=G^p\left(x,y\right)v_k$$

(3)

where the vector G_p is the ground deformation related to the p-th source and v_k is the versor associated with the k-th dataset. For DInSAR measurements, the versor coincides with the LOS direction. For the GPS data, we considered each ground motion component separately, hence, in this case, the three v_k versors are simply the east, north and up direction vectors.

In the proposed modeling approach, the geometry of the sources was fixed, which means we assumed their position and their overall size does not change with time.

The unknowns of Equations (1) and (2) are the time-varying amplitude coefficient for each source a^p. Each coefficient was discretized in n time intervals within which they were interpolated linearly. We discretized the interval from 2005 to 2012 with a time-step of one year, and the interval from 2012 to 2015 with a time-step of six months, because of the higher frequency of the data (Figure 2B and Figure 3D). Substituting Equation (3) into Equation (2) and then into Equation (1), and considering the aforementioned discretization, we obtained a big system of linear equations. Indeed, the unknowns of the inverse problem are the P × N discrete coefficients a_n^p, where P represents the total number of elementary sources and N is the number of considered time intervals. In our problem, the corresponding linear system has 1652 unknowns (p = 118 and n = 14) and about 3.6 × 10⁵ data. The solution of this overdetermined system has been obtained by Truncated Singular Value Decomposition (TSVD) [41]. The number of significant eigenvalues has been selected by using the L-curve approach [41,42].

For what concern the data error, the measurement precision for the SAR is about 1–2 mm/year for deformation velocities and 5–10 mm for the displacement time series [46]; the error for the GPS measurements is about 2.5 mm and 10 mm for horizontal and vertical component, respectively. The data covariance matrix used in the formulation of our inverse problem is a diagonal matrix whose components are the reciprocal of the data variances.

4.2. Source Geometry

To determine the best-fit geometric configuration that fits the observed ground deformation at Mauna Loa Volcano, we employed the Akaike Information Criterion (AIC) [47] statistical analysis. This criterion quantifies the balance between goodness-of-fit and complexity (i.e., number of parameters) for a given model. In the context of statistical model selection, the model possessing the lowest AIC is favored according to Occam’s razor principle. To reduce as much as possible the number of free parameters, the sources geometric configuration (i.e., geometry and position of the sources) was retrieved starting from the “a-priori” available geological and geophysical information; in particular, in our model we considered:

• the volcanic conduit, modeled by a closed pipe source [10,48,49];
• the basal decollement modeled by a dip slip dislocation [9,13,50];
• the rift zone and the associated dikes system, modeled by an opening tensile dislocation [8,11,12,51]; and
• the shallow magma chamber, modeled as a point source [52].

The presence of this magma chamber was suggested by several authors based on the inversion of geodetic data, seismological and petrological data [15,53,54,55]. We determined the position of the magma chamber by using a non-linear optimization procedure based on the [56] Simplex algorithm by exploiting the considered datasets. The performed optimization procedure located this source at a depth of about 5 km beneath the summit caldera.

In this scenario, we checked the statistical reliability of the sources geometric configurations by testing various sources combinations using the AIC (

Table 2. Akaike Information Criterion for five source combinations. Columns 2–5 indicate which sources were considered. Column 6 indicates the number of eigenvalues selected through the L-curve method. Column 7 is the residual-sum-of-squares (RSS) for each model and the last one indicates the relative AIC value. The best-fit sources geometric configuration (Model 1) is evidenced in bold.

Model	Mogi	Pipe	Dike	Decoll.	# eig.	RSS	AIC
1	N	Y	Y	Y	127	22,99275438	0
2	Y	N	Y	Y	124	26,27011658	17355.3
3	Y	Y	N	Y	133	24,57768648	8697
4	Y	Y	Y	N	132	23,11053503	675,7
5	Y	Y	Y	Y	131	22,99243673	6,2

).

For each sources combination, we performed an inversion using a TSVD technique and selecting the optimal number of eigenvalues through the L-curve method. It is evident in Table 2 that the combination having the lowest AIC involves the configuration with a volcanic pipe, the opening of a dike and slip along the basal decollement (model 1). The presence of a shallow magma chamber is not supported by the considered datasets, since its inclusion does not improve statistically the solution: actually, it increases slightly the AIC value (see difference between Models 1 and 5 in Table 2). Thus we excluded the Mogi source [57] from our analysis since its presence is not rigorously supported by the inversion of our datasets. On the other hand, the exclusion of one of the other sources (i.e., pipe, dike or decollement; Models 2–4 in Table 2) makes the residuals and the AIC significantly high, highlighting that the role of these sources is essential.

In the top part of Figure 5, we show the best-fit model configuration, displaying the geometry of the key volcano–tectonic elements used in the inverse modeling and the relevant spatial scales.

The feeding conduit located beneath the summit caldera was modeled by four consecutive sections, each of them modeled through the closed pipe model of [48].

The heights and depths of each section of the pipe conduit coincide with those used for the dike system.

The geometry and location of the basal decollement (II) beneath Mauna Loa were based on [9]; we modeled the decollement as a set of patches arranged on a regular grid over a surface dipping of about 2° toward NW and having a depth of about 10 km beneath the summit caldera. Each patch is 10 km × 10 km and the total number of patches is 34 (top panel of Figure 5). The slip on each patch occurs along the SE direction, nearly orthogonal to the rift system. The ground deformation was modeled using the Okada dip slip fault model.

The dike system (III) geometry is approximated by 80 patches (or tiles) arranged on a grid 20 km × 4 km. The length of each patch is about 6 km. The vertical width of each patch is variable, being about 0.5 km at the top and about 4 km at the bottom. The patches extend from 0.5 to 10 km depth beneath the surface (top panel of Figure 5). This is because of the decrease of resolution with depth which is known to affect geophysical imaging problems [58]. The trace of the dike extends 42 km NE and 65 km S of the summit caldera, for a total length of about 107 km, following the morphological elements marking the Mauna Loa rift [59,60,61] (Figure 1). A Poisson coefficient μ = 0.25 was assumed.

4.3. Results

Similar to what several authors [18,19,62,63,64,65,66,67] have already done, we exploited the methodology addressed to the source volume change detection through an integration of different datasets.

Volume increases, responsible for the observed ground deformation pattern, are indicated for both the central pipe and the dike system, using a color scale, in Figure 5. The slip along the basal decollement is represented by the arrows, proportional to the slip magnitude.

We estimated the “a posteriori” covariance matrix over the model 1 parameters (Table 2) [68] to estimate the average uncertainties over the volume change and slip: 0.11 × 10⁶ m³ for the volume change of each tile (patch) of the dike system, 0.15 × 10⁶ m³ for the volume change of each section of the central pipe conduit and 0.2 cm for the slip along each patch of the basal decollement.

In Figure 6, we compare the model results with the DInSAR LOS time series of the pixels evaluated on each dataset and localized the closest MLSP GPS station. We selected this station because it is close to the area of maximum ground deformation.

The comparison between data and retrieved time-varying model for GPS is shown in Figure 7 and Figure 8. In Figure 7, we compare the temporal variations on the vertical component of the GPS displacements (in black) with the final model (red). In Figure 8, we plot on a map the time-varying horizontal components, relative to each GPS station (blue dots) compared with the final model (red).

The results show a good fit between the model and vertical and horizontal displacements: the DInSAR LOS time series (Figure 6), the vertical (Figure 7) and the horizontal GPS components (Figure 8) shows a maximum difference between model and observation of only few centimeters. The overall reduction of the variance achieved is more than 80%.

5. Discussion

Figure 9 shows the volume changes of the pipe (i.e., conduit) and the tabular area (i.e., dike system), compared with the seismic moment released by the slip processes along the decollement. From 2005 to 2006, the pipe undergoes a rapid inflation, which is followed by a delayed inflation of the tabular system and by a slip along the decollement. After the initial rapid inflation, the pipe continues to slowly inflate until 2010, even if with a rate being an order of magnitude lower than the previous interval. On the other hand, the inflation of dike system continues until 2010 with a similar rate. From 2007 to 2010, the slip along the decollement stops.

Since 2010, the pipe and the tabular system has undergone a slow deflation, probably related to cooling processes, while decollement has continued to slip, although at a reduced rate with respect to 2005–2007.

The volume variation of the dike system is only one third of the total inflation along of the pipe. The seismic moment release of the decollement until 2015 can be compared to the seismic moment of an earthquake with magnitude 4.7; the retrieved maximum slip along the decollement is of the order of few centimeters.

The uncertainties over the integrated of the volume changes reported in Figure 9 are 0.57 × 10⁶ m³ for the dike system, 0.23 × 10⁶ m³ for the pipe conduit and 1.9 × 10¹⁴ Nm for the seismic moment.

Regarding the pipe conduit source, Figure 10 shows the detailed volumetric variations as function of time and depth. The result highlights that most of the conduit inflation occurs at depths between 1 and 6 km with the larger volume change at depths above 3 km. However, there is no clear evidence of volumetric vertical variation along the conduit that seems to inflate as a single body. Several authors modeled this sector of the volcanic conduit as a nearly spherical source [15].

The closed pipe model [48] gives rise to a wider deformation pattern with respect to a spherical source at the same depth. This could explain why our results found that most of the magma seems to be stored at relatively shallow depths beneath the caldera (less than 3 km), while a spherical magma chamber has been located a depth of about 5 km.

The opening of the dike system is a consequence of magma transfer from the central conduit (Figure 11). Most of the opening occurs within 5 km from the central pipe source at the depths lower than 3 km. The retrieved maximum rift opening is about 1 m in the shallowest portion of the rift. The opening propagates up to 10 km in the E branch of the tabular dike system and up to 25 km toward S (Figure 11). The relevant inflation occurred between 2005 and 2006 and additional opening episode happened from 2006 to 2007 as well.

For this scenario, in Figure 12, we show the volume change (integrated along the depth) for the whole tabular system as a function of time. In addition, in this case, there is no clear evidence of either vertical or horizontal magma propagation along the tabular system. This result could be controlled by the time scale used to discretize the time intervals (i.e., one year for 2005–2012 and six months for 2012–2015): it could be too large to image rapid magma processes. For example, during the 2014–2015 eruption at Bardarbunga (Iceland), a dike propagated for about 45 km in 14 days [4]. This timescale is compatible with what can be inferred from our current results on Mauna Loa Volcano.

6. Conclusions

Our work illustrates the potential of the use of large ground deformation datasets and the application of joint inversion techniques to provide a detailed picture of the intrusion processes occurring within volcanic systems. Since geodetic inverse problems have intrinsic non-unique solutions [69], the use of multiple large datasets can help to better constrain the inversion results.

Previous studies on ground deformation at Mauna Loa focused on the nature of the magmatic feeding system of the volcano or on the interaction between the magmatic system and the basal decollement [11,15,21,28,38]. In this work, we proposed, for the first time, a comprehensive advanced model which considers the complexity of the magmatic system and the kinematics of magma vertical ascent lateral propagation and, at the same time, highlights the role of the basal decollement in the overall deformation pattern. In this perspective, we used datasets consisting of both GPS and DInSAR ground deformation time series, for a total of about 3.6 × 10⁵ data to image the ascent of the magma along: (i) the main feeding conduit of Mauna Loa; (ii) its propagation along the rift system; and (iii) the sliding of the SE flank of the Volcano along the basal decollement.

In particular, we have shown how the modeling of the whole dataset, based on statistical inference, seems to exclude the presence of a large spheroidal magma chamber beneath Mauna Loa, favoring an alternative model. Our model indicates that magma rises along a central volcanic conduit and then propagates both southward and eastward along Northeast Rift Zone and the Southwest Rift Zone: up to 25 km toward south and about 10 km within the E branch from 2005 to 2007. We believe that future efforts should be devoted in fully exploiting the most recent DInSAR datasets with shorter revisit time (i.e., Sentinel-1) with the aim of studying short time scales (weeks) over long time interval (years). This would increase the temporal resolution of the inversion and capture the short-term phenomena as the magma propagates within conduits and dikes. Furthermore, the use of long temporal series would allow considering visco-elastic rheologies, which could provide a more accurate quantification of the volumes involved in the ground deformation sources.