Robust Bayesian Joint Inversion of Gravimetric and Muographic Data for the Density Imaging of the Puy de Dôme Volcano (France)

1 Introduction

The density structure of volcanoes is classically inferred from the inversion of gravimetric data (Camacho et al., 1997; Cella et al., 2007; Linde et al., 2014). Gravimetry provides measurements of the gravity field throughout the study area, corresponding to the integrated effect of the whole Earth, but also sensitive to the local density variations. The inversion of gravimetric data is well-known to be non-unique and ill-posed, requiring a priori geological information or a combination with other geophysical data, such as seismic travel times (Onizawa et al., 2002; Coutant et al., 2012), to constrain the models. Muography is a recent imaging method which emerged from the field of particle physics (Alvarez et al., 1970; Nagamine et al., 1995). It provides images of integrated densities of an edifice using atmospheric muons generated by the interaction of cosmic rays with the atmosphere. The muons are charged elementary particles, like electrons but about 200 times heavier, that interact with matter in a stochastic way depending on their energies and on the density and composition of the medium (Groom et al., 2001; Nagamine, 2003). At high energy, traveling along straight paths, they can penetrate up to kilometers of rocks. Their trajectories are reconstructed thanks to multi-layer muon detectors. As muons come from the atmosphere, the muon detectors can only image structures higher in elevation. This makes muography particularly suitable to image lava domes like the Puy de Dôme (Chaîne des Puys, France). Application of muography to the imaging of volcanoes has been developed in the last two decades (Tanaka et al., 2001; Lesparre et al., 2012; Kusagaya and Tanaka, 2015; Oláh et al., 2018). The rate of muons crossing the target along a given direction depends, to the first order, on the rock density integrated along this line of sight. The 3D density reconstruction by muon tomography is limited by the number of view points and the acquisition duration. For example, imaging a 1,000 m thick and 1,800 kg/m³ dense structure with a precision of 5% on the density, an angular resolution of 3° by 3° and a 1 m² detector requires an exposure of about 100 days of effective acquisition time (Cârloganu et al., 2016; Niess et al., 2016).

The combination of muography and gravimetry in a joint inversion scheme better constrains three-dimensional density models (Davis et al., 2011; Davis and Oldenburg, 2012; Okubo and Tanaka, 2012; Nishiyama et al., 2014b; Rosas-Carbajal et al., 2017; Cosburn et al., 2019). However, the joint inversion problem is still ill-posed and a regularization is needed, tuned by some parameters. These regularization parameters are determined either based on a checkerboard test (Nishiyama et al., 2014b), based on a classical L-curve (Rosas-Carbajal et al., 2017), or arbitrarily (Nishiyama et al., 2017a). However, so far inversions have shown artifacts related to the muography acquisition geometry and to the limited number of muon detectors. Based on synthetic data, Barnoud et al. (2019) designed a Bayesian inversion scheme where two regularization parameters, an a priori density standard deviation and a correlation length, can be determined in a robust way using a Cross-Validation Sum of Squares criterion, such as the Leave One Out (LOO). Using this approach, the resulting 3D density models are free of artifacts linked to the acquisition geometry, even with a limited number of muographic view points. Besides, comparing the inversions of synthetic data with one and three muographic view points, Barnoud et al. (2019) show that the resolution of the resulting density model is improved when using multiple points of view. Another difficulty in such a joint inversion is that the density estimated by muography is often lower than the density estimated by gravimetry due to the detection of 1) non-ballistic low-energy muons scattered by the volcanic edifice (Nishiyama et al., 2014a, 2016; Gómez et al., 2017; Rosas-Carbajal et al., 2017), 2) muons coming from the backward direction (Jourde et al., 2013) and 3) other charged particles (Oláh and Varga, 2017; Saracino et al., 2017). Tests on synthetic data (Barnoud et al., 2019) show that this offset leads to artifacts in the inversion results so that it should be corrected. To overcome this, the determination of a constant relative density offset can be taken into account in the inversion process (Rosas-Carbajal et al., 2017). Using synthetic data, Lelièvre et al. (2019) explore and compare various methods to automatically determine the offset and recommend a least-squares approach. In this paper, we combine the most robust and efficient techniques identified by Barnoud et al. (2019) and Lelièvre et al. (2019). We use a Bayesian formalism with a least-squares automatic computation of a potential constant density offset. The regularization parameters, namely the a priori density standard deviation and the isotropic correlation length, are determined using an LOO criterion.

We apply our method to gravimetric and muographic data acquired on the Puy de Dôme volcano to recover the 3D density structure of the edifice. The Puy de Dôme volcano is the highest volcano (1,465 m a.s.l.) of the Chaîne des Puys located in the French Massif Central (Figure 1, top left). The Chaîne des Puys is a complex field of about 80 monogenetic volcanic edifices aligned following a north-south direction. These volcanoes are built on a Variscan granitic basement whose top is situated around 800–1,000 m a.s.l. (Portal et al., 2016). The word “puy” refers to an isolated hill in the Auvergne area of the Massif Central. The Puy de Dôme is an 11,000-year-old monogenetic, composite trachytic dome (Boivin et al., 2017). The dome itself is about 400 m high and is about 1,800 m wide at its base (Figure 1). The Puy de Dôme volcano is a good natural laboratory to assess the methods for joint gravimetric and muographic inversions as high quality gravimetric and muographic datasets are available. Moreover, it is extinct and its structure does not change, ensuring a stable image during the six months of muography campaign required to get a sufficient resolution. With up to 1,600 m of rock thickness along the muon lines of sight, it is, to our knowledge, the largest lava dome presently imaged combining muography and gravimetry.

FIGURE 1. Simplified geological map of the Puy de Dôme and neighboring cones, adapted from Boivin et al. (2017). Top left: Location of the Chaîne des Puys in France. Blue square: Study area of the present paper (Figures 2, 7).

2 Joint Inversion Method

The medium is described with a three-dimensional mesh of equally spaced nodes of densities $\mathit{\rho }$. The density in the whole volume of interest is computed by trilinear interpolation of the densities assigned to the eight surrounding nodes. The gravimetric effect $\mathit{\gamma }$ produced by the volume of interest at the observation locations can be modeled using

where the matrix elements G_ij contain the modeled contribution of each node j of unitary density to each observation location i. The sensitivity matrix G is computed taking into account the topography (Coutant et al., 2012; Barnoud et al., 2016). The muographic data consist of densities averaged in conic bins of azimuth and elevation in the angle of view of the muon detector. These averaged densities μ are estimated from the flux of muons crossing the edifice (Cârloganu and the TOMUVOL Collaboration, 2018; Niess et al., 2018b) and are linearly related to the subsurface densities ρ via the sensitivity matrix M:

We use the joint linear Bayesian inversion method of Barnoud et al. (2019). To determine a constant offset between densities inferred from gravimetric data and muographic data, we adapt and implement the automatic least-squares determination proposed by Lelièvre et al. (2019). We therefore consider the following linear joint inversion problem:

where e is a vector of ones with the same length as μ_obs and c is a constant accounting for the potential density offset between the two datasets.

Following Tarantola (Tarantola, 2005), we use a Bayesian formalism where the solution minimizes the objective function

The first term of the objective function (Eq. 4) is the gravimetric data misfit

where ${\mathbf{C}}_{\gamma }$ is the gravimetric data covariance matrix. This matrix is usually diagonal as the gravimetric data errors are considered uncorrelated. The diagonal elements correspond to the data variances. The second term of the objective function (Eq. 4) is the muographic data misfit, accounting for a constant density offset c:

where C_µ is the muographic data misfit. Depending on the processing applied to obtain the averaged densities, this matrix might be diagonal (uncorrelated data) or it might be sparse and banded with non-diagonal elements. The third term of the objective function is the regularization term which consists of a distance to an a priori model with a smoothing applied:

The a priori density covariance matrix ${\mathbf{C}}_{\rho }$ is constructed with a Gaussian spatial covariance function:

where ${\sigma }_{\rho }$ is the a priori density standard deviation, D_ij is the distance between the two nodes i and j, and λ the spatial correlation length.

To obtain the density offset c, we seek the minimum of the objective function (Eq. 4) with respect to c. This is achieved by canceling its derivative with respect to c, leading to:

where we define the vector

In the particular case of a diagonal muographic data covariance matrix (uncorrelated data) with all the data variances being equal, ${\mathbf{C}}_{\mu }={\sigma }_{\mu }^2\mathbf{I}$ and the offset formulation (Eq. 9) simplifies to (Lelièvre et al., 2019):

where N_µ is the number of muographic data.

We rewrite the muographic data term of the objective function (Eq. 6) by replacing the expression of the constant offset c (Eq. 9):

where we define the matrix

Hence, solving the linear joint inverse problem with a Bayesian formalism and a constant offset (Eq. 3) is equivalent to minimizing the objective function

where A is the modified sensitivity matrix

d is the modified data vector

and C_d is the data covariance matrix

According to Tarantola (2005), minimizing the objective function (Eq. 14) leads to the estimated density

which is the center of the a posteriori Gaussian density distribution with the associated a posteriori covariance matrix (Tarantola, 2005)

The diagonal elements of ${\tilde{\mathbf{C}}}_{\rho }$ are the a posteriori variances of the estimated densities $\tilde{\mathit{\rho }}$ for each node. The modeled data computed from the estimated density $\tilde{\mathit{\rho }}$ are given by:

The inversion result (Eqs 18 and 19) is tuned by two regularization parameters, namely the a priori density standard deviation ${\sigma }_{\rho }$ and the spatial correlation length λ. We use the LOO criterion to estimate the optimal regularization parameters ${\sigma }_{\rho }$ and λ for the inversion as advocated by Barnoud et al. (2019). The inversion is performed as many times as the number of data N, one data being ignored each time. We retain the regularization parameters that minimize the error L:

3 Application to the Puy De Dôme Volcano

3.1 Data

The gravimetric data were acquired and processed by Portal et al. (2016). We use the 650 measurements covering the dome itself (Figure 2). The data processing includes the Earth tide correction using the Longman formula (Longman, 1959) and the instrumental drift. The free-air anomaly is computed using the Somigliana formula (Somigliana, 1930) and the WGS84 reference ellipsoid (Moritz, 2000). The Bouguer correction for unitary densities is computed (Coutant et al., 2012; Barnoud et al., 2016). The contribution of masses outside the area of interest and up to 160 km is corrected using a density of 2,670 kg/m³ (far-field). Close to the target, we determine the correction density with the Parasnis method (Parasnis, 1972). Using the full gravimetric dataset of Portal et al. (2016), we correct for the masses outside the volume of interest and up to about 4.5 km from the summit with a density of 2,420 kg/m³ (mid-field). For masses inside the volume of interest, we estimate the optimal density of correction of 1,890 kg/m³ (near-field). Finally, we remove the mean value of −15.97 mGal from the gravity data to account for the regional field component. The residual Bouguer anomaly, used as input data for the inversions, is shown on Figure 2. The uncertainties (i.e., data standard deviations), on the Bouguer anomaly is estimated to 0.039 mGal (Portal et al., 2016), including uncertainties from gravimetric and positioning measurements as well as from gravity corrections.

FIGURE 2. Residual Bouguer anomaly data of the Puy de Dôme volcano using correction densities of 2,670 kg/m³ (far-field), 2,420 kg/m³ (mid-field) and 1,890 kg/m³ (near-field). The gravimetric measurements are from Portal et al. (2016). The extent of this map is indicated by a blue square in Figure 1. The location of the muon detector installed at the Col de Ceyssat (CDC) is indicated (black square) as well as the angle of view of the detector (gray shaded area). The two black lines locate the slices of the inversion results shown in Figures 6, 7.

In this paper, we use as muographic data a preliminary result (Cârloganu and the TOMUVOL Collaboration, 2018; Niess et al., 2018a) from an acquisition campaign performed in 2015–2016, at a single view point. The map of average densities (Figure 3B) computed from the muon acquisition was provided by the TOMUVOL (TOmography avec des MUons atmosphériques de VOlcans) collaboration. The muon detector was installed at 1,080 m a.s.l. at the Col de Ceyssat (CDC), located at about 1,200 m from the summit of the Puy de Dôme edifice in the south-west direction (Figure 2). The TOMUVOL muon detector (Cârloganu et al., 2013; Ambrosino et al., 2015) is composed of four 1 m² layers of Gas Resistive Plate Chambers (GRPCs), similar to the ones developed for the detection of hadrons in the framework of the CALICE experiment (The CALICE Collaboration, 2016). The spacing between the outer layers is about 1.8 m. A 10 cm wall of lead is placed in the middle to deviate the trajectories of low energy particles. The campaign lasted for six months, resulting in three months of effective acquisition time (99.6 effective days of data tracking). The high resolution of the TOMUVOL detector allows tracking muons crossing up to about 1,600 m of rocks (Figure 3A). We only consider the muon path above 6° above the horizontal so that only the Puy de Dôme relief is crossed by the muons. Below this elevation, muons also cross the Petit Puy de Dôme, situated to the north-east of the Puy de Dôme (Figure 1), which is outside the region of interest in this study. We recall here the key steps of the process (Niess et al., 2016; Cârloganu and The TOMUVOL Collaboration, 2018; Niess et al., 2018a; Niess et al., 2018b; Niess et al., 2020). The muon tracks are summed up in bins of 1° by 1° in azimuth and elevation. A kernel method is used to increase the muon count per bin, hence decrease the statistical uncertainty in the estimate of the transmitted flux of muons (Cârloganu et al., 2016; Cârloganu and The TOMUVOL Collaboration, 2018). The observed rate of muons is the convolution of the transmitted flux with the effective surface of the detector. For each bin, the density best fitting the observed muon rate is retained. The muographic data used for the inversion therefore consist of a map of densities averaged along the lines of sight (Figure 3B). The overall density of the Puy de Dôme, averaged and weighted by the rock thickness (Figures 3A,B), is 1,565 kg/m³, which is 325 kg/m³ lower than the averaged density estimated from the gravimetric data. This difference can partly be explained by the fact that the muographic data do not encompass the dense lava flows to the west and to the south-east of the dome (Figure 2). The statistical and systematic errors are estimated using a bootstrap procedure (Niess et al., 2018b; Niess et al., 2018a; Niess et al., 2020). The statistical errors correspond to the measurement standard deviations. The estimated standard deviations are below 100 kg/m³ in most of the dome, with larger values close to the surface and to the south-east (Figure 3C). The systematic errors are the biases (i.e., offsets) between the measured value and the true value. The estimated biases are about one order of magnitude lower than the standard deviations (Figure 3D), with values lower than 20 kg/m³ for most of the bins. The bias is therefore underestimated compared to the 325 kg/m³ offset of the averaged density with respect to the gravimetric data.

FIGURE 3. Preliminary muographic data of the Puy de Dôme volcano as observed from the Col de Ceyssat with the TOMUVOL detector and associated errors (Cârloganu and the TOMUVOL collaboration, 2018; Niess et al., 2018a). The location of the muon detector is indicated in Figure 2. (A) Thickness of rock along the muon trajectories. (B) Averaged densities estimated from the muon count. (C) Statistical error (i.e., standard deviation). (D) Systematic error (i.e., bias).

In addition to these data, density measurements have been performed on a collection of 45 rock samples from the Puy de Dôme, including the ones from Portal (Portal, 2015). The densities for dry samples range from 1,300 kg/m³ to 2,080 kg/m³ with an average of 1,839 kg/m³ (Figure 4), in accordance to the mean density estimated from the gravimetric data. When water-saturated, the average density reaches 2,044 kg/m³ with densities up to 2,200 kg/m³. The trachytes that constitute the Puy de Dôme therefore present densities lower than usual massive trachytes, such as the ones from New Zealand for instance (Tenzer et al., 2011). This low density can be explained by the high porosity of the samples (Boudon et al., 2015; Portal, 2015).

FIGURE 4. Densities measured on dry and water-saturated samples of trachyte from the Puy de Dôme volcano. The dashed histogram corresponds to the normalized distribution of the 3D density model resulting from the joint inversion (Figure 7).

3.2 Inversion Results and Discussion

We describe the volume using a mesh of 85 × 85 × 29 density nodes with a regular spacing of 25 m. This mesh extends up to 1,000 m away from the Puy de Dôme summit and up to an elevation of 800 m, i.e., more than 600 m below the summit. The Chaîne des Puys is covered by a Lidar Digital Elevation Model (DEM) with a 50 cm spatial resolution and a 10 cm precision. We use this DEM for the computation of the gravimetric and muographic data sensitivity matrices. We present here the inversion of the gravimetric data alone, the inversion of the muographic data alone and the joint inversion of the gravimetric and muographic data with offset determination. The LOO criterion is used to determine the optimal set of correlation length λ and a priori density standard deviation ${\sigma }_{\rho }$ for each inversion (Figure 5). Synthetic tests presented by Barnoud et al. (2019) show the impact of these two regularization parameters on the density model and that the LOO is a robust criterion to determine them. This is particularly important in the present case where only one muographic view point is available. The density models resulting from the inversions using the optimal regularization parameters are shown in Figure 6 with two perpendicular slices extracted from the 3D resulting density models (see Figure 2 for locations). We show the a posteriori mean density, the a posteriori standard deviation, and a random realization within the a posteriori distribution.

FIGURE 5. Maps of the Leave One Out (LOO) criterion L (Eq. 21) for (A) the gravimetric independent inversion, (B) the muographic independent inversion and (C) the joint inversion. White dots indicate minimal values of L obtained for the optimal sets of a priori density standard deviations ${\sigma }_{\rho }$ and spatial correlation lengths λ according to the Leave One Out criterion.

FIGURE 6. Inversion results: mean of the a posteriori density distribution (blue to red color scale), standard deviation (gray scale, the darker the better resolved) and random realization of the posterior distribution (blue to red color scale). The blue to red color scales of all panels have the same range of ±400 kg/m³. (A) Independent inversion of the gravimetric data. The density color scales are centered on 1,890 kg/m³, corresponding to the optimal density correction determined with the Parasnis method from the gravimetric data. (B) Independent inversion of the muographic data. The density color scales are centered on 1,565 kg/m³, corresponding to the weighted averaged of densities from the muographic data. Note that the gray colorscale for the a posteriori standard deviation is different than for the gravimetric and joint inversion results. (C) Joint inversion. The density color scales are centered on 1,890 kg/m³, as for the gravimetric inversion, to ease the comparison between the results. The locations of the SSW-NNE and WNW-ESE sections are indicated on Figure 2.

The optimal inversion result for the gravimetric inversion is obtained for $\left({\sigma }_{\rho },\lambda \right)=\left(50\text{ kg}/{\text{m}}^3,150\text{ m}\right)$ (Figure 5A). As indicated by the a posteriori standard deviation (Figure 6A), the corresponding density model is well resolved laterally and close to the surface while the resolution decreases in the core of the dome and at depth, which is typical of the inversion of gravimetric data.

Here, the inversion of the muographic data alone is not meaningful by itself, as one cannot retrieve a 3D distribution from a single view point, but it is shown for the purpose of illustrating the specificities of the method. The optimal muographic inversion result is obtained for $\left({\sigma }_{\rho },\lambda \right)=\left(\mathrm{2,100}\text{kg}/{\text{m}}^3,40\text{m}\right)$, where the a priori density standard deviation is poorly constrained as indicated by the horizontal elongation of the area of lower values of L (Eq. 21) on the map (Figure 5B). However, with further improvements of physical simulations in muography, we expect muographic data to become closer to the true absolute density in the future. Therefore, the a priori density standard deviation could rather be determined from density distribution of the data themselves so that only the spatial correlation length could be determined with an LOO criterion. With a single view point, the inversion spreads out the average densities along the conical bins, resulting in spread out densities in the SSW-NNE slice parallel to the view direction of the detector (Figure 6B, top) and in densities looking similar to the initial map (Figure 3B) in the perpendicular WNW-ESE slice (Figure 6B, bottom).

The optimal regularization parameters for the joint inversion are $\left({\sigma }_{\rho },\lambda \right)=\left(50\text{ kg}/{\text{m}}^3,150\text{ m}\right)$ (Figure 5C), i.e., the same as the ones obtained for the gravimetric inversion. In our case study, the joint inversion is driven by the gravimetric data due to the dense coverage of these data and to the balance between the uncertainties of the gravimetric data and of the muographic data. As indicated by the mean density models (Figure 6C), the information contained in the muographic data slightly improves the imaging of the core structures of the dome compared to the gravimetric independent inversion. We can note that a low density structure at around 1,000 m of elevation on the WNW-ESE profile appears in the joint inversion and that the high density core below the summit is better delineated at depth. However, standard deviations indicate that this improvement is very subtle. This can be attributed to the availability of only one muographic view point. This is supported by tests on synthetic data with the topography of the Puy de Dôme (Barnoud et al., 2019; Lelièvre et al., 2019). The addition of view points greatly helps constraining the model close to the surface (Barnoud et al., 2019), in this way improving the imaging at depth (Lelièvre et al., 2019).

The automatic least-squares determination of a constant offset allows to easily invert residual Bouguer anomaly data jointly with averaged densities from muographic data. The constant offset determined in the joint inversion from Eq. 9 is c = 1,650 kg/m³. We recall that the residual Bouguer anomalies have been obtained by substracting a constant density of 1,890 kg/m³, while the muographic data correspond to absolute densities. Consequently, the offset c corresponds to a mean difference of 240 kg/m³ between the muographic and the gravimetric data. This is lower than the difference of 325 kg/m³ observed from the data indicating that the joint inversion helps the resulting density model to be consistent with the two datasets. Note that we have assumed here a constant density offset, taken into account using a mathematical model. To take into account a non constant density offset, a comprehensive physical model would be necessary to fully estimate the bias contained in the data. To ease the comparison with the gravimetric inversion result, we choose the average gravity density of 1,890 kg/m³ as central value for the joint inversion result. The statistical distribution of the densities resulting from the joint inversion is shown in Figure 4. This distribution is closer to the distribution of the dry rock sample densities (average of 1,839 kg/m³) than to the distribution of the water-saturated rock sample densities (averaged of 2,044 kg/m³). Even if the joint inversion result was centered on the average density from the muographic data (1,565 kg/m³) or on an intermediate density value, the whole densities would decrease by a few hundreds of kg/m³, still more consistent with the distribution of the dry rock samples than with the distribution of the water-saturated samples. This supports the idea of a dry dome, highly permeable, so that water quickly flows at depth, in accordance with the literature (Lecoq, 1831; Boudon et al., 2015).

The 3D density model resulting from the joint inversion allows to identify several key structures of the Puy de Dôme volcano (Figures 6C, 7). In addition to the SSW-NNE and WNW-ESE cross-sections of densities and associated standard-deviations (Figure 6C), the 3D density structure of the Puy de Dôme is shown using isosurfaces of 1,690 kg/m³ and 2,090 kg/m³, corresponding to the averaged density of 1,890 kg/m³ ± 200 kg/m³ (Figure 7). These isosurfaces highlight low and high density zones inside the edifice but should not be considered as sharp limitations of the structures. Instead, the extent of the structures (especially the depth extent) is better shown by the random realization of densities from the a posteriori distribution (Figure 6C). We interpret the main identified structures (Figure 7) with regards to the geological map (Figure 1) and the available literature on the area (Miallier et al., 2010; Portal et al., 2016; Boivin et al., 2017; Portal et al., 2019; Deniel et al., 2019). Three main low density anomalies are identified (noted 1, 2 and 3 on Figure 7): 1) to the north-east of the Puy de Dôme, toward the Petit Puy de Dôme edifice, 2) below the north-west slope of the Puy de Dôme at depth, in the Corneboeufs area (Miallier et al., 2010) and 3) at the location of Puy Redon. These three low density anomalies correspond to scoria and cinder Strombolian cones erected before the construction of the Puy de Dôme (Boivin et al., 2017). To the north-west of the Puy de Dôme, higher density anomalies are found in the area covered by basaltic rocks (Boivin et al., 2017) (lava flow noted 4 on Figure 7) which are known to be denser than trachytes. The model also highlights the trachytic core of the dome (noted 5 on Figure 7), denser than the surrounding trachytic breccia, and an associated feeding conduit located to the NNW of the summit. This structure is typical of volcanic domes, with a compact core surrounded by a clastic talus resulting from rockfalls (Hale et al., 2009), as at Soufriere Hills volcano in Montserrat, Lesser Antilles (Sparks and Young, 2002). Another example is the Showa-Shinzan lava dome, Usu, Japan, where the joint gravimetric and muographic inversion results of Nishiyama et al. (2017b) show a massive cylindrical lava block of about 300 m in diameter. We also identify similar higher densities in the area of the Cratère Kilian (noted 6 on Figure 7) which presents trachyte rocks similar to those of the Puy de Dôme (Boivin et al., 2017). Finally, the model brings out a linear high density structure to the south-east, at the bottom of the edifice (noted 7 on Figure 7). This structure is oriented in a north 30° direction and might correspond to an old underlying trachy-basaltic or basaltic lava flow (Figure 1) (Boivin et al., 2017). As shown by the standard-deviations in Figure 6C, the whole density model is better resolved in the first 200–300 m below the surface. The dense trachytic core of the Puy de Dôme 5) and the low density material of the Puy Redon 3) are better delineated at depth thanks to the addition of the muographic data, as indicated by the standard-deviations show in Figures 6A,C. However, structures such as the high density body of the Cratère Kilian 6) or the buried lava flow 7) are located outside the angle of view of the muon detector and are mostly constrained by the gravimetric data. These structures are therefore well constrained laterally, but their depth extents are poorly resolved. Note that, in general, even when some structures are situated outside the angle of view of the muon detector, their imaging is likely to be improved by the joint inversion. Indeed, the gravimetric and muographic datasets are linked in the inversion process, so that the model improvement is global. Nevertheless, the gain in resolution depends on the distribution of gravimetric data (lateral extent and density of coverage), hence on the resolving power of the gravimetric data at depth.

FIGURE 7. 3D density model of the Puy de Dôme volcano resulting from the joint inversion of muographic and gravimetric data, with geological map wrapped on topography. The lateral extent of the model is indicated by the blue square on Figure 1. The averaged density is 1,890 kg/m³. Iso densities of 1,690 and 2,090 kg/m³, corresponding to the averaged density ±200 kg/m³, are chosen to highlight low and high density zones in blue and red respectively. The SSW-NNE and WNW-ESE vertical sections are the same as in Figure 6 (locations indicated on Figure 2). The numbered features indicate the structures discussed in the article. The black square and dashed lines indicate the location of the muon detector at Col de Ceyssat (CDC) and its angle of view.

4 Conclusion

We have presented a robust Bayesian joint inversion method developed in order to reconstruct the 3D density structure of a 1,600 m large volcanic dome. The determination of the two optimal inversion regularization parameters with an LOO algorithm overcomes the artifacts linked to the acquisition geometry. The automatic determination of a constant density offset accounts for a possible bias between the density inferred from the two datasets.

For the case study of the Puy de Dôme volcano, the final 3D density structure is free of data acquisition artifacts and the recovered structures are in accordance with previous geological and geophysical studies. In particular, our results confirm that the Puy de Dôme is dry and permeable, with low density porous trachytic rocks. The muographic data acquired at one location bring a slight improvement to the model compared to the gravimetric data alone. The trachytic core of the dome and its surrounding less dense talus are better delineated. The recovered structure is typical of a lava dome such as Soufriere Hills, Montserrat. The resolution improvement could be amplified thanks to muon acquisitions from complementary view points as shown with synthetic tests (Barnoud et al., 2019; Lelièvre et al., 2019).

The joint inversion method presented in this paper, combined with improved high resolution muon detectors and multiple acquisition view points, has the potential to better image the structure and plumbing system of kilometer size volcanic edifices. Beyond imaging, the method could be used to monitor slow processes including fracturing and flanks instabilities.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author Contributions

AB developed the codes for the inversion of the muographic and gravimetric data, produced the results and wrote the paper. VC and PGL participated in the development of the inversion methodology. AP and LG acquired the gravimetric data. AB and AP processed the gravimetric data. PL participated in the acquisition of the muographic data. AP, PL, and PB provided the density measurements of rock samples. PB provided the geological map. All authors participated in the interpretation of the results and reviewed the manuscript.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We warmly thank Valentin Niess, Eve Le Ménédeu, and Didier Miallier for fruitful exchanges, collaboration and support. The authors thank the TOMUVOL collaboration for providing the preliminary muographic data (www.obs.univ-bpclermont.fr/tomuvol/). We thank Olivier Coutant for the forward gravimetric computation code. We thank Andrew Harris for discussions on lava dome structures. We thank Yosuke Aoki and Valerio Acocella from the editorial board for their help in the review process and two anonymous reviewers for their comments and corrections that helped improve the article. The LIDAR data used in this study derive from a collective project driven by the Center Régional Auvergnat de l’Information Géographique (CRAIG), financially supported by the Conseil Départemental du Puy-de-Dôme, the European Regional Development Fund and the University of Clermont-Ferrand. Computations have been performed on the supercomputer facilities of the Mésocentre Clermont Auvergne University. This research was financed by the French Government Laboratory of Excellence initiative n° ANR-10-LABX-0006, the Region Auvergne and the European Regional Development Fund. This is Laboratory of Excellence ClerVolc contribution number 430.

References