# Improvement of density models of geological structures by fusion of gravity data and cosmic muon radiographies

Kevin Jourde, Dominique Gibert, Jacques Marteau

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

## To cite this version:

Kevin Jourde, Dominique Gibert, Jacques Marteau. Improvement of density models of geological structures by fusion of gravity data and cosmic muon radiographies. Geoscientific Instrumentation, Methods and Data Systems, 2015, 4, pp.177-188. 10.5194/gi-4-177-2015. insu-03579527

HAL Id: insu-03579527

https://insu.hal.science/insu-03579527

Submitted on 18 Feb 2022

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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>

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

# Improvement of density models of geological structures by fusion of gravity data and cosmic muon radiographies

K. Jourde¹, D. Gibert²,³, and J. Marteau⁴

¹Institut de Physique du Globe de Paris, Sorbonne Paris Cité, CNRS – UMR7154, Univ. Paris Diderot, Paris, France

²OSUR – Géosciences Rennes, CNRS – UMR6118, Univ. Rennes 1, Rennes, France

³Volcano Observatories, Institut de Physique du Globe de Paris, Paris, France

⁴Institut de Physique Nucléaire de Lyon, CNRS – UMR5822, Univ. Claude Bernard, Lyon, France

*Correspondence to:* D. Gibert (gibert@univ-rennes1.fr) and K. Jourde (jourde@ipgp.fr)

Received: 25 February 2015 – Published in Geosci. Instrum. Method. Data Syst. Discuss.: 2 April 2015

Revised: 6 August 2015 – Accepted: 6 August 2015 – Published: 25 August 2015

**Abstract.** This paper examines how the resolution of small-scale geological density models is improved through the fusion of information provided by gravity measurements and density muon radiographies. Muon radiography aims at determining the density of geological bodies by measuring their screening effect on the natural flux of cosmic muons. Muon radiography essentially works like a medical X-ray scan and integrates density information along elongated narrow conical volumes. Gravity measurements are linked to density by a 3-D integration encompassing the whole studied domain. We establish the mathematical expressions of these integration formulas – called acquisition kernels – and derive the resolving kernels that are spatial filters relating the true unknown density structure to the density distribution actually recovered from the available data. The resolving kernel approach allows one to quantitatively describe the improvement of the resolution of the density models achieved by merging gravity data and muon radiographies. The method developed in this paper may be used to optimally design the geometry of the field measurements to be performed in order to obtain a given spatial resolution pattern of the density model to be constructed. The resolving kernels derived in the joined muon-gravimetry case indicate that gravity data are almost useless for constraining the density structure in regions sampled by more than two muon tomography acquisitions. Interestingly, the resolution in deeper regions not sampled by muon tomography is significantly improved by joining the two techniques. The method is illustrated with examples for the La Soufrière volcano of Guadeloupe.

## 1 Introduction

Determining the density distribution inside geological structures is of major importance in many domains of Earth sciences. The development of space-borne techniques to either determine the geoid shape or directly measure the planetary gravity field dramatically improved the quality of the data available to perform studies at the global and regional scales (e.g. Ménard et al., 2003; Tapley et al., 2005). This significantly boosted new areas of research in hydrology (e.g. Llubes et al., 2004; Longuevergne et al., 2013), erosion (e.g. Mouyen et al., 2012), and climate change (e.g. Chen et al., 2006; Wouters et al., 2014; Song et al., 2015).

For studies performed at local scales, from kilometre down to decametre, classical gravimetry methods remain the main approach used to recover the density distribution underground. Despite huge improvements to gravity meters, either relative or absolute (Lederer, 2009; Riccardi et al., 2001), gravity surveying remains a lengthy, costly and difficult task, especially on volcanoes and rough topography (e.g. Carbone et al., 2003). Meanwhile, the progress and evolution of research domains demand even more challenging capabilities of geophysical imaging of the density distribution inside the Earth. Conventional gravity surveying may quickly reach its limits in applications – hydrology, volcanology, civil engineering, archaeology – where high resolution is mandatory or to monitor density changes underground (e.g. Christiansen et al., 2011; Creutzfeldt et al., 2014). Measurement points generally remain sparse, especially when absolute gravity meters are employed, and are not suitable for producing high-
resolution models of the density structure. This situation is further complicated by the strong non-uniqueness that characterizes the gravimetric inverse problem and by sophisticated signal processing methods needed to isolate the relevant information (e.g. Crossley et al., 2012).

In the domain of natural hazards, determining the volume of potentially unstable rock masses – cliffs, volcanoes, steep landscapes – is of primary importance in identifying the exposed areas and estimating the risk level. In volcanic regions, low-density unconsolidated materials are particularly subject to destabilization due to their low strength and their high fluid content (e.g. Le Friant et al., 2006). If these vulnerable materials are located on active volcanoes, their destabilization may trigger further damage to the edifice due to the rapid decompression of shallow hydrothermal reservoirs. Elaborating models of the density structure of lava domes subject to an intense hydrothermal alteration is of primary importance to better constrain hazard models.

Muon radiography is a new method that allows one to recover the density distribution inside rock volumes at kilometre scale by measuring their screening effect on the natural cosmic muon flux crossing rocks. Since the pioneering works by Nagamine (1995, 2003), Nagamine et al. (1995) and Tanaka et al. (2001), recent studies have illustrated the interest of the method in imaging spatial and temporal variations of the density inside mountains and volcanoes (Tanaka et al. 2005, 2007a, b, 2008, 2009a, b, 2013; Gibert et al., 2010; Cârloganu et al., 2013; Lesparre et al., 2010, 2012c; Shinohara and Tanaka, 2012; Eppelbaum and Khesin, 2012; Portal et al., 2013; Carbone et al., 2014). Muon radiography is a straight-ray transmission method involving a Radon transform that markedly differs from the 3-D integrative gravity method. As will be discussed in the next section, muon radiography uses the flux of muons “coming from above”, and it is limited to the imaging of shallow structures located above the particle detector. Also, muon radiography only provides information on parts of the density structure that are crossed by the rays, contrarily to the gravity method that brings information on the whole density distribution. These differences between both methods motivate the joining of both types of data to elaborate density models of complex geological structures.

Studies combining muon data and gravity measurements remain scarce, and we emphasize the early study by Caffau et al. (1997), who compared muon tomography with gravity measurements. More recently, Davis and Oldenburg (2012) and Nishiyama et al. (2014) presented joined inversions of gravity data and muon tomography using a straightforward linear regularized inversion based on block models. In the present study, we develop a quantitative methodology to examine how information brought by gravimetry data and by muon radiography may be joined to improve the resolution of density models of highly heterogeneous structures like altered active volcanoes. Fusion of information is studied through an approach based on the resolving kernels which

provide a way to quantitatively evaluate the resolution of the resulting density model independently of any particular parametrization (e.g. block discretization). Resolving kernels only depend on the geometrical properties of the data acquisition (i.e. locations of measurement points and telescope acceptance functions), and they allow one to perform prior analysis to evaluate the model improvements that may be expected by joining additional gravimetric data or muon radiographies. In this way, field measurements may be optimized with respect to the characteristics aimed at for the resulting density model. We begin by establishing the relationships between the density structure and both muon tomography and gravity data. Next, we derive the resolving kernels respectively corresponding to individual gravity and muon inversions and to joint gravity–muon inversion. The resolving kernels translate the information contained in the data into information concerning the density structure.

In order to give the reader a practical insight, the theoretical developments of general interest that are derived in the paper are illustrated with examples taken from real field experiments conducted on La Soufrière of Guadeloupe to emphasize the practical interest of combining muon radiographies and gravity measurements. La Soufrière of Guadeloupe is an active volcano located in the Lesser Antilles arc. The last magmatic eruption occurred in AD 1530 when the present lava dome formed, and the last phreatic eruption occurred in 1976. Recent field measurements in its vents (Allard et al., 2014) and sources (Villemant et al., 2014) show a significant regain of activity in the 2006–2012 period. More recently, new vents appeared during our muon tomography experiments. Developing 3-D density models of the lava dome is of primary interest to assess the structure of the edifice and to better constrain the upper hydrothermal system and its related hazards. The muon tomography experiments were already described in various articles (Lespare et al., 2012c; Jourde et al., 2013). Three sites, called Ravine Sud, Rocher Fendu and Savane à Mulets, were explored and are represented in Fig. 1.2. The results showed important density heterogeneities in the volcanic dome well correlated with the surface vent positions. The gravimetry survey is currently running. For the purpose of this article, we simulated 100 measurements regularly spaced on a grid that covers the dome (Fig. 1.1).

## 2 Basic principles of cosmic muon radiography

The small cross section of muons in ordinary matter (Barrett et al., 1952) allows the hard component of the muon spectrum (Dorman, 2004; Tang et al., 2006; Gaisser and Stanev, 2008) to cross hectometres, and even kilometres, of rock. Most muons crossing the rock volume have a negligible scattering relative to the instrument angular resolution and travel along straight trajectories, ranking muon tomography among the class of straight-ray scanning imaging methods
<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 1.** (1) The red dots represent the positions at which we simulated gravimetry measurements on the Soufrière of Guadeloupe; (2) muon tomography data coverage. The lines represent the observation axes of the muon telescope when located at Rocher Fendu (red), Ravine Sud (yellow) and Savane à Mulets (green).</figcaption></figure>

(Malmqvist et al., 1979; Marteau et al., 2011). In practice, muon tomography is performed by using a series of pixelated particle detectors that allow one to determine the trajectories of the muons passing through the rock body. Portable field telescopes presently used sample hundreds of directions and allow one to scan an entire volcano from a single viewpoint in a couple of weeks (Fig. 1.2). By counting the number of muons passing through the target, the attenuation onto the incident muon flux is determined for each sampled direction and is used to produce a radiography of the object opacity (expressed in $g cm^{-2}$) or of average density along ray paths if the object geometry is known.

Lesparre et al. (2010) establish a feasibility formula where the achievable density resolution is related to the measurement duration (i.e. time resolution), the total apparent rock thickness (i.e. total opacity) and the telescope acceptance (i.e. the detection capacity of the matrices). The feasibility formula is written as an inequality and gives practical hints for designing field experiments and evaluating which density heterogeneities can be resolved inside a given geological target, for a given amount of time and a given telescope. In a more recent study, Jourde et al. (2013) present experimental evidence of a flux of upward-going particles that occurs

in certain field conditions. These particles have trajectories parallel to those of the muons emerging from the rock body to radiography, but they travel through the telescope from rear to front. These upward-going particles may constitute a huge Poissonian noise that could strongly alter the radiographies. Jourde et al. (2013) give practical recommendations for choosing experimental sites likely to give the best possible signal-to-noise ratio, and they also put strong constraints on the time resolution of the electronic detection chain necessary to statistically recognize particles coming from the rear face of the telescopes.

### 3 The sampling of the density distribution by muon tomography and gravimetry

Here, we recall the main formula relating the density distribution to the data, i.e. fluxes of muons and gravity measurements. In the inverse problem framework, these formulas describe the forward problem for each method. In the remainder of this paper, we suppose that the muon data have been cleaned of perturbing effects such as upward-going fluxes of particles as described in Jourde et al. (2013).

#### 3.1 Muon tomography

The primary information used in muon tomography consists in cosmic muon flux attenuation measurements resulting from the screening produced by the geological volume to be scanned. Attenuation is measured by counting the number of muons emerging from the volume for each observation axis, $s_m = (r_m, P_m(\varphi, \theta))$, of the telescope (Fig. 1.2). $r_m$ represents the position of the telescope, and $P_m$ the observation axis acceptance pattern which depends on $(\varphi, \theta)$, the azimuth and zenith angles referenced at $r_m$ (see Fig. 2). Note that $r_m$ is the same for all the observation axes on a given site. $P_m$ depends on the telescope geometry and angular orientation on the site. Our standard field telescopes count $31 \times 31$ observation axes and, in a field experiment where the telescope successively occupies several places around the target, the number, $M$, of data may easily reach several hundredths. For example, if we use the muon tomography data from the three Soufrière sites, $M = 3 \times 31 \times 31$. In practice, $M$ is lower, as many axes point downward or above the volcano.

The $P_m$ ($cm^2$) shape depends on the detection matrices' structure (see Lesparre et al. (2012a, b) for further details). It has a steep peak centered on a small solid angle region $\Omega_m$. It is identically null outside $\Omega_m$ (Figs. 2 and 3). Observe that $P_m$ must not be confounded with the integrated pixel acceptance $T_m$ ($cm^2 sr^{-1}$) used for instance in Lesparre et al. (2010),
<figure><img src="image_7.png" /><figcaption>**Figure 2.** Muon tomography reference frame and notations. $R_A$ and $R_T$ respectively are the absolute orthonormal and the instrument reference frames. An observation axis $\mathbf{s}_m = (\mathbf{r}_m, \mathcal{P}_m(\varphi, \theta))$ is represented with $\mathbf{r}_m$ the vector that localizes the telescope position and $\mathcal{P}_m$ its acceptance pattern (restrained to the solid angle $\Omega_m$). The spherical coordinates $(\varphi, \theta)$ here are localizing the steep acceptance peak mentioned in Sect. 3.1.</figcaption></figure>

$$ \mathcal{T}_m = \int_0^{2\pi} \int_0^{\pi} \mathcal{P}_m(\varphi, \theta) \times \sin(\theta) d\theta d\varphi = \int_{\Omega_m} \mathcal{P}_m(\varphi, \theta) d\Omega. \quad (1) $$

The number of muons attributed to a given line of sight actually corresponds to all muons detected in $\Omega_m$. Inside the geological volume the trajectories of these muons describe a conical volume whose apex is located at the telescope, $\mathbf{r}_m$.

The attenuation of the muon flux caused by the rock screen depends on the amount of matter encountered by the particles along their trajectories. For a given straight trajectory $\mathbf{t} = (\mathbf{r}, \varphi, \theta)$ ($\mathbf{r}$ is a telescope site and $(\varphi, \theta)$ the azimuth and zenith angles referenced at $\mathbf{r}$), it is quantified by the density line integral along $\mathbf{t}$ and the opacity

$$ \rho = \int_{t} \rho(\xi) d\xi = L \times \bar{\rho}, \quad (2) $$

where $\rho$ is the density, $L$ the particle path length, and $\bar{\rho}$ is the average density along $\mathbf{t}$. The differential flux associated with $\mathbf{t}$ may be expressed as a function $\delta\phi_t = \frac{\partial^3\phi}{\partial\Omega\partial S}(\varphi, \varphi, \theta)$ ($\mathbf{s}^{-1} \text{cm}^{-2} \text{sr}^{-1}$) that accounts for the muon flux that reaches the instrument. Then the measured flux $\phi_m$ for the $m$th line of sight relates to the opacity via the relationship

$$ \phi_m = \mathcal{T}_m^{-1} \int_{\Omega_m} \mathcal{P}_m(\varphi, \theta) \times \delta\phi_t(\varrho, \varphi, \theta) d\Omega. \quad (3) $$

Note that the integration is restricted to the small solid angle $\Omega_m$ because of the compact support of the observation axis acceptance $\mathcal{P}_m$ (Fig. 3).

<figure><img src="image_8.png" /><figcaption>**Figure 3.** Representation of the acceptance. The horizontal bars labelled $M_1$, $M_2$ and $M_3$ represent the pixelated detection matrices of the telescope. We draw $\mathcal{T}_m$ and $\mathcal{P}_m/\mathcal{T}_m$ for 15 observation axes $\mathbf{s}_m$ symmetrically distributed on the left and right sides of the main axis of the telescope (arbitrarily indexed for $m$ going from 1 to 15; we do not represent all the observation axes, for clarity purposes). The red and black sawtooth-like curves are a merging of the normalized intra-pixel acceptance $\mathcal{P}_m/\mathcal{T}_m$ for even and odd lines of sights respectively. The blue dots represent the 15 discrete values of the integrated pixel acceptance $\mathcal{T}_m$ obtained by integrating each $\mathcal{P}_m$ on the unit sphere (Eq. 1). Observe that the solid angle $\Omega_m$ associated with a given $\mathbf{s}_m$ (here represented for $m=11$) overlaps the solid angle of the neighbouring lines of sight. The dashed lines are plotted along the acceptance steep peaks discussed in Sect. 3.1.</figcaption></figure>

$\delta\phi_t$ is not linearly related to $\rho$; however, for small opacity fluctuations we assume that $\delta\phi_t$ may be approximated by its first-order development around the local average density, $\rho_0(\mathbf{r})$. $\rho_0(\mathbf{r})$ is the prior density model of the geological structure. For a given path $\mathbf{t}$ it reads as

$$ \delta\phi_t(\rho) = \delta\phi_t(\rho_0) + (\rho - \rho_0) \times \frac{d\delta\phi_t}{d\rho}(\rho_0) + o(\rho), \quad (4) $$

where $\rho_0 = \int_t \rho_0(\xi) d\xi$. Rearranging the terms and letting $\alpha_t = \frac{d\delta\phi_t}{d\rho}(\rho_0)$, we obtain

$$ \delta\phi_t(\rho) - \delta\phi_t(\rho_0) \approx \alpha_t \int_t [\rho(\xi) - \rho_0(\xi)] d\xi. \quad (5) $$

Inserting Eq. (5) into Eq. (3), we get the approximate equation

$$ \phi_m - \phi_0 \approx \int_{\Omega_m} d\Omega \int_t \frac{\mathcal{P}_m(\varphi, \theta)}{\mathcal{T}_m} \times \alpha_t \times [\rho(\xi) - \rho_0(\xi)] d\xi, \quad (6) $$

where $\phi_0 = \phi_m(\rho_0)$ is the flux corresponding to the prior density model $\rho_0(\mathbf{r})$ and $\mathbf{t} = (\mathbf{r}_m, \varphi, \theta)$.
In the remainder of the present paper, we shall use the centred and normalized flux,

$$ \tilde{\phi}_m = \frac{\phi_m - \phi_m(\rho_0)}{\phi_m(\bar{\rho}_{\min}) - \phi_m(\bar{\rho}_{\max})} = \frac{\phi_m - \phi_0}{C_{\phi,m}}, \quad (7) $$

where $\bar{\rho}_{\min}$ and $\bar{\rho}_{\max}$ are expected extreme values of the density.

## 3.2 Gravimetry

Gravimetry aims to estimate the gravity field generated by surrounding objects measuring locally the vertical acceleration they produce. The vertical acceleration $g$ is directly related to the density spatial distribution through the Newton law:

$$ g_n = G \int_V \frac{(\mathbf{r}_n - \mathbf{r}) \cdot \mathbf{e}_z}{\| \mathbf{r}_n - \mathbf{r} \|^3} \times \rho(\mathbf{r}) d\mathbf{r}, \quad (8) $$

where the vector $\mathbf{r}_n$ represents the location of the nth measurement point (in our example, $n$ runs from 1 to 100). As for muon tomography, we use the normalized gravity anomaly $\tilde{g}_n$ defined as

$$ \tilde{g}_n = \frac{g_n - g_n(\rho_0)}{|g_n(\bar{\rho}_{\min}) - g_n(\bar{\rho}_{\max})|} = \frac{g_n - g_0}{C_{g,n}}. \quad (9) $$

# 4 Resolving kernel approach

## 4.1 The acquisition kernels

We define $X$, the space that contains the set of continuous $L^2$ functions going from $\mathbb{R}^3$ into $\mathbb{R}$. The 3-D density distribution $\rho$ belongs to $X$ and is related to the muon flux measurements, $\tilde{\phi}_m$, and to the gravity data, $\tilde{g}_n$, through the action of acquisition kernels $\mathcal{G}$ and $\mathcal{M}$ which also belong to $X$. This reads

$$ \tilde{\phi}_m = \langle \mathcal{M}_m, \rho - \rho_0 \rangle_X, \quad m = 1, \dots, M \quad (10) $$

$$ \tilde{g}_n = \langle \mathcal{G}_n, \rho - \rho_0 \rangle_X, \quad n = 1, \dots, N, \quad (11) $$

where $\langle \cdot, \cdot \rangle_X$ is the $X$ inner scalar product, and $M$ and $N$ are respectively the number of muon tomography and gravimetry data. From Eqs. (6)–(9), we obtain explicit expressions for $\mathcal{M}$ and $\mathcal{G}$,

$$ \mathcal{M}_m(\mathbf{r}) = \frac{P_m(\varphi, \theta)}{T_m} \times \frac{\alpha_t}{C_{\phi,m} \xi^2}, \quad (12) $$

$$ \mathcal{G}_n(\mathbf{r}) = \frac{G}{C_{g,n}} \times \frac{(\mathbf{r}_n - \mathbf{r})}{\| \mathbf{r}_n - \mathbf{r} \|^3} \cdot \mathbf{e}_z. \quad (13) $$

Observe that the $1/\xi^2$ term in Eq. (12) comes from the spherical coordinate elementary volume expression, $\xi^2 d\xi d\Omega$, inserted into Eq. (6). Examples of acquisition kernels are plotted in Fig. 4 and discussed in Sects. 5.1 and 5.2. The $X$ structure allows one to introduce prior information

<figure><img src="image_9.png" /><figcaption>(1) gravimetry acquisition kernel, $\mathcal{G}$</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(2) tomography acquisition kernel, $\mathcal{M}$</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 4.** (1) Acquisition kernel of a gravimetric measurement; (2) acquisition kernel for a single observation axis of the muon tomograph. The kernels are normalized with reference to their maximum value and printed in a $\log_{10}$ scale.</figcaption></figure>

into the problem. For instance, the classical inner product of $L^2$ continuous functions,

$$ \langle f, g \rangle_X = \int_V f(\mathbf{r}) \times g(\mathbf{r}) d\mathbf{r}, \quad (14) $$

can be replaced by the weighted inner product,

$$ \langle f, g \rangle_X = \int_V \int_V w(\mathbf{r}', \mathbf{r}'') \times f(\mathbf{r}') \times g(\mathbf{r}'') d\mathbf{r}' d\mathbf{r}'' \quad (15) $$

where the weight function $w(w(\mathbf{r}', \mathbf{r}') > 0, w(\mathbf{r}', \mathbf{r}'') = w(\mathbf{r}'', \mathbf{r}')$ plays the role of a covariance function that may be used to neglect the impact of the free air zone around the studied structure for gravimetry and muon tomography (see Eq. 25 and its comment below). It may also serve to introduce a correlation length for the density variations.

## 4.2 The resolving kernel

The 3-D density distribution, $\hat{\rho}(\mathbf{r})$, obtained by solving the set of linear equations (Eqs. 10, 11), is a degraded version of the true density distribution, $\rho(\mathbf{r})$, both because of the limited number of data available and because of the filtering (i.e. blurring) effect of the acquisition kernels. In the remainder of the paper, we shall use the set of undifferentiated acquisition kernels
$$ \{\zeta_k\} = \{G_n\} \cup \{M_m\} \quad k = 1, \dots, K = M + N, \quad (16) $$

and the set of undifferentiated normalized data

$$ \{d_k\} = \{\tilde{g}_n\} \cup \{\tilde{\phi}_m\} \quad k = 1, \dots, K = M + N. \quad (17) $$

We now formulate the inverse problem in the framework of functional spaces where the family of acquisition kernels constitutes a set of generating functions of a subspace $X_K$ of $X$ of dimension $K$ (Tarantola and Nercessian, 1984; Bertero et al., 1985). This implicitly assumes that the $\zeta_k$ are linearly independent with respect to the retained inner product; i.e. no acquisition kernel can be written as a linear combination of the other kernels. The noticeable instances where this important assumption is not satisfied correspond to situations where several data have been acquired identically, i.e. either at the same location for gravity measurements or with the same position and orientation of the telescope for muon tomography. In such cases the dimension of $X_K$ is reduced since the redundant data may be merged (i.e. averaged) into a single one.

The best density distribution that can be recovered through the inversion process (it is the best because it takes all the information contained in the data and makes the fewest hypotheses about the $X_K$ complementary subspace) is a linear combination of the generating functions,

$$ \hat{\rho}(\mathbf{r}) - \rho_0(\mathbf{r}) = \sum_{k=1}^{K} a_k \times \zeta_k(\mathbf{r}). \quad (18) $$

The components $a_k$ of Eq. (18) are obtained by minimizing the quadratic distance $\epsilon_Y$ between the data and the corresponding values given by the density model

$$ \epsilon_Y = \| \{\zeta_k\}, \rho - \rho_0 \|_X - \{d_k\} \|_Y \\ = \sum_{k=1}^{K} W_k \times (\langle \zeta_k, \rho - \rho_0 \rangle_X - d_k)^2. \quad (19) $$

$Y$ is the weighted Euclidean space that contains the measurements. The weights $W_k$ permit one to introduce prior information about the measurement quality. It is possible to introduce crossed terms $W_{ij}$ if the measurements are not independent, but it is not the case here. We get

$$ a_k = \sum_{j=1}^{K} W_j \times S^{k,j} \times \langle \zeta_j, \rho - \rho_0 \rangle_X, \quad (20) $$

where $S^{k,j}$ is the $(k, j)$ component of the Gram matrix inverse defined as $S_{k,j} = W_j \times \langle \zeta_k, \zeta_j \rangle_X$. Using Eq. (20), Eq. (18) becomes

$$ \hat{\rho}(\mathbf{r}) - \rho_0(\mathbf{r}) = \sum_{k=1}^{K} \zeta_k(\mathbf{r}) \sum_{j=1}^{K} W_j \times S^{k,j} \times \langle \zeta_j, \rho - \rho_0 \rangle_X. \quad (21) $$

The presence of $\langle \zeta_j, \rho \rangle_X$ in the right-hand part of this equation indicates that the density distribution actually recovered, $\hat{\rho}(\mathbf{r})$, is assembled from projections of the true unknown density, $\rho(\mathbf{r})$, onto the acquisition kernels. The recovered density is a filtered version of the true density distribution, and the filter (i.e. the resolving kernel) depends on the data. This can be made more explicit by rewriting Eq. (21) as (Bertero et al., 1985)

$$ \hat{\rho}(\mathbf{r}) - \rho_0(\mathbf{r}) = \int_V \Delta(\mathbf{r}, \mathbf{r}') \times (\rho(\mathbf{r}') - \rho_0(\mathbf{r}')) \mathrm{d}\mathbf{r}', \quad (22) $$

where we introduce the resolving kernel

$$ \Delta(\mathbf{r}, \mathbf{r}') = \sum_{j=1}^{K} b_j(\mathbf{r}) \times \tilde{\zeta}_j(\mathbf{r}'), \quad (23) $$

with

$$ b_j(\mathbf{r}) = \sum_{i=1}^{K} W_i \times \zeta_i(\mathbf{r}) \times S^{i,j}, \quad (24) $$

$$ \tilde{\zeta}_j(\mathbf{r}') = \int_V w(\mathbf{r}', \mathbf{r}'') \times \zeta_j(\mathbf{r}'') \mathrm{d}\mathbf{r}'' \quad (25) $$

$\tilde{\zeta}_j$ is the acquisition kernel $\zeta_j$ modulated by the prior information represented by the $w$ function. For instance, $w$ may be an indicator function used to limit the support of the acquisition kernels to the volume of interest.

## 5 Characterization of the resolving kernels

A resolving kernel, $\Delta(\mathbf{r}, \mathbf{r}')$, is a function defined in the whole space that plays the role of a spatial filter. When applied to the true density distribution, it gives the reconstructed density. The amplitude and the shape of $\Delta$ render the achievable resolution of the reconstructed density structure. According to Eq. (23), it is a linear combination of the acquisition kernels. $\Delta$ may be characterized in different ways by using several properties to quantify its resolution and anisotropy. These properties should be as independent as possible from a specific resolving kernel and allow the user to easily appreciate the resolution and its eventual bias. In the present study, we simply compare the resolving kernels against the ideal kernel represented by a Dirac distribution $\delta(\mathbf{r} - \mathbf{r}')$. This is achieved through the projection $\gamma$ of $\Delta(\mathbf{r}, \mathbf{r}')$ onto a $\delta(\mathbf{r} - \mathbf{r}')$,

$$ \gamma(\mathbf{r}) = \langle \Delta(\mathbf{r}, \mathbf{r}'), \delta(\mathbf{r} - \mathbf{r}') \rangle_X = \int_V w(\mathbf{r}, \mathbf{r}') \times \Delta(\mathbf{r}, \mathbf{r}') \mathrm{d}\mathbf{r}'. \quad (26) $$

Here the $X$ scalar product is defined by Eq. (15). Here, the Hamming function is used as the weight function

$$ w(\mathbf{r}, \mathbf{r}') = \frac{H_L(||\mathbf{r} - \mathbf{r}'||)}{K_w} \times \left[ 1 + \cos\left(\frac{2\pi \times ||\mathbf{r} - \mathbf{r}'||}{L}\right) \right], \quad (27) $$
where $K_w$ is a normalizing constant, $H_L$ is a rectangular pulse that restricts $w$ to $|\mathbf{r}-\mathbf{r}'| \in [-L/2; L/2]$, and $L=25$ m. If $\mathbf{r}$ or $\mathbf{r}'$ is located outside the volcano, we take $w(\mathbf{r}, \mathbf{r}')=0$. This choice is explained in Sect. 5.4.

We now display resolving kernels corresponding to the data acquisition shown in Fig. 1.1 for muon tomography and in Fig. 1.2 for gravimetry. Muon radiographies are taken from three sites equidistantly located along the southern edge of the volcano. Gravity measurements are assumed to be made on a regular grid over the entire lava dome.

Accounting for the fact that the acquisition kernels $\zeta_k$ are either for gravity or for muon tomography (Eq. 16), we successively consider the case of resolving kernels obtained for muon tomography alone, for gravity data alone, and for a combination of muon tomography and gravity data. We compute $\Delta(\mathbf{r}, \mathbf{r}')$ for two positions $\mathbf{r} = \{\mathbf{r}_1; \mathbf{r}_2\}$ located along a vertical line that goes through the centre of the dome (Fig. 5). Points $\mathbf{r}_1$ and $\mathbf{r}_2$ are respectively inside and below the volume of the lava dome spanned by the lines of sight of the telescopes (Fig. 1.1). The parameter $\gamma$ is computed and plotted on the four characteristic slices represented in Fig. 5.

## 5.1 Gravimetry kernels

Figure 4.1 shows a gravimetry acquisition kernel $G$. Remember that the data used in the present study are normalized relative to a reference model with density $\rho_0$ (Eq. 9). The gravimetry acquisition kernels are very sensitive to density fluctuations close to the measurement point because of its $1/r^2$ term. The gravity data are actually the component of the gravity field anomaly taken along the local vertical, and the acquisition kernel becomes less and less sensitive as we get closer to the horizontal plane that contains the measurement point.

The gravimetry inverse problem is systematically ill-posed (e.g. Al-Chalabi, 1971) because no matter the number of measurements, the resolving kernel mostly integrates information around the measurement positions, i.e. near the surface. An illustration of this problem is given for the resolving kernels of $\mathbf{r}_1$ and $\mathbf{r}_2$ (Figs. 6.1 and 7.1). For gravimetry inversions it is more realistic to model $\rho(\mathbf{r})$ by a function that depends on a few discrete parameters (even if it means losing the linearity between the data and the measurements) rather than trying a continuous inversion.

Observe that $G$ integrates the density over the entire volume and provides information for point $\mathbf{r}_2$ located below the lines of sight of the telescope.

## 5.2 Muon tomography kernels

Figure 4.2 shows a typical muon tomography acquisition kernel $M$ (Eq. 12). It has a conical shape whose aperture angle depends on the distance between the front and rear detection matrices of the telescope. The apex of the kernel is located at the telescope, and the kernel widens as we move away from

<figure><img src="image_12.png" /><figcaption>(1) vertical observation planes</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(2) horizontal observation planes</figcaption></figure>

Figure 5. 3-D views of the cross sections used to represent the resolving kernels in Figs. 6 and 7. The resolving kernels are computed at points 1 and 2. Point 1 is located at a level Z1 in the part of the lava dome scanned by the lines of sight of the muon telescope (Fig. 1.2) and point 2 is located at Z2 below the ray coverage of the telescope.

the telescope; thus, the local sensitivity is decreasing. Moreover, the triangular shape of the intra-pixel acceptance $P_m$ (see Figs. 2 and 3) makes the sensitivity maximum along the main line of sight $s_m$.

Figure 4.2 shows we are as sensitive to a density change occurring on a few tenths of metres in front of the telescope as to the same change happening on a few hundred metres beside the volcano. It reveals how deterministic the telescope position is. If one desires to image or monitor a specific region belonging to a bigger structure, the measurement will be much more sensitive if the telescope is in front of it. The important heterogeneities inside the muon tomography acquisition kernels forbid us from using the Radon transform mathematical corpus. For an equivalent resolution and scanning, the kernels can be regularized by taking the telescope away from the volcano and reducing the angular aperture. We then get into the typical experimental conditions of a medical X-ray tomography. But, the consequences are a weaker particle flux (a longer acquisition time) and a greater sensitivity to potential noises. So, a compromise has to be found, but
<figure><img src="image_14.png" /><figcaption>(1) Gravimetry data</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(2) Ravine Sud tomography data</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(3) All sites tomography data</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>(4) Gravimetry and tomography data</figcaption></figure>

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

**Figure 6.** Resolving kernel at point $\mathbf{r}_1$ (Fig. 5): (1) the gravity data alone; (2) Ravine Sud muon radiography alone; (3) combination of three muon radiographies; and (4) joined muon and gravity data sets. See Fig. 1 for the locations of gravity measurements and the three sites for muon radiographies. The resolving kernel's absolute value is normalized with reference to the value computed at point $\mathbf{r}_1$ and represented with a $\log_{10}$ scale.

<figure><img src="image_19.png" /><figcaption>(1) Gravimetry data</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>(2) Gravimetry and tomography data</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>(3) Gravimetry and tomography data + prior information</figcaption></figure>

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

**Figure 7.** Resolving kernel at point $\mathbf{r}_2$ (Fig. 5): (1) gravity data alone; (2) joined muon radiographies and gravity. (3) is obtained by joining muon radiographies, gravity, and some prior information about the density spatial correlation. See Fig. 1 for the locations of gravity measurements and the three sites for muon radiographies. The resolving kernel's absolute value is normalized with reference to the value computed at point $\mathbf{r}_2$ and represented with a $\log_{10}$ scale.

the actual lack of understanding of the noises and the already very long acquisition times we are facing lead us to take the telescope the closest we can to the volcano.

We draw the reader's attention to the fact that, despite their compact support, the acquisition kernels overlap each other for neighbouring main lines of sight (see Fig. 3). As will be seen below, this characteristic is fundamental to understanding the shape of the resolving kernels.

A muon tomography resolving kernel is a linear combination of muon tomography acquisition kernels $M$ (Eq. 23), and the inversion process optimizes the $b_i$ coefficients to obtain the best density model (Eq. 24). Figure 6.2 and 6.3 show
the resolving kernel for point **r**₁ for different combinations of tomography data sets.

When using data acquired from a single place located at the southernmost edge of the volcano (Ravine Sud; see Fig. 1.2), the resolving kernel (Fig. 6.2) encompasses lines of sight spanning a limited range of azimuths. Consequently, the filtering effect of Δ(**r**₁, **r**') integrates ρ along a long narrow cone to give the estimated density ρ̂(**r**₁). The fact that Δ(**r**₁, **r**') does not have a compact support like the M kernels comes from the overlapping of neighbouring acquisition kernels that produces a transfer of information among lines of sight.

When simultaneously using all three muon tomography sites, the resolving kernel includes acquisition kernels that span a wider range of sight azimuths. Consequently, the resolving kernel is more localized onto point **r**₁ (Fig. 6.3). However, the number of radiographies remains small, and the kernel has a spider shape visible in the horizontal slice of Fig. 6.3.

Observe that the resolving kernel Δ(**r**₂, **r**) equals 0, since all acquisition kernels are null in this part of the volcano.

## 5.3 Joined muon tomography and gravimetry kernels

We now consider resolving kernels computed by using both muon *M* and gravity *G* acquisition kernels.

Gravimetry does not improve significantly the inversion process at **r**₁, and the resolving kernel Δ(**r**₁, **r**') (Fig. 6.4) looks very similar to the one obtained for the muon radiographies alone (Fig. 6.3). The information provided by muon tomography is dominant relative to gravimetry, except at the immediate vicinity of the gravity measurement points.

The situation is very different for point **r**₂, where the resolving kernels obtained by joining muon radiographies and gravity data (Fig. 7.2) appear very different from the gravity kernel (Fig. 7.1). The most conspicuous effect is that muon data compensates for the great sensitivity of gravimetry at near-surface locations by shifting the centre of the mass of the resolving kernel downward. This considerably improves the vertical resolution achievable in the deepest parts of the volcano.

The conclusions are different if only one tomography acquisition is available. In that case the gravimetry measurements have an impact on the upper part of the dome because they contribute to resolving the ambiguity about the anomaly spatial depth relative to the acquisition position. But, then the zone below the dome will lack data to be properly constrained.

## 5.4 Impact of prior information

The choice made for the *X* and *Y* weight functions w(**r', **r'') (Eq. 15) and Wᵢ=₁...K (Eq. (19)) has an important influence on the obtained resolving kernels Δ(**r**, **r'**).

In the *X* space, the diagonal term w(**r', **r'**) = w(**r', **r'** = **r'**) permits one to adjust the local degree of prior knowledge of ρ(**r**). For instance, w(**r**, **r**) = 0 in regions where ρ(**r**) is assumed to be sufficiently well known to have no impact on our measurements. This corresponds to situations where ρ₀(**r**) = ρ(**r**) and where the concerned regions do not have to be accounted for in the inversion process. In our case, we use it to cancel the free-air impact on muon tomography and gravimetry, but we can also constrain it to incorporate direct field measurements of the density.

The non-diagonal part w(**r', **r'** ≠ **r'**) may be used to introduce a spatial correlation in ρ. This can be done through ζ, the convolution of ζ with w (Eq. 25). Here, we use a simple Hamming function with a 25 m correlation length everywhere in the dome (Eq. 27), and ζ is a smoothed version of ζ which attenuates the 1/r² effect previously mentioned (for muon tomography, it permits one to get closer to the X-ray tomography experimental conditions previously detailed). The correlation introduced by the Hamming function increases the acquisition kernel sensitivity further from the measurement point toward the central and northern parts of the dome. This produces a better localization of the resolving kernel at **r**₂, as can be checked by comparing Fig. 7.3 with Fig. 7.2, where no spatial correlation was applied. The counterpart of this effect is a de-sharpening of the kernel at point **r**₂. w is a regularizing low-pass filter that removes spurious short-wavelength fluctuations in the density model and reduces the ill-conditioning of the inverse problem (e.g. Bertiero et al., 1988).

The choice of w is problem-dependent and must be sustained by prior knowledge. The Hamming function acts as a low-pass filter with a limited support compatible with the large homogeneous zones observed in the field: massiveandesite, hydrothermally altered material and possibly large cavities.

In the *Y* space, the weights Wᵢ=₁...K allow one to assign different quality factors to the available data at one inversion location. For instance, in muon tomography, the *Ws* permit one to account for the fact that not all observation axes have the same integrated acceptance Tm (Fig. 3). The quality of the gravity data strongly depends on the ground stability (i.e. tilt stability during measurement sequences) and the presence of wind (i.e. vibrations of the gravity meter). The non-diagonal terms Wᵢⱼ,ᵢⱼ are null, as the measurements are independent of the ones from the others.

## 5.5 γ maps

The γ(**r**) index defined in Eq. (26) may be used to estimate the resolution achievable everywhere in the volcano. Figure 8 shows slices of the γ function obtained for the gravity data (Fig. 8.1) and by joining the three muon tomography data sets together with the gravity measurements (Fig. 8.2). The gravimetry γ slices clearly reveal the important sensitivity of the data to density variations located in the immediate vicin-
<figure><img src="image_23.png" /><figcaption>**Figure 8.** Representation of $|\gamma|$ on the four slices defined in Fig. 5. The results are represented with a $\log_{10}$ scale for (1) only gravimetry and (2) gravimetry and tomography data.</figcaption></figure>

ity of the measurement points and the very low sensitivity to the density structure located deeper in the lava dome.

The muon-gravimetry $\gamma$ slices confirm the results obtained for points $\{r_1; r_2\}$ and show the considerable improvement of the resolution obtained when jointly using the muon and gravity data sets. They also reproduce the asymmetric resolution due to the conical shape of the muon acquisition kernels $M$. Since the places occupied by the telescope are located along the southern edge of the volcano, a finer resolution is obtained for the southern part of the lava dome. This corresponds to the dark-red circular sector visible in the left horizontal slice in Fig. 8.2. As expected, the resolution is coarser in the central and northern parts of the dome. The same slice shows that the north-eastern region of the volcano is resolved by the gravity data alone since no lines of sight of the telescope cross this part of the dome.

The $\gamma$ map is a useful tool to plan an acquisition survey. One can easily compute how $\gamma$ is changed with different possible measurement campaigns and select the most pertinent one depending on the region of interest, the available time in the field and the accessibility of the site. This choice is critical as muon tomography acquisitions are long (a few weeks) and gravity measurements delicate. It can also be used to design a mesh for the problem. The meshing element density can roughly follow the $\gamma$ map fluctuations.

We emphasize that other definitions may be used for $\gamma$ and that a single index may prove insufficient to characterize the shape of resolving kernels. Consequently, we recommend performing a 3-D examination of individual resolving kernels at locations of particular importance (i.e. like detecting places where density changes occur).

## 6 Conclusions

The resolving kernel analysis discussed in the present study allows one to quantitatively assess the way in which gravity data and density muon radiographies may be joined to improve the spatial resolution of density models of geological structures.

Thanks to the compact support of muon acquisition kernels, high resolution is achievable in parts of the density model sampled by muon radiography. The resolution actually obtained depends on the number and geometrical arrangement of the radiographies available for the model construction (compare Fig. 6.2 and 6.3 for the La Soufrière example). In parts of the model sampled by muon radiography, the fusion of gravity data and muon radiographies does not lead to significant improvement of the density model (compare Fig. 6.3 and 6.4 for the La Soufrière example).

A main result concerns the improvement of the resolution obtained in the deeper parts of the density model when joining muon and gravity data, and despite the fact that these parts are not directly sampled by muon tomography. Actually, a fraction of information brought by the muon data is transferred to the deep regions of the model through the long-range coupling of the gravity acquisition kernels used to construct the joined resolving kernels (compare Fig. 7.1 and 7.2 for the La Soufrière example).

The muon tomography acquisition kernel has a conic shape, and noise considerations lead us to put the cone apex just in front of the studied structure (see Fig. 4.2 for the La Soufrière example). It results in a large decrease in the sensitivity between the front and rear of the volcano. This problem adds to the heterogeneous tomography sampling and forbids us from using the standard Radon transform usually adopted in X-ray tomography medical experiments to inverse the density. It also shows how deterministic the telescope position is if one desires to image or monitor a specific region belonging to a bigger body.
The positive weight function w of the inner product (Eq. 15) may be used to introduce prior information both by limiting the support of the density distribution to reconstruct and by introducing a spatial correlation smoothing the $1/r^2$ effect of muon tomography and gravimetry acquisition kernels. This extends the range of sensitivity of the measurements and results in a solution with a more homogeneous quality (compare Fig. 7.2 and 7.3 for the La Soufrière example). The γ maps give an overview of the resolving kernel geometry in the density model and may be used to optimally plan future acquisition surveys to improve the resolution in parts of the model (see Fig. 8 for the La Soufrière example).

*Acknowledgements.* This study is part of diaphane project anr-14-ce04-0001. We acknowledge the financial support from the UnivEarthS Labex program of Sorbonne Paris Cité (anr-10-labx-0023 and anr-11-idex-0005-02). This study also received funding from Swisstopo through the Mont-Terri MD project. This is IPGP contribution 3623.

Edited by: L. Eppelbaum

## References

Al-Chalabi, M.: Some studies relating to non-uniqueness in gravity and magnetic inverse problems, Geophysics, 36, 835–855, 1971.

Allard, P., Aiuppa, A., Beauducel, F., Calabrese, S., Di Napoli, R., Crispi, O., Gaudin, D., Parello, F., Hammouya, G., and Tamburrello, G.: Steam and gas emission rate from La Soufrière volcano, Guadeloupe (Lesser Antilles): Implications for the magnetic supply during degassing unrest, Chem. Geol., 348, 76–93, doi:10.1016/j.chemgeo.2014.06.019, 2014.

Barrett, P. H., Bollinger, L. M., Cocconi, G., Eisenberg, Y., and Greisen, K.: Interpretation of Cosmic-Ray Measurements Far Underground, Rev. Modern Phys., 24, 133–178, 1952

Bertero, M., De Mol, C., Pike, E. R., Bertero, M., De Mol, C., and Pike, E. R.: Linear inverse problems with discrete data, I. General formulation and singular system analysis, Inv. Probl., 1, 301, 1985.

Bertero, M., De Mol, C., Pike, E. R., Bertero, M., De Mol, C., and Pike, E. R.: Linear inverse problems with discrete data: II. Stability and regularisation, Inv. Probl., 4, 573, 1988.

Caffau, E., Coren, F., and Giannini, G.: Underground cosmic-ray measurement for morphological reconstruction of the Grotta Gigante natural cave, Nucl. Instrum. Meth. A, 385, 480–488, 1997.

Carbone, D., Budetta, G., Greco, F., and Rymer, H.: Combined discrete and continuous gravity observations at Mount Etna, J. Volcanol. Geoth. Res., 123, 123–135, 2003.

Carbone, D., Gibert, D., Marteau, J., Diament, M., Zuccarello, L., and Galichet, E.: An experiment of muon radiography at Mt Etna (Italy), Geophys. J. Int., 196, 633–643, 2014.

Cârloganu, C., Niess, V., Béné, S., Busato, E., Dupieux, P., Fehr, F., Gay, P., Miallier, D., Vulpescu, B., Boivin, P., Combaret, C., Labazuy, P., Laktineh, I., Lénat, J.-F., Mirabito, L., and Portal, A.: Towards a muon radiography of the Puy de Dôme, Geosci. Instrum. Method. Data Syst., 2, 55–60, doi:10.5194/gi-2-55-2013, 2013.

Chen, J. L., Wilson, C. R., and Tapley, B. D.: Satellite gravity measurements confirm accelerated melting of Greenland ice sheet, Science, 313, 1958–1960, 2006.

Christiansen, L., Lund, S., Andersen, O. B., Binning, P. J., Rosbjerg, D., and Bauer-Gottwein, P.: Measuring gravity change caused by water storage variations: Performance assessment under controlled conditions, J. Hydrol., 402, 60–70, 2011.

Creutzfeldt, B., Troch, P. A., Güntner, A., Ferré, T., Graeff, T., and Merz, B.: Storage-discharge relationships at different catchment scales based on local high-precision gravimetry, Hydrol. Process., 28, 1465–1475, 2014.

Crossley, D., De Linage, C., Hinderer, J., Boy, J. P., and Famiglietti, J.: A comparison of the gravity field over Central Europe from superconducting gravimeters, GRACE and global hydrological models, using EOF analysis, Geophys. J. Int., 189, 877–897, 2012.

Davis, K. and Oldenburg, D. W.: Joint 3D of muon tomography and gravity data to recover density, ASEG Extend. Abstr., 1, 1–4, 2012.

Dorman, L. I.: Cosmic Rays in the Earth’s Atmosphere and Underground, Kluwer Academic Publishers, Dordrecht, Boston, London, 2004.

Eppelbaum, L. V. and Khesin, B. E.: Geophysical Studies in the Caucasus, Springer, Berlin, Heidelberg, p. 411, 2012.

Gaisser, T. and Stanev, T.: Cosmic Rays, Phys. Lett. B, 667, 254–260, doi:10.1016/j.physletb.2008.07.028, 2008.

Gibert, D., Beauducel, F., Déclais, Y., Lesparre, N., Marteau, J., Nicollin, F., and Tarantola, A.: Muon tomography: Plans for observations in the Lesser Antilles, Earth Planets Space, 62, 153–165, 2010.

Jourde, K., Gibert, D., Marteau, J., de Bremond d'Ars, J., Gardien, S., Girerd, C., Ianigro, J.-C., and Carbone, D.: Experimental detection of upward-going cosmic particles and consequences for correction of density radiographies of volcanoes, Geophys. Res. Lett., 40–24, 6334–6339, 2013.

Lederer, M.: Accuracy of the relative gravity measurement, Acta Geodyn. Geomater., 6, 383–390, 2009.

Le Friant, A., Boudon, G., Komorowski, J. C., Heinrich, P., and Semet, M. P.: Potential flank-collapse of Soufrière volcano, Guadeloupe, Lesser Antilles? Numerical simulation and hazards, Nat. Hazards, 39, 381–393, 2006.

Lesparre, N., Gibert, D., Marteau, J., Déclais, Y., Carbone, D., and Galichet, E.: Geophysical muon imaging: feasibility and limits, Geophys. J. Int., 183, 1348–1361, 2010.

Lesparre, N., Gibert, D., and Marteau, J.: Bayesian Dual Inversion of Experimental Telescope Acceptance and Integrated Flux for Geophysical Muon Tomography, eophys. J. Int., 188, 490–497, doi:10.1111/j.1365-246X.2011.05268.x, 2012a.

Lesparre, N., Marteau, J., Déclais, Y., Gibert, D., Carlus, B., Nicollin, F., and Kergosien, B.: Design and operation of a field telescope for cosmic ray geophysical tomography, Geosci. Instrum. Method. Data Syst., 1, 33–42, doi:10.5194/gi-1-33-2012, 2012b.

Lesparre, N., Gibert, D., Marteau, J., Komorowski, J.-C., Nicollin, F., and Coutant, O.: Density Muon Radiography of La Soufrière of Guadeloupe Volcano: Comparison with Geological, Electrical Resistivity and Gravity data, Geophys. J. Int., 190, 1008–1019, 2012c.
Llubes, M., Florsch, N., Hinderer, J., Longuevergne, L., and Amalvict, M.: Local hydrology, the Global Geodynamics Project and CHAMP/GRACE perspective: some case studies, J. Geodynamics, 38, 355–374, 2004.

Longuevergne, L., Wilson, C. R., Scanlon, B. R., and Crétaux, J. F.: GRACE water storage estimates for the Middle East and other regions with significant reservoir and lake storage, Hydrology. Earth Syst. Sci., 17, 4817–4830, doi:10.5194/hess-17-4817-2013, 2013.

Malmqvist, L., Jönsson, G., Kristiansson, K., and Jacobsson, L.: Theoretical studies of in-situ rock density determinations using underground cosmic-ray muon intensity measurements with application in mining geophysics, Geophysics, 44, 1549–1569, 1979.

Marteau, J., Gibert, D., Lesparre, N., Nicollin, F., Noli, P., and Giacoppo, F.: Muons tomography applied to geo-sciences and volcanology, Nucl. Instrum. Meth. A, 695, 23–28, doi:10.1016/j.nima.2011.11.061, 2011.

Ménard, Y., Fu, L.-L., Escudier, P., Parisot, F., Perbos, J., Vincent, P., Desai, S., Haines, B., and Kunstmann, G.: The Jason-1 mission special issue: Jason-1 calibration/validation, Mar. Geodesy, 26, 131–146, 2003.

Mouyen, M., Masson, F., Hwang, C., Cheng, C. C., Le Moigne, N., Kao, R., and Hsieh, W. C.: Erosion effects assessed by repeated gravity measurements in southern Taiwan, Geophys. J. Int., 198, 187–197, 2012.

Nagamine, K.: Geo-tomographic observation of inner-structure of volcano with cosmic-ray muons, J. Geogr., 104, 998–1007, 1995.

Nagamine, K.: Introductory Muon Science, Cambridge University Press, Cambridge, UK, 208 pp., 2003.

Nagamine, K., Iwasaki, M., Shimomura, K., and Ishida, K.: Method of probing inner-structure of geophysical substance with the horizontal cosmic-ray muons and possible application to volcanic eruption prediction, Nucl. Instrum. Meth. A, 356, 585–595, 1995.

Nishiyama, R., Tanaka, Y., Okubo, S., Oshima, H., Tanaka, H., and Maekawa, T.: Integrated processing of muon radiography and gravity anomaly data toward the realization of high-resolution 3D density structural analysis of volcanoes: Case study of Showa-Shinzan lava dome, Usu, Japan, J. Geophys. Res.-Solid Ea., 119, 699–710, 2014.

Portal, A., Labazuy, P., Lénat, J.-F., Béné, S., Boivin, P., Busato, E., Cárloganu, C., Combaret, C., Dupieux, P., Fehr, F., Gay, P., Laktineh, I., Miallier, D., Mirabito, L., Niess, V., and Vulpescu, B.: Inner structure of the Puy de Dôme volcano: cross-comparison of geophysical models (ERT, gravimetry, muon imaging), Geosci. Instrum. Method. Data Syst., 2, 47–54, doi:10.5194/gi-2-47-2013, 2013.

Riccardi, U., Rosat, S., and Hinderer, J.: Comparison of the Microg LaCoste gPhone-054 spring gravimeter and the GWR-C026 superconducting gravimeter in Strasbourg (France) using a 300-day time series, Metrologia, 48, 28, 2011.

Shinohara, H. and Tanaka, H.: Conduit magma convection of a rhyolitic magma: Constraints from cosmic-ray muon radiography of Iwodake, Satsuma-Iwojima volcano, Japan, Earth Planet. Sc. Lett., 349–350, 87–97, 2012.

Song, C., Ke, L., Huang, B., and Keith, S. R.: Can mountain glacier melting explains the GRACE-observed mass loss in the south-

east Tibetan Plateau: from a climate perspective?, Global Planet. Change, 124, 1–9, 2015.

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

Tanaka, H., Nagamine, K., Nakamura, S. N., and Ishida, K.: Radiographic measurements of the internal structure of Mt. West Iwate with near horizontal cosmic ray muons and future developments, Nucl. Instrum. Meth. A, 555, 164–172, 2006.

Tanaka, H., Nakano, T., Takahashi, S., Yoshida, J., Ohshima, H., Maekawa, T., Watanabe, H., and Niwa, K.: Imaging the conduit size of the dome with cosmic ray muons: the structure beneath Showa Shinzan Lava Dome, Japan, Geophys. Res. Lett., 34, L22311, doi:10.1029/2007GL031389, 2007a.

Tanaka, H., Nakano, T., Takahashi, S., Yoshida, J., Takeo, M., Oikawa, J., Ohminato, T., Aoki, Y., Koyama, E., Tsuji, H., and Niwa, K.: High resolution imaging in the inhomogeneous crust with cosmic-ray muon radiography: The density structure below the volcanic crater floor of Mt. Asama, Japan, Earth Planet. Sc. Lett., 263, 104–113, 2007b.

Tanaka, H., Nakano, T., Takahashi, S., Yoshida, J., Takeo, M., Oikawa, J., Ohminato, T., Aoki, Y., Koyama, E., Tsuji, H., Ohshima, H., Maekawa, T., Watanabe, H., and Niwa, K.: Radiographic imaging below a volcanic crater floor with cosmic-ray muons, Am. J. Sci., 308, 843–850, 2008.

Tanaka, H., Uchida, T., Tanaka, M., Shinohara, H., and Taira, H.: Cosmic-ray muon imaging of magma in a conduit: De-gassing process of Satsuma-Iwojima Volcano, Japan, Geophys. Res. Lett., 36, L01304, doi:10.1029/2008GL036451, 2009a.

Tanaka, H., Uchida, T., Tanaka, M., Takeo, M., Oikawa, J., Ohminato, T., Aoki, Y., Koyama, E., and Tsuji, H.: Detecting a mass change inside a volcano by cosmic-ray muon radiography (muography): First results from measurements at Asama volcano, Japan, Geophys. Res. Lett., 36, L17302, doi:10.1029/2009GL039448, 2009b.

Tanaka, H. K.: Subsurface density mapping of the earth with cosmic ray muons, Nucl. Phys. B-, 243, 239–248, 2013.

Tang, A., Horton-Smith, G., Kudryavtsev, V. A., and Tonazzo, A.: Muon simulations for Super-Kamiokande, KamLAND, and CHOOZ, Phys. Rev. D, 74, 053007, doi:10.1103/PhysRevD.74.053007, 2006.

Tapley, B., Ries, J., Bettadpur, S., Chambers, D., Cheng, M., Condi, F., Gunter, B., Kang, Z., Nagel, P., Pastor, R., Pekker, T., Poole, S., and Wang, F.: GGM02 – An improved Earth gravity field model from GRACE, J. Geodesy, 79, 467–478, 2005.

Tarantola, A. and Nercessian, A.: Three-dimensional inversion without blocks, Geophys. J. Int., 76, 299–306, 1984.

Villemant, B., Komorowski, J. C., Dessert, C., Michel, A., Crispi, O., Hammouya, G., Beauducel, F., and De Chabalier, J. B.: Evidence for a new shallow magma intrusion at La Soufrière of Guadeloupe (Lesser Antilles), insights from longterm geo-chemical monitoring of halogen-rich hydrothermal fluids, J. Volcanol. Geoth. Res., 285, 247–277, doi:10.1016/j.jvolgeores.2014.08.002, 2014.

Wouters, B., Bonin, J. A., Chambers, D. P., Riva, R. E. M., Sassen, I., and Wahr, J.: GRACE, time-varying gravity, Earth system dynamics and climate change, Rep. Prog. Phys., 77, 116801, doi:10.1088/0034-4885/77/11/116801, 2014.