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

Anne Barnoud, Valérie Cayol, Peter G Lelièvre, Angélie Portal, Philippe Labazuy, Pierre Boivin, Lydie Gailler

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

Anne Barnoud, Valérie Cayol, Peter G Lelièvre, Angélie Portal, Philippe Labazuy, et al.. Robust Bayesian Joint Inversion of Gravimetric and Muographic Data for the Density Imaging of the Puy de Dôme Volcano (France). Frontiers in Earth Science, 2021, 8, 10.3389/feart.2020.575842. hal-03108803

HAL Id: hal-03108803

https://uca.hal.science/hal-03108803

Submitted on 13 Jan 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

Check for
updates

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

Anne Barnoud¹,²*, Valérie Cayol², Peter G. Lelièvre³, Angélie Portal²,⁴, Philippe Labazuy², Pierre Boivin² and Lydie Gailler²

¹LPC, Université Clermont Auvergne, CNRS/IN2P3, Clermont-Ferrand, France, ²Laboratoire Magmas et Volcans, Université Clermont Auvergne, CNRS, IRD, OPGC, Clermont-Ferrand, France, ³Department of Mathematics and Computer Science, Mount Allison University, Sackville, NB, Canada, ⁴BRGM, DRP/IGT, Orléans, France

Imaging the internal structure of volcanoes helps highlighting magma pathways and monitoring potential structural weaknesses. We jointly invert gravimetric and muographic data to determine the most precise image of the 3D density structure of the Puy de Dôme volcano (Chaîne des Puys, France) ever obtained. With rock thickness of up to 1,600 m along the muon lines of sight, it is, to our knowledge, the largest volcano ever imaged by combining muography and gravimetry. The inversion of gravimetric data is an ill-posed problem with a non-unique solution and a sensitivity rapidly decreasing with depth. Muography has the potential to constrain the absolute density of the studied structures but the use of the method is limited by the possible number of acquisition view points, by the long acquisition duration and by the noise contained in the data. To take advantage of both types of data in a joint inversion scheme, we develop a robust method adapted to the specificities of both the gravimetric and muographic data. Our method is based on a Bayesian formalism. It includes a smoothing relying on two regularization parameters (an a priori density standard deviation and an isotropic correlation length) which are automatically determined using a leave one out criterion. This smoothing overcomes artifacts linked to the data acquisition geometry of each dataset. A possible constant density offset between both datasets is also determined by least-squares. The potential of the method is shown using the Puy de Dôme volcano as case study as high quality gravimetric and muographic data are both available. Our results show that the dome is dry and permeable. Thanks to the muographic data, we better delineate a trachytic dense core surrounded by a less dense talus.

**Keywords:** joint inversion method, Bayesian inversion, gravimetry, muography, volcanology, lava dome, Puy de Dôme

## 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

**OPEN ACCESS**

**Edited by:**
Yosuke Aoki,
The University of Tokyo, Japan

**Reviewed by:**
Mousumi Roy,
University of New Mexico,
United States
László Oláh,
The University of Tokyo, Japan

**\*Correspondence:**
Anne Barnoud
anne.barnoud@gmail.com

**Specialty section:**
This article was submitted to Volcanology, a section of the journal Frontiers in Earth Science

**Received:** 24 June 2020
**Accepted:** 09 October 2020
**Published:** 11 January 2021

**Citation:**
Barnoud A, Cayol V, Lelièvre PG, Portal A, Labazuy P, Boivin P and Gailler L (2021) Robust Bayesian Joint Inversion of Gravimetric and Muographic Data for the Density Imaging of the Puy de Dôme Volcano (France). Front. Earth Sci. 8:575842
doi: 10.3389/feart.2020.575842
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.

## 2 JOINT INVERSION METHOD

The medium is described with a three-dimensional mesh of equally spaced nodes of densities $\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 $\gamma$ produced by the volume of interest at the observation locations can be modeled using

$$G\rho = \gamma \qquad (1)$$

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 $\mu$ are estimated from the flux of muons
<figure><img src="image_4.png" /><figcaption>**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).</figcaption></figure>

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:

$$ M\rho = \mu. \qquad (2) $$

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:

$$ \begin{bmatrix} G \\ M \end{bmatrix} \rho = \begin{bmatrix} \gamma_{\text{obs}} \\ \mu_{\text{obs}} - c e \end{bmatrix}. \qquad (3) $$

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

$$ \varphi(\rho) = \|G\rho - \gamma_{\text{obs}}\|_{\gamma}^{2} + \|M\rho - \mu_{\text{obs}} + c e\|_{\mu}^{2} + \|\rho - \rho_{\text{prior}}\|_{\rho}^{2}. \qquad (4) $$

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

$$ \|G\rho - \gamma_{\text{obs}}\|_{\gamma}^{2} = (G\rho - \gamma_{\text{obs}})^{T} C_{\gamma}^{-1} (G\rho - \gamma_{\text{obs}}) \qquad (5) $$

where Cγ 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:

$$ \|M\rho - \mu_{\text{obs}} + c e\|_{\mu}^{2} = (M\rho - \mu_{\text{obs}} + c e)^{T} C_{\mu}^{-1} (M\rho - \mu_{\text{obs}} + c e) \qquad (6) $$

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:

$$ \|\rho - \rho_{\text{prior}}\|_{\rho}^{2} = (\rho - \rho_{\text{prior}})^{T} C_{\rho}^{-1} (\rho - \rho_{\text{prior}}). \qquad (7) $$

The a priori density covariance matrix Cρ is constructed with a Gaussian spatial covariance function:

$$ C_{\rho,ij} = \sigma_{\rho}^{2} \exp(-D_{ij}^{2} / \lambda^{2}) \qquad (8) $$

where σρ is the a priori density standard deviation, Dij 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:

$$ c = - \frac{e^{T} C_{\mu}^{-1} (M\rho - \mu_{\text{obs}})}{e^{T} C_{\mu}^{-1} e} = b^{T} (M\rho - \mu_{\text{obs}}) \qquad (9) $$
where we define the vector

$$b = - \frac{C_{\mu}^{-1} e}{e^{T} C_{\mu}^{-1} e} \qquad (10)$$

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

$$c = - \frac{e^{T} (M\rho - \mu_{\text{obs}})}{N_{\mu}} \qquad (11)$$

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):

$$\begin{aligned} & \|M\rho - \mu_{\text{obs}} + ce\|_{\mu}^{2} \\ &= [(I + eb^{T})(M\rho - \mu_{\text{obs}})]^{T} C_{\mu}^{-1} [(I + eb^{T})(M\rho - \mu_{\text{obs}})] \\ &= \|BM\rho - B\mu_{\text{obs}}\|_{\mu}^{2} \end{aligned} \qquad (12)$$

where we define the matrix

$$B = I + eb^{T} = I - \frac{ee^{T} C_{\mu}^{-1}}{e^{T} C_{\mu}^{-1} e} \qquad (13)$$

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

$$\begin{aligned} \varphi(\rho) &= \|A\rho - d\|_{d}^{2} + \|\rho - \rho_{\text{prior}}\|_{\rho}^{2} \\ &= (A\rho - d)^{T} C_{d}^{-1} (A\rho - d) + (\rho - \rho_{\text{prior}})^{T} C_{\rho}^{-1} (\rho - \rho_{\text{prior}}) \end{aligned} \qquad (14)$$

where A is the modified sensitivity matrix

$$A = \begin{bmatrix} G \\ BM \end{bmatrix}, \qquad (15)$$

d is the modified data vector

$$d = \begin{bmatrix} \gamma_{\text{obs}} \\ B\mu_{\text{obs}} \end{bmatrix}, \qquad (16)$$

and $C_d$ is the data covariance matrix

$$C_d = \begin{bmatrix} C_{\gamma} & 0 \\ 0 & C_{\mu} \end{bmatrix}. \qquad (17)$$

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

$$\tilde{\rho} = \rho_{\text{prior}} + C_{\rho} A^{T} (A C_{\rho} A^{T} + C_{d})^{-1} (d - A \rho_{\text{prior}}) \qquad (18)$$

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

$$\tilde{C}_{\rho} = C_{\rho} - C_{\rho} A^{T} (A C_{\rho} A^{T} + C_{d})^{-1} A C_{\rho}. \qquad (19)$$

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

$$d_{\text{cal}} = \begin{bmatrix} G\tilde{\rho} \\ M\tilde{\rho} + ce \end{bmatrix}. \qquad (20)$$

The inversion result (Eqs 18 and 19) is tuned by two regularization parameters, namely the a priori density standard deviation $σ_ρ$ and the spatial correlation length $λ$. We use the LOO criterion to estimate the optimal regularization parameters $σ_ρ$ 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:

$$L(\sigma_{\rho}, \lambda) = \frac{1}{N} \sum_{l=1}^{N} \frac{[(A\tilde{\rho})_{l} - d_{l}]^{2}}{\sigma_{d,l}^{2}}. \qquad (21)$$

# 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.

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
<figure><img src="image_5.png" /></figure>

**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**.

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

**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).

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
<figure><img src="image_10.png" /><figcaption>**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**).</figcaption></figure>

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.

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).

## 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 $\lambda$ and a priori density standard deviation $\sigma_p$ 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.

The optimal inversion result for the gravimetric inversion is obtained for ($\sigma_p, \lambda$) = (50 kg/m³, 150 m) (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
<figure><img src="image_11.png" /><figcaption>**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 *σ*<sub>p</sub> and spatial correlation lengths λ according to the Leave One Out criterion.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**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**.</figcaption></figure>

illustrating the specificities of the method. The optimal muographic inversion result is obtained for ($\sigma_p, \lambda$) = (2, 100 kg/m³, 40 m), 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.
<figure><img src="image_13.png" /></figure>

**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.

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 ($\sigma_\rho, \lambda$) = (50 kg/m³, 150 m) (**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 subtracting 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.

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; Lelè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.

REFERENCES

Alvarez, L. W., Anderson, J. A., Bedwei, F., Burkhard, J., Fakhry, A., Girgis, A., et al. (1970). Search for hidden chambers in the pyramids. *Science* 167, 832–839.
doi:10.1126/science.167.3919.832

Ambrosino, F., Anastasio, A., Bross, A., Béné, S., Boivin, P., Bonechi, L., et al. (2015). Joint measurement of the atmospheric muon flux through the puy de Dôme volcano with plastic scintillators and resistive plate chambers detectors. *J. Geophys. Res.* 120, 7290–7307. doi:10.1002/2015JB011969

Barnoud, A., Cayol, V., Niess, V., Cárloganu, C., Lelièvre, P., Labazuy, P., et al. (2019). Bayesian joint monographic and gravimetric inversion applied to volcanoes. *Geophys. J. Int.* 218, 2179–2194. doi:10.1093/gji/ggz300

Barnoud, A., Coutant, O., Bouligand, C., Gunawan, H., and Deroussi, S. (2016). 3-D linear inversion of gravity data: method and application to Basse-Terre volcanic island, Guadeloupe, Lesser Antilles. *Geophys. J. Int.* 205, 562–574. doi:10.1093/gji/ggw030

Boivin, P., Besson, J.-C., Briot, D., Deniel, C., Gourgaud, A., Labazuy, P., et al. (2017). Volcanology of the Chaîne des Puys Massif Central Français. *Aydat, France: Parc Naturel Régional des Volcans d'Auvergne, 199.*

Boudon, G., Balcone-Boissard, H., Villemant, B., and Morgan, D. J. (2015). What factors control superficial lava dome explosivity? *Sci. Rep.* 5, 14551. doi:10.1038/srep14551

Cárloganu, C.; The TOMUVOL collaboration (2018). "Density imaging of puy de Dôme volcano with atmospheric muons in French Massif central as a case study for volcano muography," in EGU general assembly conference abstracts, vol. 20, 18420.

Cárloganu, C., and Le Ménédeu, E.; The TOMOVOL Team (2016). "Density imaging of puy de Dôme volcano with atmospheric muons in French Massif Central as a case study for volcano muography," in Geophysical research abstracts, EGU General Assembly, Vol. 18.

Cárloganu, C., Niess, V., Béné, S., Busato, E., Dupieux, P., Fehr, F., et al. (2013). Towards a muon radiography of the Puy de Dôme. *Geosci. Instrum. Methods Data Syst.* 2, 55–60. doi:10.5194/gi-2-55-2013

Camacho, A. G., Montesinos, F. G., and Vieira, R. (1997). A three-dimensional gravity inversion applied to São Miguel island (Azores). *J. Geophys. Res.* 102, 7717–7730. doi:10.1029/96JB03667

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 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.

Cella, F., Fedi, M., Florio, G., Grimaldi, M., and Rapolla, A. (2007). Shallow structure of the Somma-Vesuvius volcano from 3D inversion of gravity data. *J. Volcanol. Geoth. Res.* 161, 303–317. doi:10.1016/j.jvolgeores.2006.12.013

Cosburn, K., Roy, M., Guardincerri, E., and Rowe, C. (2019). Joint inversion of gravity with cosmic ray muon data at a well-characterized site for shallow subsurface density prediction. *Geophys. J. Int.* 217, 1988–2002. doi:10.1093/gji/ggz127

Coutant, O., Bernard, M. L., Beauducel, F., Nicollin, F., Bouin, M. P., and Roussel, S. (2012). Joint inversion of P-wave velocity and density, application to La Soufrière of Guadeloupe hydrothermal system. *Geophys. J. Int.* 191, 723–742. doi:10.1111/j.1365-246X.2012.05644.x

Davis, K., and Oldenburg, D. W. (2012). *Joint 3D of muon tomography and gravity data to recover density*. CSIRO, 1–4.

Davis, K., Oldenburg, D. W., Kaminski, V., Pilkington, M., Bryman, D., Bueno, J., et al. (2011). "Joint 3D inversion of muon tomography and gravity data," in International workshop on gravity, electrical & magnetic methods and their applications 2011: international workshop on gravity, electrical & magnetic methods and their applications, October 10–13, 2011 (Beijing, China: Society of Exploration Geophysicists), 81.

Deniel, C., Boivin, P., Miallier, D., and Gerbe, M. C. (2019). Multi-stage growth of the trachytic lava dome of the Puy de Dôme (Chaîne des Puys, France). Field, geomorphological and petro-geochemical evidence. *J. Volcanol. Geoth. Res.* 396, 106749. doi:10.1016/j.jvolgeores.2019.106749

Gómez, H., Gibert, D., Goy, C., Jourde, K., Karyotakis, Y., Katsanevas, S., et al. (2017). Forward scattering effects on muon imaging. *J. Instrum.* 12, P12018. doi:10.1088/1748-0221/12/12/P12018

Groom, D. E., Mokhov, N. V., and Striganov, S. I. (2001). Muon stopping power and range tables 10 MeV–100 TeV. *Atomic Data Nucl. Data Tables* 78, 183–356. doi:10.1006/adnd.2001.0861

Hale, A. J., Calder, E. S., Wadge, G., Loughlin, S. C., and Ryan, G. A. (2009). Modelling the lava dome extruded at Soufrière Hills volcano, Montserrat, August 2005–May 2006: Part I: dome shape and internal structure. *J. Volcanol. Geoth. Res.* 187, 53–68. doi:10.1016/j.jvolgeores.2009.08.023

Jourde, K., Gibert, D., Marteau, J., de Bremond d'Ars, J., Gardien, S., Girerd, C., et al. (2013). Experimental detection of upward going cosmic particles and consequences for correction of density radiography of volcanoes. *Geophys. Res. Lett.* 40, 6334–6339. doi:10.1002/2013GL058357
Kusagaya, T., and Tanaka, H. K. M. (2015). Development of the very long-range cosmic-ray muon radiographic imaging technique to explore the internal structure of an erupting volcano, Shinmoe-dake, Japan. *Geosci. Instrum. Methods Data Sys.* 4, 215–226. doi:10.5194/gi-4-215-2015

Lecoq, H. (1831). *Itinéraire de Clermont au Puy de Dôme, ou description de cette montage et de la vallée de Royat et Fontanat*. Clermont-Ferrand, France: Imprimerie de Thibaud-Landriot.

Lelièvre, P. G., Barnoud, A., Niess, V., Cârloganu, C., Cayol, V., and Farquharson, C. G. (2019). Joint inversion methods with relative density offset correction for muon tomography and gravity data, with application to volcano imaging. *Geophys. J. Int.* 218, 1685–1701. doi:10.1093/gji/ggz251

Lesparre, N., Gibert, D., and Marteau, J. (2012). Density muon radiography of La Soufrière of Guadeloupe volcano: comparison with geological, electrical resistivity and gravity data. *Geophys. J. Int.* 185, 1–14. doi:10.1111/j.1365-246X.2012.05546.x

Linde, N., Baron, L., Finizola, A., Revil, A., Muccini, F., Cocchi, L., et al. (2014). 3-D density structure and geological evolution of Stromboli volcano (Aeolian Islands, Italy) inferred from land-based and sea-surface gravity data. *J. Volcanol. Geoth. Res.* 273, 58–69. doi:10.1016/j.jvolgeores.2014.01.006

Longman, I. M. (1959). Formula for computing the tidal accelerations due to the Moon and the Sun. *J. Geophys. Res.* 64, 2351–2355. doi:10.1029/JZ064i012p02351

Miallier, D., Boivin, P., Deniel, C., Bourgaud, A., Lanos, P., Sforna, M., et al. (2010). The ultimate summit eruption of Puy de Dôme volcano (Chaîne des Puys, French Massif Central) about 10,700 years ago. *Compt. Rendus Geosci.* 342, 847–854. doi:10.1016/j.crte.2010.09.004

Moritz, H. (2000). Geodetic reference system 1980. *J. Geodes.* 74, 128–140. doi:10.1007/S001900050278

Nagamine, K. (2003). *Introductory muon science*. Cambridge, UK: Cambridge University Press, 208.

Nagamine, K., Iwasaki, M., Shimomura, K., and Ishida, K. (1995). Method of probing inner-structure of geophysical substance with the horizontal cosmic-ray muons and possible application to volcanic eruption prediction. *Nucl. Instrum. Methods Phys. Res. Sect. A Accel. Spectrom. Detect. Assoc. Equip.* 356, 585–595. doi:10.1016/0168-9002(94)01169-9

Niess, V., Barnoud, A., Cârloganu, C., and Le Ménédeu, E. (2018a). “Methods and tools for transmission muography,” in EGU general assembly conference abstracts, vol. 20, 18110

Niess, V., Barnoud, A., Cârloganu, C., and Le Ménédeu, E. (2018b). Backward monte-carlo applied to muon transport. *Comput. Phys. Commun.* 229, 54–67. doi:10.1016/j.cpc.2018.04.001

Niess, V., Barnoud, A., Cârloganu, C., and Martineau-Huynh, O. (2020). TURTLE: a C library for an optimistic stepping through a topography. *Comput. Phys. Commun.* 247, 106952. doi:10.1016/j.cpc.2019.106952

Niess, V., Barnoud, A., Le Ménédeu, E., Cayol, V., and Cârloganu, C. (2016). “Density imaging of Puy de Dôme volcano by joint inversion of muographic and gravimetric data,” in AGU fall meeting (AGU).

Nishiyama, R., Ariga, A., Ariga, T., Käser, S., Lechmann, A., Mair, D., et al. (2017a). First measurement of ice-bedrock interface of alpine glaciers by cosmic muon radiography. *Geophys. Res. Lett.* 44, 6244–6251. doi:10.1002/2017gl073599

Nishiyama, R., Miyamoto, S., and Naganawa, N. (2014a). Experimental study of source of background noise in muon radiography using emulsion film detectors. *Geosci. Instrum. Method. Data Syst.* 3, 29–39. doi:10.5194/gi-3-29-2014

Nishiyama, R., Miyamoto, S., Okubo, S., Oshima, H., and Maekawa, T. (2017b). 3D density modeling with gravity and muon-radiographic observations in Showa-Shinzan lava dome, Usu, Japan. *Pure Appl. Geophys.* 174, 1061–1070. doi:10.1007/s00024-016-1430-9

Nishiyama, R., Taketa, A., Miyamoto, S., and Kasahara, K. (2016). Monte Carlo simulation for background study of geophysical inspection with cosmic-ray muons. *Geophys. J. Int.* 206, 1039–1050. doi:10.1093/gji/ggw191

Nishiyama, R., Tanaka, Y., Okubo, S., Oshima, H., Tanaka, H. K. M., and Maekawa, T. (2014b). Integrated processing of muon radiography and gravity anomaly data toward the realization of high-resolution 3-D density structural analysis of volcanoes: case study of Showa-Shinzan lava dome, Usu, Japan. *J. Geophys. Res. Solid Earth* 119, 699–710. doi:10.1002/2013JB010234

Okubo, S., and Tanaka, H. K. M. (2012). Imaging the density profile of a volcano interior with cosmic-ray muon radiography combined with classical gravimetry. *Meas. Sci. Technol.* 23, 042001. doi:10.1088/0957-0233/23/4/042001

Oláh, L., Tanaka, H. K. M., Ohminato, T., and Varga, D. (2018). High-definition and low-noise muography of the Sakurajima volcano with gaseous tracking detectors. *Sci. Rep.* 8, 3207. doi:10.1038/s41598-018-21423-9

Oláh, L., and Varga, D. (2017). Investigation of soft component in cosmic ray detection. *Astropart. Phys.* 93, 17–27. doi:10.1016/j.astropartphys.2017.06.002

Onizawa, S. y., Mikada, H., Watanabe, H., and Sakashita, S. (2002). A method for simultaneous velocity and density inversion and its application to exploration of subsurface structure beneath Izu-Oshima volcano, Japan. *Earth Planets Space* 54, 803–817. doi:10.1186/BF03352074

Parasnis, D. S. (1972). *Principles of applied geophysics*. London, UK: Chapman and Hall, 402.

Portal, A. (2015). Étude géophysique de la structure interne d'un dôme volcanique : le Puy de Dôme et son environnement (Chaîne des Puys, France). PhD thesis. Clermont-Ferrand (France): Université Blaise Pascal.

Portal, A., Fargier, Y., Labazuy, P., Lénat, J.-F., Boivin, P., and Miallier, D. (2019). 3D electrical imaging of the inner structure of a complex lava dome, Puy de Dôme volcano (French Massif Central, France). *J. Volcanol. Geoth. Res.* 373, 97–107. doi:10.1016/j.jvolgeores.2019.01.019

Portal, A., Gailler, L.-S., Labazuy, P., and Lénat, J.-F. (2016). Geophysical imaging of the inner structure of a lava dome and its environment through gravimetry and magnetism. *J. Volcanol. Geoth. Res.* 320, 88–99. doi:10.1016/j.jvolgeores.2016.04.012

Rosas-Carbajal, M., Jourde, K., Marteau, J., Deroussi, S., Komorowski, J.-C., and Gibert, D. (2017). Three-dimensional density structure of La Soufrière de Guadeloupe lava dome from simultaneous muon radiographies and gravity data. *Geophys. Res. Lett.*, 44, 6743–6751. doi:10.1002/2017GL074285

Saracino, G., Ambrosino, F., Bonechi, L., Bross, A., Cimmino, L., Ciaranfi, R., et al. (2017). The MURAVES muon telescope: technology and expected performances. *Ann. Geophys.* 60, 0103. doi:10.4401/ag-7378

Somigliana, C. (1930). *Sul campo gravitazionale esterno del geode ellisoidico*. Atti della Reale Academia Nazionale dei Lincei. Classe di scienze fisiche, matematiche e naturali. Rendiconti, 6, 237–243.

Sparks, R. S. J., and Young, S. R. (2002). The eruption of Soufrière Hills Volcano, Montserrat (1995–1999): overview of scientific results. *Geol. Soc. London Memoirs* 21, 45–69. doi:10.1144/gsl.mem.2002.021.01.03

Tanaka, H., Nagamine, K., Kawamura, N., Nakamura, S. N., Ishida, K., and Shimomura, K. (2001). Development of the cosmic-ray muon detection system for probing internal-structure of a volcano. *Hyperfine Interact.* 138, 521–526. doi:10.1023/A:1020843100008

Tarantola, A. (2005). *Inverse problem theory and methods for model parameter estimation*. Philadelphia, PA: SIAM, 342.

Tenzer, R., Sirguey, P., Rattenbury, M., and Nicolson, J. (2011). A digital rock density map of New Zealand. *Comput. Geosci.* 37, 1181–1191. doi:10.1016/j.cageo.2010.07.010

The CALICE collaboration (2016). First results of the CALICE SDHCAL technological prototype. *J. Instrum.* 11, P04001. doi:10.1088/1748-0221/11/04/P04001

**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.

Copyright © 2021 Barnoud, Cayol, Lelièvre, Portal, Labazuy, Boivin and Gailler. *This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*