Robust Bayesian Joint Inversion of Gravimetric and Muographic Data for the Density Imaging of the Puy de Dôme Volcano (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.


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 3 dense structure with a precision of 5% on the density, an angular resolution of 3°by 3°and a 1 m 2 detector requires an exposure of about 100 days of effective acquisition time 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(Nishiyama et al., , 2016Gó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  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 leastsquares 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.

JOINT INVERSION METHOD
The medium is described with a three-dimensional mesh of equally spaced nodes of densities ρ. 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 γ 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 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 C 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: 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 C ρ is constructed with a Gaussian spatial covariance function: where σ ρ 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, C μ σ 2 μ I and the offset formulation (Eq. 9) simplifies to : 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): 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 densitỹ which is the center of the a posteriori Gaussian density distribution with the associated a posteriori covariance matrix (Tarantola, 2005) The diagonal elements ofC ρ are the a posteriori variances of the estimated densitiesρ for each node. The modeled data computed from the estimated densityρ are given by: 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: 3 APPLICATION TO THE PUY DE DÔME VOLCANO

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 3 (farfield). 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 3 (mid-field). For masses inside the volume of interest, we estimate the optimal density of correction of 1,890 kg/m 3 (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 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 2 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  Frontiers in Earth Science | www.frontiersin.org January 2021 | Volume 8 | Article 575842 5 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°a bove 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 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 3 , which is 325 kg/m 3 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 3 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 3 for most of the bins. The bias is therefore underestimated compared to the 325 kg/m 3 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 3 to 2,080 kg/m 3 with an average of 1,839 kg/m 3 (Figure 4), in accordance to the mean density estimated from the gravimetric data. When water-saturated, the average density reaches 2,044 kg/m 3 with densities up to 2,200 kg/m 3 . 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).

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 σ ρ 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 (σ ρ , λ) (50 kg/m 3 , 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 illustrating the specificities of the method. The optimal muographic inversion result is obtained for (σ ρ , λ) (2, 100 kg/m 3 , 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 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 3 . (A) Independent inversion of the gravimetric data. The density color scales are centered on 1,890 kg/m 3 , 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 3 , 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 3 , 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.
Frontiers in Earth Science | www.frontiersin.org January 2021 | Volume 8 | Article 575842 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 (σ ρ , λ) (50 kg/m 3 , 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 Lelièvre et al., 2019). The addition of view points greatly helps constraining the model close to the surface , in this way improving the imaging at depth .
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 3 . We recall that the residual Bouguer anomalies have been obtained by substracting a constant density of 1,890 kg/m 3 , while the muographic data correspond to absolute densities. Consequently, the offset c corresponds to a mean difference of 240 kg/m 3 between the muographic and the gravimetric data. This is lower than the difference of 325 kg/m 3 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 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 3 . Iso densities of 1,690 and 2,090 kg/m 3 , corresponding to the averaged density ±200 kg/m 3 , 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.
Frontiers in Earth Science | www.frontiersin.org January 2021 | Volume 8 | Article 575842 8 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 3 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 3 ) than to the distribution of the water-saturated rock sample densities (averaged of 2,044 kg/m 3 ). Even if the joint inversion result was centered on the average density from the muographic data (1,565 kg/m 3 ) or on an intermediate density value, the whole densities would decrease by a few hundreds of kg/m 3 , 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 standarddeviations ( Figure 6C), the 3D density structure of the Puy de Dôme is shown using isosurfaces of 1,690 kg/m 3 and 2,090 kg/m 3 , corresponding to the averaged density of 1,890 kg/m 3 ± 200 kg/m 3 (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.

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 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 Frontiers in Earth Science | www.frontiersin.org January 2021 | Volume 8 | Article 575842 9 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.