Cosmic muon flux measurement and tunnel overburden structure imaging
Ran Han Qian Yu Zhiwei Li Jingtai Li Yaping Cheng Bin Liao Lixiang Jiang Sidao Ni Tianfang Liu Zheng Wang
physics.ins-det astro-ph.IM hep-ex

GBKsong





Corresponding authorhanran@ncepu.edu.cn

Science and Technology on Reliability and Environmental Engineering Laboratory, Beijing Institute of Spacecraft Environment Engineering, Beijing 100094, China

Corresponding author zwli@apm.ac.cn
State Key Laboratory of Geodesy and Earth's Dynamics, Innovation for Precision Measurement Science and Technology, Chinese Academy of Sciences, Wuhan 430077, China



Science and Technology on Reliability and Environmental Engineering Laboratory, Beijing Institute of Spacecraft Environment Engineering, Beijing 100094, China


College of nuclear science and technology, Beijing Normal University,100875,Beijing, China


Science and Technology on Reliability and Environmental Engineering Laboratory, Beijing Institute of Spacecraft Environment Engineering, Beijing 100094, China



State Key Laboratory of Geodesy and Earth's Dynamics, Innovation for Precision Measurement Science and Technology, Chinese Academy of Sciences, Wuhan 430077, China


School of Nuclear Science and Engineering, North China Electric Power University, Beijing 102206, China



  We present a cosmic ray muon tomographic experiment for measuring the muon flux and imaging the tunnel overburden structures in Changshu, China. The device used in this study is a tracking detector based on the plastic scintillator with SiPM technology, which can be conveniently operated in field works. The compact system with sensitive area of 6400 cm^2 can measure the angular distribution of cosmic muons. It's able to image the overburden density length from the surface of overburden to the detector along the muon tracks. The open sky muon flux measurement outside the tunnel has a good agreement with the modified Gassier Formula model. The distributions of muon flux at three positions inside the tunnel are very similar to that of open sky. Assuming the average density of overburden compact sandstone is 2.65 g/cm^3, the overburden thickness can be obtained from the density length derived from the difference of muon flux inside and outside the tunnel. Moreover, for known penetrated lengths (i.e., topography of overburden), the density anomalies of the overburden can also been obtained. This study suggests a potential application for imaging and detecting subsurface structures in civil engineering, tunnels or caverns with the cosmic ray muon telescope.
    

    

  Keywords: muon flux, plastic scintillator telescope, muon tomography, overburden structure imaging





Cosmic muon flux measurement and tunnel overburden structure imaging
    Zheng WANG
    December 30, 2023
====================================================================








§ INTRODUCTION



Muons are secondary particles of galactic cosmic rays. Naturally, galactic cosmic rays rain steadily from all directions to the Earth. The interaction of these high energy ions and the atmosphere nuclei produce a cascade of secondary particles. Benefit from the isotropic distribution, muons can be detected at anywhere on the earth in principle with a good approximation. The muon flux reaching the surface of the Earth is about 1 muons/min/cm^2<cit.>.  Muon can penetrate hundreds of meter rocks and can be used to imaging the subsurface structures of large targets because of the high energy and large mass properties <cit.>. This is known as the muon tomography, which is an effective and powerful technique to image the density structures of geological bodies. This kind of measurement is based on the energy loss of high energy muons in matter. Since the muon flux is steeply decreased as a function of the energy, therefore passed matter above the detector modifies the threshold energy for the detectable particles and correspondingly the flux. Due to this fact, the measured muon flux correlates with the density length of traversed material, which is the key issue in the imaging of the internal anomaly structure of large objects<cit.><cit.><cit.><cit.>.

Subsurface structures of the Earth is one of the most important aim in geophysics. However, density structure is still a challenging problem in Geophysics.  With seismological and magnetotelluric measurements, only wavespeed and conductivity structures of the Earth can be obtained. Gavity data can provide constraints on the dnesity structure but with low resolution. Fortunately, Muon is sensitive to the subsurface density structures with high-resolution, which can be used to imaging the density structures reliably. Muon tomography technology has been used in geological studies since 1950 and was first used to explore the overburden density structure above a mountain tunnel <cit.>.  In 2017,the measurement performed by the frame of the METROPOLIS project obtained the evidence of the existence of unexpected cavities in Mt. Echia in Naples  <cit.>. In the same year, Charlotte Rowe and collaborators recovered a 3D tomography of a tunnel beneath Los Alamos, New Mexico <cit.>. The recent study performed by the same group illustrated the potential of inverting for shallow geologic structures by combining gravity and muon measurements <cit.>. Also, this method has been widely applied in the study of the internal structure of volcanoes,such as <cit.><cit.><cit.><cit.><cit.>


Muon tomography is also a very promising technique for subsurface resource exploration and monitoring. In order to facilitate the detection of heavy metal bodies such as uranium, gold or lead, muon measurement experiments have been performed on existing mines<cit.><cit.>. Monitoring tunnel voids can be time-consuming and may involve many personnel working in potentially hazardous environments.  As a non-invasive technique, muon tomography can detecting overburden density anomalies reliablely and accurately<cit.>.


For the future development of muon tomography in the subsurface resource exploration and monitoring applications, we conducted a new muon tomography experiment in a tunnel in Changshu, China. This study introduces the basic structure of the muon telescope, a portable and durable tracking-detector based on plastic scintillator. We described the approach of muon flux measurement and the overburden structure imaging. The detector was placed at three different positions inside the tunnel and one position outside the tunnel, and all the muon flux  measurements results at each point are given. Based on the difference of muon flux between inside and outside the tunnel (i.e. muon flux ratio), we get the overburden structures above the three points inside the tunnel. The experiments well recover the overburden structures with this non-invasive technique, and can obtain the density structures with known overburden thickness distribution, or obtain the overburden thickness with known density distribution.



§ THE EXPERIMENT SITE AND MUON TELESCOPE


The tunnel experiment is located in the northwest of Changshu, Jiangsu Province. The tunnel has a length of about 250m , a height of 2.35m and a width of 2.5m. The temperature inside of the tunnel is around 17.3 degrees Celsius all year round. The tunnel bedrock is compact quartz sandstone with the average density of 2.65 g· cm^-3, which is a dense and hard rock. The average overburden thickness is around 20m, and ranges from 5m(at entrance) to 40m.  The picture of tunnel entrance are shown in Figure. <ref>.





The muon telescope is a three-layer plastic scintillator detector as shown in Figure. <ref>. We designed it in a light-weighted and portable way to use it in outdoor situations. Each layer contains two planes, each plane is made up of 16 plastic scintillator strips, so there are 96 strips in total. The two planes in each layer is placed orthogonally, to provide x and y coordinates.

The size of the plastic scintillator strip is 80cm × 5cm × 1cm, which makes the total effective area is 80cm × 80cm. Each plastic scintillator is covered by highly-reflective enhanced specular reflector (ESR) to improve the scintillation light collection and two silicon photo-multiplier (SiPM) at the end of each strip to collect the light.
The distance between each layer can be adjusted and the maximum distance between the top and bottom layer is 100cm. For the measurement at Changshu tunnel, the distance between the each layer is set to be 50cm.


The model of plastic scintillator is HND-S2 made by Chinese factory, which has relatively high light yield, 50%-60% compared with paranaphthalene. The wavelength of the emission scintillation photon is typically in the 395 nm to 425 nm range. The attenuation length is larger than 2m and the decay time is around 2.8 ns and the refractive index of this scintillator is 1.58, which makes it totally compatible with regular SiPM. Other character such as the H/C ratio is close to 1.1 , the density is 1.05 g· cm^-3  and the softening temperature is 348K - 353K.  The SiPM used is MicroFC-30035-SMT provided by SensL company. The wavelength at the highest quantum efficiency is 420 nm, which matches the emission spectral of the plastic scintillator. It also has other nice features, such as low operating voltage, excellent temperature stability, robustness, compactness, output uniformity, etc.






Data from each strip is processed using the front-end readout electronics (FEE). FEE is equipped with 6 discriminators and 1 FPGA. The discriminator is triggered logical OR signals from the two SiPMs at the ends of strips. After the signal was trigged, the discriminator amplifies the signal and if it exceeds the 32mV threshold, the strip will be considered as fired one, we keep the event with more than 3 strips fired from different layers. The 32 channel signals of each layer array plane are processed by one discriminator. FPGA receives 32× 6=192 digital trigger signals from 6 discriminators and the data is transmitted to a computer via a mini-USB cable.



§ MUON FLUX MEASUREMENT


The muon detector has been deployed from April 22 to July 8, 2019 at the tunnel. Apart from the handing and debugging time for the detector, the tunnel data-taking campaign consisted of up to 73 days. The system was operated continually in the site for everyday. Before operation in the tunnel one open-sky measurement was taken for later detector calibration. Data taking was performed in runs of about 1 hour with at least 3 strips fired, storing all the buffer events.  We successfully took 71 days data in the tunnel and 2 days data at the entrance of the tunnel for the open sky condition. Approximatively 3.01× 10^7 events were obtained in the tunnel and 9.06× 10^6 events outside the tunnel. Table. <ref> presents the data acquired at each of the detector positions. In this section, the analysis procedure to get muon flux at each point is introduced.





 §.§ Data Analysis


The data analysis starts with the reconstruction of the muon tracks event by event. Muon track can be reconstructed in cartesian coordinates (x, y, z) in the detector system, where x and y represent the coordinates in the detection plane while z axis is perpendicular to the detection plane.  In the following analysis, the muon track deflection angle caused by scattering is ignored. The track reconstruction uses only events with coincident hits in all six layers of the detector, and keep events with only one hit per layer.

The reconstruction algorithm preselects coincident hits in a 20 ns time window in all six detection layers, from which tracks are obtained analytically through a straight line procedure. The muon track reconstructed by six layers is much cleaner which could eliminate the fake muon event of other particles or electronic noise striking on two matrices at the same time. In this sample, we select events containing just a single muon track which will lose some of real muon events, but we kept all those marked as “golden track".

We define the pre-selection efficiency as the number of events that satisfy the straight line procedure divided by the number of events with 3 strips fired. At the end of this selection, we estimate an average data selection efficiency around 21.0%.  During the whole data taking period to the open sky, the golden track rate is 11.34 Hz, for the other 3 points inside the tunnel, the rate as follows: 1.30 Hz, 1.04 Hz, 0.85 Hz, respectively. The correspond golden events number are listed in the last column of Table.<ref>. Figure.<ref> gives the map of muon event rate of golden track measured at each location. The muon event rate of selected golden track with θ_x and θ_y distribution as Figure.<ref>. θ_x or θ_y denotes projection of a golden muon track zenith angle in the x-z or y-z telescope plane.







From the selected golden muon track to muon flux measurement, the detector performance and geometrical factor still need to be considered.  A useful quantity for characterizing the behavior of a detector meant for measuring muon flux is the effective surface, S_eff, which takes into account the detection performance and the geometrical acceptance of the detection setup. The geometrical acceptance A(θ_x,θ_y) includes the effective area of the detector, solid angle and the geometrical factor. The detection performance ϵ(θ_x,θ_y) includes the detector efficiency and selection efficiency. Having the number N_g of golden muons detected for each direction (θ_x,θ_y) during a period Δ T , the integrated muon flux I will be Equation.<ref> as


    I(θ_x,θ_y) = N_g/S_eff(θ_x,θ_y)×Δ T(sr^-1· cm^-2· s^-1 )


where


    S_eff(θ_x,θ_y) =  A(θ_x,θ_y) ×ϵ(θ_x,θ_y)



The use of detector matrices formed with a limited number of pixels of finite size causes geometrical acceptance A(θ_x,θ_y). The acceptance of the telescope is a key parameter since the goal of the project is to assign an opacity and therefore a density to the target from an absolute measurement of flux. The acceptance of the telescope includes the area of the detector plane, solid angle of each pixel Ω(θ_x,θ_y) and the geometrical factor Λ(θ_x,θ_y).

 The quantization of solid angel Ω(θ_x,θ_y) depends on the number of pixels N_x× N_y, on their size d, and on the distance D of separating the two detector matrices. In our case, N_x=N_y=16,d=5cm and D=100cm. Typical solid angel curves in this experiment are shown in the left of Figure.<ref>. Taken as an angular resolution, solid angle can privde the spatial resolution of one telescope. For example, with D = 100 cm and for the (0,0) vertical direction, the solid angle is about 2.5×10^-3 sr^-1, this allowing to detect heterogeneities with a size of about 10cm at a distance of 5m, also means the spatial resolution of our telescope can reach 10cm at 5m distance to the observation.

The geometrical factor just in the case of the muon telescopes made of parallel layers of finite dimensions, the particles are
required to physically cross the outer layers in order to induce coincident signals. This factor can be analytically calculated for a perfect
telescope with square layers of size d^2 and a distance D between the two outermost layers. In the right of Figure.<ref>, we show the geometrical acceptance of the detector for each detection element. From this figure we can see that the largest effective area corresponds to the vertical area for every detection element which can be as large as ≃ 16 cm^2· sr^-1. However, the smallest effective areas  corresponds to the body diagonal directions.
For large (θ_x,θ_y), A(θ_x, θ_y) is close to 0, which means only a few muon events can be recorded in these regions. In other words, the effective detection region lie in the range between -36 deg and 36 deg in the detector frame.



The detection performance ϵ(θ_x,θ_y) includes the detector efficiency and selection efficiency.
The detector efficiency D_eff calculated strip by strip. Only the vertical events were used to do the detector efficiency calibration. Figure.<ref> are the efficiency with the strip number, except the 40th strip with a little bit low efficiency, all the others are good enough and the average efficiency is around 93.5%. The final straight line efficiency will be get from all the linear six hits, this gives the average detection efficiency is around 64.3% for the muon track. About the selection efficiency, benefit from that the detector are stable in time based on scintillator, an average trigger efficiency of 80% is estimated.






 §.§ Muon Flux Measurement Results 


 Based on Equation.<ref>, as a function of (θ_x,θ_y), the muon flux is estimated by dividing the number of golden tracks in each angular bin to S_eff(θ_x,θ_y) and a period Δ T in Table.<ref>. The measured muon flux, averaged over bins in the control regions is represented as a function of the (θ_x,θ_y) in Figure. <ref>.


The muon flux in vertical direction(θ_x=θ_y=0 direction, i.e. Zenith angle θ=0) of four positions are shown as below:

  I_(0,0)=9.98× 10^-3cm^-2· s^-1· sr^-1  (Open sky)

  I_(0,0)=1.03× 10^-3cm^-2· s^-1· sr^-1  (P1)

  I_(0,0)=0.83× 10^-3cm^-2· s^-1· sr^-1  (P2)

  I_(0,0)=0.71× 10^-3cm^-2· s^-1· sr^-1  (P3)



  Muon flux distribution over the zenith angle θ is also taken into account at each point, dependency on azimuth angle ϕ being uniform. Results of our experiment are shown in Figure.<ref>. It is seen from Figure.<ref> that the angular distributions is satisfactorily coincident with the modified Gassier Formula <cit.>. The bar in each point gives the range of muon flux with the same zenith angle θ but different azimuth angle ϕ. Zenith angle θ dependency at different overburden structures, which is of interest by itself<cit.>. As it is seen from results, the distribution of position 1, position2 and position 3 turn out to be very similar to the distribution at the open sky data only with less muon flux number, contrary to the frequent expectation that they should become narrower, since the absorption length increases with angle.




§ OVERBURDEN STRUCTURE RESULTS


 The main aim of our measurements in the tunnel system was to verify the applicability of the muon telescope for underground structure, i.e., the surface reconstruction and/or estimate rock inhomogeneity. Here we investigated an underground structure method using cosmic muons. In this section, we will show the approach how to get the overburden density length structure,ρ· L from cosmic muon measurement. In the three measurements, the detector is located at the same elevation inside the tunnel. It is feasible to estimate the average density of the overburden in different directions along a specific line of sight of the detector with known overburden thickness. Moreover, if the average density is known, the overburden thickness distribution can also be obtained.

 The relation between the muon flux and density length is derived from the theory equation firstly. By comparing the experiment results with this relation,the corresponding density length in different zenith and azimuth angle for each experiment site can be calculated. Here we introduce the method to get the relation between muon flux and overburden thickness.
 The number of muons expected to be acquired inside the tunnel N(θ) can be written as Equation.(<ref>)


    N_in(θ) =  ∫^∞_E_minΦ(θ, E)  dE · S_eff·Δ T_in,

where Φ(θ, E) is the differential muon flux as a function of the zenith angle θ and the muon energy E, this paper using modified Gaisser Formula <cit.> model, E_min is the minimal kinetic energy that a muon must possess to reach the detector location without being absorbed by the overburden. Δ T_in is the data taking time and  S_eff is the calibration factors as described in the above section. By comparing the measured muon flux at different positions, one can get the E_min at each point. Using the relation between E_min and density length in Table.EIV-6 of Ref. <cit.> to get the density length of the overburden.

Usually, it is difficult to get precise measurement muon flux, since from measurement muon number to muon flux, we need to take care calibration factors S_eff very carefully. To eliminate the detector's performance influence, the muon flux ratio will be often used in the relation derivation. This ratio K^a is defined as the muon flux in the tunnel to the muon flux in the open sky as,

    K^a=N^ν_in(θ)/Δ T^ν_in/N^ν_out(θ)/Δ T^ν_out,

where subscript in, out denote data selected inside the tunnel and outside the tunnel(the open sky data). From Equation.<ref> we can see that the S_eff is canceled, only need to normalized to the data taking time. The muon flux ratio with different density length <cit.> are plotted in Figure.<ref>.





 Compare the measurement K^a in our experiment and Figure.<ref>, the measurement density length results of each points at the tunnel are shown in Figure. <ref> in detector system.




We expressed the same density length distribution but with (X,Y,Z) in geodetic coordinate system as the top of Figure.<ref>. With the known topography of the tunnel, one can get the density distribution from the density length. In turn assuming an average rock density, the thickness of the overburden will be got. The bottom of the Figure.<ref> show the thickness distribution assuming the average rock density is 2.65g/cm^3. From the figure we can see that the overburden thickness above each position does not change rapidly. The range of the overburden thickness above three positions are around 20.00 ± 2.95 m, 23.03± 2.75 m, 25.18 ± 4.44 m respectively. The standard derivation shows the variation range of the measured overburden thickness in different zenith and azimuth angles. Our results are consistent with the mountain topography above the tunnel.
  



§ CONCLUSION

 

The cosmic ray muon experiments inside and outside the tunnel have been conducted to imaging the overburden structure. Cosmic ray muon tomography has been proved to be a feasible and powerful technique for the identification of void space and significant density anomalies of subsurface structures. Three measurements have been performed inside the tunnel of Changshu area with overburden thickness of 15∼30 meter rock equivalent. Muon tomography can also be performed at larger depths deep to a few hundred meters, at the cost of increased measurement time due to the limited acceptance and strong absorption of muon. The current muon telescope is based on scintillator and SiPM technology with reliable tracking performance and low power consumption. The study demonstrates the feasibility for subsurface structure imaging with cosmic ray muon. The muon tomographic approach is promising in underground tombs investigation, subway tunnels engineering, geophysical, or industrial applications. With finer segmented scintillator detector, and observations from different positions and angles around the study target, a tomographic inversion with higher resolution can be conducted to obtain the three-dimensional density structures.




This work is supported by National Natural Science Foundation of China(Grant No. U1865206), State Key Laboratory of Geodesy and Earth's Dynamics and National Key Laboratory Foundation (Grant No. 6142004180203).





0

Hagiwara:2002

Hagiwara, Kaoru, et al, Review of Particle Physics: Particle data group,Physical Review D 66.1 I (2002) 100011-10001958.

Bogdanova:2006
Bogdanova L N, Gavrilov M G, et al, Cosmic muon flux at shallow depths underground, Physics of Atomic Nuclei 69(8)  (2006) 1293-1298.

Bryman:2014
Bryman D, Bueno J, et al, Muon geotomography bringing new physics to orebody imaging, Society of Economic Geologists (2014).

Lesparre:2010
Lesparre, N., et al. Geophysical muon imaging: feasibility and limits, Geophysical Journal International 183.3 (2010) 1348-1361.
Tanaka:2001
Tanaka, H., Nagamine, K., et al.,  Development of the cosmic-ray muon detection
system for probing internal-structure of a volcano, Hyperfine Interact.,138 (2001) 521-C526
Tanaka:2003
Tanaka, H., Nagamine, K., et al, Development of a two-fold segmented detection
system for near horizontally cosmic-ray muons to probe the internal structure of a volcano, Nucl. Instrum. Methods A507 (2003) 657-C669.

Ceorge:1955
 E. P. George,Cosmic reys measure of overburden of tunnel,Common wealth Engineer,(1955) 455-1955.

Saracino:2017
Saracino, G., et al. Imaging of underground cavities with cosmic-ray muons from observations at Mt. Echia (Naples), Scientific reports,7.1 (2017) 1181.

Guardincerri:2017
Guardincerri, Elena, et al. 3D cosmic ray muon tomography from an underground tunnel, Pure and Applied Geophysics174.5 (2017) 2133-2141.


Cosburn:2019
 Cosburn, Katherine, et al. Joint inversion of gravity with cosmic ray muon data at a well-characterized site for shallow subsurface density prediction, Geophysical Journal International 217.3 (2019) 1988-2002.

Beauducel:2008
Beauducel, F., et al. Muon radiography of volcanoes and the challenge at Mt. Vesuvius presentation (2008).


Marteau:2012
Marteau, J., et al. Muons tomography applied to geosciences and volcanology, Nucl. Instrum. Methods 695 (2012) 23-28.


Fehr:2012
Fehr, Felix, and Tomuvol Collaboration. Density imaging of volcanos with atmospheric muons, Journal of Physics: Conference Series, 375, No. 5(2012).

Tanaka:2009a
Tanaka, H., et al, Cosmic-ray muon imaging of magma in a conduit: degassing process
of Satsuma-Iwojima Volcano, Japan, Geophys. Res. Lett., 36, L01304 (2009) doi:10.1029/2008GL036451.

Tanaka:2009b
Tanaka, H. et al.,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(2009)doi:10.1029/2009GL039448.

Bryman:2015
Bryman D, Bueno J, Jansen J. Blind test of muon geotomography for mineral exploration, ASEG Extended Abstracts, (2015) 1-3.

Baccani:2019
Baccani G, Bonechi L, Bongi M, et al, Muon Radiography of Ancient Mines: The San Silvestro Archaeo-Mining Park (Campiglia Marittima, Tuscany), Universe  5(1) (2019) 34.

Thompason:2019
Thompson L F, et al.,he application of muon tomography to the imaging of railway tunnels. arXiv:1906.05814.

Lesparre:2012a
Lesparre, N., et al. Design and operation of a field telescope for cosmic ray geophysical tomographyGeosci. Instrum. Method. Data Syst.1 (2012) 33-42.


Guan:2015
Guan M, Chu M C, Cao J, et al. A parametrization of the cosmic-ray muon flux at sea-level,arXiv:1509.06176, 2015.
Groom:2001
Groom, Donald E., Nikolai V. Mokhov, and Sergei I. Striganov. Muon stopping power and range tables 10 MeV-100 TeV, Atomic Data and Nuclear Data Tables 78.2 (2001) 183-356.















