3D Cosmic Ray Muon Tomography from an Underground Tunnel
Elena Guardincerri, Charlotte Rowe, Emily Schultz-Fellenz, Mousumi Roy, Nicolas George, Christopher Morris, Jeffrey Bacon, Matthew Durham, Deborah Morley, Kenie Plaud-Ramos, Daniel Poulson, Alain Bonneville, Richard Kouzes
33D Cosmic Ray Muon Tomography from an Underground Tunnel
Elena Guardincerri *, Charlotte Rowe , Emily Schultz-Fellenz , Mousumi Roy , Nicolas George , Christopher Morris , Jeffrey Bacon , Matthew Durham , Deborah Morley , Kenie Plaud-Ramos , Daniel Poulson , Alain Bonneville , Richard Kouzes Affiliations: Los Alamos National Laboratory University of New Mexico Pacific Northwest National Laboratory *Correspondence to: [email protected]
Abstract : We present an underground cosmic ray muon tomographic experiment imaging 3D density of overburden, part of a joint study with differential gravity. Muon data were acquired at four locations within a tunnel beneath Los Alamos, New Mexico, and used in a 3D tomographic inversion to recover the spatial variation in the overlying rock-air interface, and compared with a priori knowledge of the topography. Densities obtained exhibit good agreement with preliminary results of the gravity modeling, which will be presented elsewhere, and are compatible with values reported in the literature. The modeled rock-air interface matches that obtained from LIDAR within 4 m, our resolution, over much of the model volume. This experiment demonstrates the power of cosmic ray muons to image shallow geological targets using underground detectors, whose development as borehole devices will be an important new direction of passive geophysical imaging. ne Sentence Summary:
3D density modeling using an underground cosmic ray muon detector points the way to new borehole methods for passive shallow geophysical imaging.
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 ( ). 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 (
2, 3 ) 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 ( ), 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 ( ) and geologic mapping ( ). 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. xperiment Site Our data were acquired at four locations within a tunnel, excavated horizontally under the main Los Alamos town-site mesa ( ) (Figure 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 o 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.
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.
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. The tunnel entrance is ndicated by a red star. Figure 1b illustrates the configuration of the tunnel and ancillary structures.
At the tunnel site (Figure 2) the Tshirege Member (upper unit) of the Bandelier Tuff is exposed ( ). This member can be further subdivided into cooling units, which have demonstrable mineralogic and physical variations due to the episodic nature of the Tshirege Member eruptive sequence.
Figure 2:
Field site exposure and stratigraphic column. A) Tunnel entrance and first tier of overburden B) Stratigraphy of Bandelier Tuff ( ). Visible here is only Unit 1 of the Tshirege Member, overlying the uppermost Ottowi. Stratigraphy in b) has been aligned with exposed units visible in a). Method
The processes contributing to the muon attenuation in matter are well known ( ), 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 travelling along a path j, then (1) R j = L ji ρ ii ∑ (A) (B) here each L ji is the path length of muons along path j through the i-th matter element, and ρ i is the density of the i-th matter element. 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 ( ). 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: (2) If 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: (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 ( ) for standard rock to model the dependence E min (R) and the probability distribution functions provided by the Cosmic-Ray Shower Monte Carlo software ( ) to model the energy and angular distribution of cosmic-ray muons at the 2100 m elevation of Los Alamos. E min = E min ( R ) f µ = f µ ( θ , E ) N surv N tot = f ( θ , E ) E min ∞ ∫ f ( θ , E ) ∞ ∫ Data acquisition and analysis
The Los Alamos National Laboratory (LANL) Mini Muon Tracker (MMT) ( ), shown in Figure 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 m . Figure 3 : The Mini Muon Tracker (MMT) deployed inside the tunnel . Table 1 presents the data acquired at each of the detector positions in the tunnel, as well as outside for obtaining true background flux. Position
Table 1:
Muon Data Acquisition
Muon tracks obtained from the data at each location were divided in 968 angular bins each one having Δθ = Δϕ = π /44 in the range θ ∈ [0, π /4], ϕ ∈ [ −π , π ] . As a result, 968 x 4 = 3872 muon paths were considered. The detector’s acceptance is limited to a solid angle of ~ 45 degrees about the vertical axis and becomes marginal for larger inclination angles. Since we used the ratios (see Equation 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. To solve equation (1) we parameterized the 5.4 m 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. 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 Megan Lewis. (A) (B) (C) ifferent regularization methods exist ( ); here we adopt a smoothing constraint on the rock density through an exponential covariance ( ): (4) where σ ρ is the a priori error on the density, λ is a correlation length and |r i -r j | is the distance between the i-th and the j-th cells. A generalized χ can then be defined as: (5) where R is given by Equation (1), R are the density lengths obtained from the data using Equation (2) and Equation (3), ρ 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. Equation (5) can be minimized ( ) with the constraint from equation (1) and the solution is: (6) The solution of Equation (6) depends on three parameters, ρ σ ρ and λ . In order 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/cm to the rock, and 0 g/cm 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. C reg ( i , j ) = σ ρ e − r i − r j / λ Χ gen = Χ d + Χ reg = R − R o ( ) C d − R − R o ( ) + ( ρ − ρ ) C reg − ( ρ − ρ ) ρ = ρ + C reg ⋅ L T ⋅ ( C d + L ⋅ C reg ⋅ L T ) − ⋅ L ⋅ C reg The value of ρ was chosen equal to 1.5 g/cm , the average value for the cells in our density model. The values of σ ρ and λ were chosen in order to maximize the agreement between the density model used and the result while keeping the ratio Χ d2 /NDF close to 1, NDF being 3872, the number of degrees of freedom for the un-regularized problem: ρ − ρ synth i = N vox ∑ = min Χ d / NDF = ( R − R ( ) C d − ( R − R )) / NDF = ⎧⎨⎪⎪⎩⎪⎪ (7) The optimal values of λ and σ ρ were found respectively equal to 188 m and 1.2 g/cm . For these values Χ d2 /NDF = 1.001 and ρ − ρ synth i = N vox ∑ = 22,862 (g/cm ) . Figure 5 shows the dependence of the solutions on the two parameters λ and σ ρ when ρ = 1.5 g/cm , and in particular shows that the solution is quite stable with respect to variations of λ and σ ρ . Figure 6 shows how the quantity ρ − ρ synth i = N vox ∑ depends on the distance of the cells considered to the straight line running along the center of tunnel, where the muon detector was deployed. 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. Figure 5: Solution dependence on regularization and density parameters. A) Dependence of Χ d2 on λ and σ ρ for r = 1.5 g/cm . B) Dependence of ρ − ρ synth ∑ on λ and σ ρ for r = 1.5 g/cm . The white star indicates the values of λ and σ ρ that satisfy Equations (7). 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. (m) λ
40 60 80 100 120 140 160 180 200 ) ( g / c m ρ σ × Χ (m) λ
40 60 80 100 120 140 160 180 200 ) ( g / c m ρ σ × synth ρ - ρ ∑ (A) (B) igure 6: Dependence of solution quality on distance of cells from tunnel axis. ρ − ρ synth i = N vox ∑ N 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. Results
The data were fitted using the method described above, with ρ = 1.5 g/cm σ ρ =1.2 g/cm and λ = 188 m. We present the results in Figure 7. Only those cells having density larger than 0.4 g/cm 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/cm is compatible with the range [1.28 g/cm , 1.84 g/cm ] for the density of rhyolitic tuffs based on the measurements by (16). The average density determined for air cells is (0.0±0.3)g/cm . 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/cm . Topographic contour maps of the mesa for our model vs. LIDAR appear in Figures 7E and 7F. Figure 8 shows the difference between the elevation contours obtained respectively from the inversion of the muon data and from the LIDAR scan. The difference is, over most of the range within the detector’s acceptance, smaller than 4 m (the size of the cells, our intrinsic resolution). It becomes larger where the cliff is steeper, and therefore small offsets in the horizontal direction can produce large errors in the vertical direction, and at the edge of the acceptance region, where the density of the muon tracks is lower and the problem is less constrained. Figure 7: Tomographic results and comparison to LIDAR rock/air interface. Areas outside detector’s acceptance are not shown. (A) and (B) Vertical cross section through model and LIDAR, respectively. (C) and (D) Horizontal cross sections through model and LIDAR, respectively. (E) and (F) Model estimation, and LIDAR validation, respectively, for topographic elevation contours in map view.
Figure 8: Difference between the elevation contours obtained from the actual muon data, shown in Figure 7(E) and 7(F)
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 (17) and prove promising for this nascent technology.
References and Notes: K.A. Olive, S. Golwala, P. Vogel, R.Y. Zhu, The review of particle physics.
Chinese Phys. C , doi: 10.1088/1674-1137/38/9/090001 (2014). 2. H.K.M. Tanaka, T. Nakano, S. Takahashi, J. Yoshida, M. Takeo, J. Oikawa, T. Ohminato, Y. Aoki, E. Koyama, H. Tsuji, K. Niwa, 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. Sci. Lett. , 104-113 (2007). 3.
H. K. M. Tanaka, H. Taira, T. Uchida, M. Tanaka, M. Takeo, T. Ohminato, Y. Aoki, R. Nishitama, D. Shoji, H. Tshuiji, Three-dimensional computational axial tomography scan of a volcano with cosmic ray muon radiography,
J. Geoph. Res. (2010). 4.
R. Nishiyama, Y. Tanaka, S. Okubo, H. Oshima, H. K. M. Tanaka, T. Maekawa, Integrated processing of muon radiography and gravity anomaly data toward the realization of high-resolution 3-D density structural analysis of volcanoes: Case study of Showa-Shinzan lava dome, Usu, Japan.
J. Geophys. Res.
D. Bryman, J. Bueno, K. Davis, V. Kaminski, Z. Liu, D. Oldenburg, M. Pilkington, R. Sawyer, “Muon Geotomography — Bringing new physics to ore-body imaging,” in
Building exploration capability for the 21st century , K. D. Kelley, H. C. Golden, Eds., Society of Economic Geologists, Inc, Boulder, CO. (2014). 6.
H. K. M. Tanaka, Muographic mapping of the subsurface density structures in Miura, Boso and Izu peninsulas, Japan.
Scientific Reports , 8305 (2015). . C. A. Rowe, E. Guardincerri, M. Roy, M. Dichter, Joint Tomographic Imaging of 3--D Density Structure Using Cosmic Ray Muons and High--Precision Gravity Data, presented at the American Geophysical Union Fall Meeting, San Francisco, December 2015. 8.
A. C. Lavine, C. J. Lewis, D. K. Katcher, J. N. Gardner, J. Wilson, “Geology of the North-central to Northeastern Portion of Los Alamos National Laboratory, New Mexico,”
Los Alamos National Laboratory report LA-14043-MS (2003). 9.
C. J. Lewis, J.N. Gardner, E. S. Schultz-Fellenz, A. Lavine, S. L. Reneau, S. Olig, Fault interaction and along-strike variation in throw in the Pajarito fault system, Rio Grande rift, New Mexico,
Geosphere ,. 252-269, (2009). 10. W. R. Leo,
Techniques for Nuclear and Particle Physics Experiments . Springer-Verlag, (1987). 11.
D. E. Groom, N. V. Mokhov, S. I. Striganov, Muon stopping power and range tables 10 MeV – 100 TeV.
At. Data Nucl. Data Tables
C. Hagmann, D. Lange, J. Verbeke, D. Wright, “Cosmic-Ray Shower Library (CRY),” Lawrence Livermore National Laboratory document UCRL-TM-229453 (2012). 13.
J. O. Perry, “Advanced Applications of Costmic Ray Muon Radiography ,” Ph.D. Thesis, The University of New Mexico (2013). 14.
W. H. Press, S. A. Teukolsky, W. T. Vetterling, B. F. Flannery,
Numerical Recipes in C . Cambridge Univ. Press, New York (1992). 15.
A. Tarantola, B. Valette, Generalized nonlinear inverse problems solved using the least squares criterion.
Rev. Geophys. Sp. Phys. , 219-232 (1982). 6. W. D. Purtynum,
Geology and Physical Properties of the Near-Surface Rocks at Mesita de Los Alamos , Los Alamos County, New Mexico. U.S. Geological Survey Open-File Report 0180 (1967). 17.
Bonneville, A., et al., 2016, submitted 18.
R. Pordes, D. Petravick, W. Kramer, D. Olson, M. Livny, A. Roy, P. Avery, K. Blackburn, T. Wenaus, F. Würthwein, I. Foster, R. Gardner, M. Wilde, A. Blatecky, J. McGee, R. Quick, The open science grid.
J. Phys. Conf. Ser . , 012057. doi:10.1088/1742-6596/78/1/012057 (2007). 19. I. Sfiligoi, C.D. Bradley, B. Holzman, P. Mhashilkar, S. Padhi, F. Wurthwein, The pilot way to grid resources using glidein WMS. In: , 428-432, doi:10.1109/CSIE.2009.950 (2009). Acknowledgments
This work was supported by the Department of Energy (DOE) Subsurface Technology and Engineering Research, Development, and Demonstration program and by the LANL Center for Space and Earth Science. We used resources provided by the Open Science Grid ( ), which is supported by the National Science Foundation and the U.S. DOE Office of Science. We thank Megan O. Lewis for kindly providing the images in Figure 4. This is Los Alamos Publication LA-UR-16-XXXXXX.), which is supported by the National Science Foundation and the U.S. DOE Office of Science. We thank Megan O. Lewis for kindly providing the images in Figure 4. This is Los Alamos Publication LA-UR-16-XXXXXX.