1 Introduction

Cosmic ray muons are naturally produced in the atmosphere by the interactions of primary cosmic rays with the nuclei present in the atmosphere itself. They reach the surface of the Earth at a rate of about 1/cm2/min with a well-known, broad, continuous energy distribution and a mean energy of ~4 GeV (Olive et al. 2014). The attenuation of cosmic ray muons passing through matter can be used to estimate the density of objects through which they pass. Compared to background muon flux rates, the fraction of muons detected along a given path provides the integrated density length along that path.

The first three-dimensional density image of a geophysical object using cosmic ray muons (Tanaka et al. 2007, 2010) imaged a volcano using low-angle muon data from a detector located at multiple positions around it. Later, the same team obtained a tomographic image of a volcano combining muon attenuation data from two different angles with gravity measurements (Nishiyama et al. 2014), using a topographic map of the volcano as a prior. Some applications with deployment of detectors in tunnels have also been conducted in ore exploration (Bryman et al. 2014) and geologic mapping (Tanaka 2015). We present muon tomography of the Los Alamos mesa obtained using a detector placed at multiple positions within a tunnel beneath the target overburden. Our 3D inversion recovers well the overlying air/rock interface and the bulk density of the overburden, without invoking any prior constraint on the topography of the mesa.

2 Experiment Site

Our data were acquired at four locations within a tunnel, excavated horizontally under the main Los Alamos town-site mesa (Rowe et al. 2015) (Fig. 1). The tunnel has a length of about 100 m and an overburden ranging between 0 m (at the entrance) and 100 m at its innermost point. The long axis of the tunnel is oriented N5°W. Now decommissioned and de-classified from its Cold War status, the tunnel was used in this study to generate a tomographic density image of the mesa by inverting the measured attenuation of cosmic ray muons through the overburden.

Fig. 1
figure 1

Field experiment site and layout. a Map showing position within New Mexico (inset) of the field area and a topographic map of the Los Alamos Canyon site with tunnel entrance marked by a red star. Contour intervals are in feet above sea level. b Plan view of the tunnel (after McGehee et al. 2004)

Full size image

The mesa is comprised of the Quaternary Bandelier Tuff, a sequence of ash-flow tuffs deposited during the most recent major caldera-forming eruptions of the nearby Jemez volcano, located approximately 10 km west of the site. Figure 1a is a location map (inset) for the target with detailed topography of Los Alamos Canyon and the mesa above it. A red star indicates the tunnel entrance. Figure 1b illustrates the configuration of the tunnel and ancillary structures.

At the tunnel site (Fig. 2) the Tshirege Member (upper unit) of the Bandelier Tuff is exposed (Lavine et al. 2003; Lewis et al. 2009). This member can be further subdivided into cooling units, which have demonstrable mineralogical and physical variations due to the episodic nature of the Tshirege Member eruptive sequence.

Fig. 2
figure 2

Field site exposure and stratigraphic column. a Tunnel entrance and first tier of overburden; b stratigraphy of Bandelier Tuff (8). Visible here is only Unit 1 of the Tshirege Member, overlying the uppermost Ottowi. Stratigraphy in (b) (after Lavine et al. 2003) has been aligned with exposed units visible in (a)

Full size image

3 Method

The processes contributing to the muon attenuation in matter are well known (Leo 1987), so that a measurement of muon attenuation can be used to determine the amount of matter, or “range” R (density times length) traversed by the particles along their path. If R j is the range of a muon traveling along a path j, then

$$R_{j} = \sum\limits_{i} {L_{ji} \rho_{i} }$$
(1)

where each L ji is the path length of muons along path j through the ith matter volume, and ρ i is the density of the ith matter volume. Equation (1) is leveraged to determine the three-dimensional density distribution of an object from the range of muons traveling along different paths through it, via an inverse linear problem.

Both the stopping power and the range R (density times length) for muons have been calculated and tabulated as a function of the initial energy of the muons for different materials and compounds (Groom et al. 2001). The relationship between the minimum energy E min that a muon must possess to traverse a material without being absorbed and its range R is tabulated for a set of elements and materials:

$$E_{\hbox{min} } = E_{\hbox{min} } (R).$$
(2)

If \(f_{\mu } = f_{\mu } (\theta ,E)\) is the energy distribution of cosmic ray muons at the Earth’s surface as a function of their zenith angle θ and of their kinetic energy E, then the ratio of the number of muons N surv surviving the passage through an object to the number N tot of incident muons is given by:

$$\frac{{N_{\text{surv}} }}{{N_{\text{tot}} }} = \frac{{\int_{{E_{\hbox{min} } }}^{\infty } {f(\theta ,E)} }}{{\int\nolimits_{0}^{\infty } {f(\theta ,E)} }}.$$
(3)

Equations (2) and (3), together with the knowledge of f(θ,E) can be used to obtain the range R of the muons along a path from their measured attenuation. We used values tabulated by Groom et al. (2001) for standard rock to model the dependence E min(R) and the probability distribution functions provided by the Cosmic ray Shower Monte Carlo software (Hagmann et al. 2012) to model the energy and angular distribution of cosmic ray muons at the 2100 m elevation of Los Alamos.

4 Data Acquisition and Analysis

The Los Alamos National Laboratory (LANL) Mini Muon Tracker (MMT) (Perry, 2013), shown in Fig. 3, was used for this experiment. The detector consists of two modules made of 576 sealed, aluminum drift tubes arranged in planes. Each of the two modules consists of six horizontal planes of tubes; tubes in each plane are oriented orthogonally to those of its immediate neighbor. The detector is capable of tracking muons with 2.5 milliradians angular resolution across a surface of 1.4 m2.

Fig. 3
figure 3

The Mini Muon Tracker (MMT) deployed inside the tunnel

Full size image

The majority of cosmic rays in the atmosphere are muons and electrons (Olive et al. 2014; Grieder 2001). The flux of muons at sea level is approximately 1 cm−2 min−1, and the electron flux is about 35–40% of the muons. A small hadronic component, consisting of mainly protons and neutrons, is also present and it amounts to ~10% of the total flux at the 2100 m elevation of the Los Alamos tunnel, where the measurement was conducted. Cosmic ray muons cross the 4.3 cm aluminum of the drift tubes with minimal attenuation and produce straight trajectories; we consequently discriminated them from the few electrons surviving the passage through the detector based on their lower track residuals.

Table 1 presents the data acquired at each of the detector positions in the tunnel, as well as outside for obtaining true background flux.

Table 1 Muon data acquisition
Full size table

Muon tracks obtained from the data at each location were divided into 968 angular bins each one having ∆θ = ∆φ = π/44 in the range θ ∈ [0, π/4], φ ∈ [−π,π]. As a result, 968 × 4 = 3872 muon paths were considered.

The detector’s acceptance is limited to a solid angle of ~45° about the vertical axis and becomes marginal for larger inclination angles. Since we used the ratios (see Eq. 3) between the number of muons recorded underground and the number of muons recorded on the surface in each angular bin, any detector acceptance artifacts cancel and do not affect the data.

Figure 4 shows three views of the mesa with the four detector locations and the 45 acceptance cones above them.

Fig. 4
figure 4

Muon acceptance cones for 4 positions of the detector. Red cones indicate limits of the MMT acceptance for muon tracks. a Map view. b View along strike of canyon. c View perpendicular to strike of canyon. Axes labeled in NAVD88 coordinates. Figure courtesy of Alex Johnson

Full size image

To solve Eq. (1) we parameterized the 5.4 × 106 m3 volume above the tunnel as a rectilinear model comprised of 84,672 cubic cells with 4 m sides. Cells residing outside of the detector’s acceptance at all locations were excluded from the inversion to reduce the computational time required to solve the problem; hence, the number of cells considered was 38,546. Equation (1) cannot simply be inverted: the matrix L is singular since this problem is underdetermined in part of the parameter space considered. We, therefore, apply linear regularization.

Different regularization methods exist (Tarantola and Valette 1982; Press et al. 1992); here we adopt a smoothing constraint on the rock density through an exponential covariance (Tarantola and Valette 1982):

$$C_{\text{reg}} (i,j) = \sigma_{\rho }^{2} e^{{ - \left| {r_{i} - r_{j} } \right|/\lambda }}$$
(4)

where σ ρ is the a priori error on the density, λ is a correlation length and |r i −r j | is the distance between the ith and the jth cells. A generalized χ 2 can then be defined as:

$$X_{\text{gen}}^{2} = X_{\text{d}}^{ 2} + X_{\text{reg}}^{ 2} = (R - R_{0} )C_{\text{d}}^{ - 1} (R - R_{0} ) + (\rho - \rho_{0} )C_{\text{reg}}^{ - 1} (\rho - \rho_{0} )$$
(5)

where R is given by Eq. (1), R 0 are the density lengths obtained from the data using Eq. (2) and Eq. (3), ρ 0 is an initial guess on the density of each cell and C d is a diagonal matrix whose (j,j) entry is the standard deviation of the muon range along path j based on the statistics collected along that path. From Tarantola and Valette (1982) Eq. (5) can be minimized with the constraint from Eq. (1) and the solution is:

$$\rho = \rho_{0} + C_{\text{reg}} \cdot L^{\text{T}} \cdot (C_{\text{d}} + L \cdot C_{\text{reg}} \cdot L^{\text{T}} )^{ - 1} \cdot L \cdot C_{\text{reg}} .$$
(6)

The solution of Eq. (6) depends on three parameters, ρ 0, σρ and λ. To optimize them we used synthetic data.

We obtained a LIDAR map of the Los Alamos mesa to reproduce its contours and assigned a uniform density of 1.8 g/cm3 to the rock, and 0 g/cm3 to the air. We then employed a forward model to calculate the expected muon rates and angular distributions at the actual detectors locations, and subsequently applied the inversion algorithm described above.

The value of ρ 0 was chosen equal to 1.5 g/cm3, the average value for the cells in our density model. The values of σ ρ and λ were chosen to maximize the agreement between the density model used and the result while keeping the ratio X 2d /NDF close to 1, NDF being 3872, the number of degrees of freedom for the un-regularized problem:

$$\left\{ {\begin{array}{*{20}c} {\sum\limits_{i = 1}^{{N_{\text{vox}} }} {\left| {\rho - \rho_{\text{synth}} } \right|^{2} }} = {\text {min}} \\ {X_{\text{d}}^{ 2} /{\text{NDF}} = (\left( {R - R_{0} } \right)C_{\text{d}}^{ - 1} (R - R_{0} ))/{\text{NDF}} = 1} \\ \end{array} } \right.$$
(7)

The optimal values of λ and σ ρ were found, respectively, equal to 188 m and 1.2 g/cm3. For these values Χ 2d /NDF = 1.001 and \(\sum\nolimits_{i = 1}^{{N_{\text{vox}} }} {\left| {\rho - \rho_{\text{synth}} } \right|^{2} }\) = 22,862 (g/cm3)2. Figure 5 shows the dependence of the solutions on the two parameters λ and σ ρ when ρ 0 = 1.5 g/cm3, and in particular shows that the solution is quite stable with respect to variations of λ and σ ρ. Figure 6 shows how the quantity \(\sum\nolimits_{i = 1}^{{N_{\text{vox}} }} {\left| {\rho - \rho_{\text{synth}} } \right|^{2} }\) depends on the distance of the cells considered to the straight line running along the center of tunnel, where the muon detector was deployed.

Fig. 5
figure 5

Solution dependence on regularization and density parameters. a Dependence of Χ 2d on λ and σ ρ for ρ 0 = 1.5 g/cm3. b Dependence of \(\sum {\left| {\rho - \rho_{\text{synth}} } \right|}^{2}\) on λ and σ ρ for ρ 0 = 1.5 g/cm3. The white star indicates the values of λ and σ ρ that satisfy Eq. (7)

Full size image
Fig. 6
figure 6

Dependence of solution quality on distance of cells from tunnel axis. \(\frac{{\sum\nolimits_{i = 1}^{{N_{\text{vox}} }} {\left| {\rho - \rho_{\text{synth}} } \right|^{2} } }}{{N_{\text{vox}} }}\) for cells at distances from a line running along the center of the tunnel. Only those cells inside a cylinder of given radius r, shown on the X-axis, were used to calculate the quantity shown on the Y-axis

Full size image

Muon trajectories often did not cross in those cells along the tunnel or immediately around it, so that the agreement between the input model and the solution obtained from it is sub-optimal for regions along, and close to, the tunnel. When more distant cells are considered, the agreement improves and reaches an optimal value for distances in the neighborhood of 65 m from the tunnel. For larger distances, the trajectories of muons become increasingly parallel to one another, thus unable to resolve density anomaly positions.

5 Results

The data were fitted using the method described above, with ρ 0 = 1.5 g/cm3, σ ρ = 1.2 g/cm3 and λ = 188 m. We present the results in Fig. 7. Only those cells having density larger than 0.4 g/cm3 and inside the detectors’ acceptance are shown for clarity. Cross sections and map views are compared against the corresponding LIDAR values. The profile of the mesa is clearly visible in the tomographic results, and the average value obtained for the rock densities of (1.2 ± 0.9) g/cm3 is compatible with the range [1.28 g/cm3, 1.84 g/cm3] for the density of rhyolitic tuffs based on the measurements by Purtymum (1967). The average density determined for air cells is (0.0 ± 0.3) g/cm3. For the purpose of calculating this value, a cell’s assignment to rock vs. air was made by comparing the density profile determined inverting the muon data against the threshold value of 0.4 g/cm3. Topographic contour maps of the mesa for our model vs. LIDAR appear in Fig. 7e, f.

Fig. 7
figure 7

Tomographic results and comparison to LIDAR rock/air interface. Areas outside detector’s acceptance are not shown. a, b Vertical cross sections through model and LIDAR, respectively. c, d Horizontal sections through model and LIDAR, respectively. e, f Model estimation, and LIDAR validation, respectively, for topographic elevation contours in map view

Full size image

6 Conclusions

We demonstrate 3D tomography of a geological structure obtained by an unconstrained inversion of cosmic ray muon data acquired underground. The image reproduces well the profile of the overburden, and the density obtained for its interior is compatible with independent estimates from gravity as well as standard published densities for this lithology.

Data were acquired at four locations along a straight line within the tunnel, mimicking the configuration of a string of small borehole detectors in a horizontal borehole beneath a target reservoir or other body; we thus demonstrate the potential of borehole muon geophysics for imaging subsurface targets. A prototype borehole detector was developed in this project and tested in the tunnel; its measurements, although of somewhat lower resolution, show fidelity to those from the MMT (Bonneville et al. 2017) and prove promising for this nascent technology.