11 Experimental density radiography of Wudalianchi volcano with cosmic ray muons
Y. Cheng a , R. Han a , Z.Li a , J.Li a , J.Li b , W.Gu b , X.Yang a , X.Ouyang b , B.Liao c a Science and Technology on Reliability and Environmental Engineering Laboratory, Beijing Institute of Spacecraft Environment Engineering, Beijing 100094, China b School of Nuclear Science and Engineering, North China Electric Power University, 2 BeiNong Road, Beijing 102206, China c College of nuclear science and technology, Beijing Normal University, 19 Xinjiekou outer St., Beijing 100875, China
E-mail: [email protected] S U M M A R Y Muon radiography is a promising technique to image the internal density structures upto a few hundred meters scale, such as tunnels, pyramids and volcanos, by measuring the flux attenuation of cosmic ray muons after trvaling through these targets. In this study, we conducted an experimantal cosmic ray muon radiography of the Wudalianchi volcano in northeast China for imaging its internal density structures. The muon detector used in this study is made of plastic scintillator and silicon photomultiplier. After about one and a half month observation for the Laoheishan volcano cone in the Wudalianchi volcano, from September 23rd to November 10th, 2019, more than 3 million muon tracks passing the data selection criteria are obtained. Based on the muon observations and the high-esoluiton topography from aerial photogrammetry by unmanned aerial vehicle, the relative density image of the Laoheishan volcano cone is obtained. The experiment in this study is the first muon radiography of volcano performed in China, and the results suggest the feasibility of radiography technique based on plastic scintillator muon detector. As a new passive geophysical imaging method, cosmic ray muon radiography could become a promising method to obtain the high-resoution 2-D and 3-D density structures for shallow geological targets. K
EYWORDS : Muon radiography; Volcanoes; Telescope – 1 – Introduction
Volcanic hazard assessment and risk management remain two very important scientific subjects with heavy implications both on the population safety and economic development. Anticipating future activity of volcanoes requires monitoring of their activity as well as the information on their internal structures. In this contribution, we would like to use a new method which is known as muon radiography to investigate volcano inner structures. The basic idea of muon radiography is that muons will lose energy via ionization, radiation etc when they pass through mediums. Low energy muons will be terminated in the medium if all their kinetic energy is lost. The survival muon rate can be used to infer the average density in the pathway of the muon track. This proposal was brought up about 50 years ago (Luis 1970, Zhou 1987 ). As a results of the technology improvements in muon trackers, especially the plastic scintillator development, muon radiography is frequently applied in area of architecture, tunnel exploration, volcano imaging (Daniele 2014, Ambrosino, 2015, Morishima, 2017, Kunihiro, 2017, Ran 2020, ELENA, 2017, Hiroyuki 2009). We carry out the first volcano muon radiography in Wudalianchi area in China. The Wudalianchi Volcano Field (WDF) is situated in Northeast China. It’s about 1800 km away from the Pacific plate. The two neighbouring volcanoes is Changbaishan and Jingbohu volcanoes, 400 and 600 km away respectively. The most recent eruption of the WDF is around 300 years ago at the Laoheishan and Huoshaoshan volcano cones (Li , – 2 – volcano by muon radiography. A density map reflecting the average density in the muon pathway will be demonstrated. Method × × Error! Reference source not found. ). A high-precision 3-D image of the volcano is derived, the horizontal – 3 – resolution of the drone image is ~1 cm, and the vertical resolution is 10 cm. From the 3D image of the volcano, we calculated the path-length of the muon trajectory based on a stand-alone code. The path length was obtained from the detector point of view, see Figure 6. Figure 1 Left: Drone scan map of the Laoheishan region. The detector is placed close to the sightseeing road illustrated as a red square. Right: Elevation map of the Laoheishan volcano. We can see the relative height seen at the detector sight is ~200m. Also, an elevation image is derived. The elevation of the detector is 400 m. The highest part of the volcano is ~ 600 m. So the relative height of the volcano is ~200 m, not a very huge volcano. Based on the elevation image, we can get the density length image given in the geological coordinate frame whose horizontal axis denotes the azimuth angle and the vertical axis is the elevation angle, see the right plot of Figure 1. The path lengths are represented by different colours. This image presents the full range (elevation and azimuth) path length of the volcano, which provides us a precise estimation and baseline for thickness imaging of the volcano. Furthermore, according to this image, the detector – 4 – is lifted up 20 degrees to do better scanning of the volcano. In this way, the acceptance of the detector is maximized, therefore, the data-taking time can be shortened. There are three stages in the whole data-taking procedure. First the orientation of the detector is towards northwest with 10 degrees to do volcano imaging. And then we rotate the detector to face the open sky where there is no volcano as baseline calibration run. The third rotation is after a quick look at the available data, and we find part of the image of the volcano is missing, therefore, we turn the detector to again to northeast 30 degrees to get a full image of the volcano. Data analysis 3.1
Data Selection
Approximately 3 million muon tracks were acquired at Laoheishan site. Data taking was performed in time period of about one day, storing all the buffer events. In our measurement, the following selection method is used to determine the muon tracks among all the buffer events. Supposing that the coordinates are (x1, y1), (x2, y2), (x3, y3) for plane 1, plane 2 and plane 3, these event lines can be determined by standard fitting procedure with the line cut, |(𝑥 − 𝑥 ) − (𝑥 − 𝑥 )| < 1 |(𝑦 − 𝑦 ) − (𝑦 − 𝑦 )| < 1 where 1 denotes the id difference of the scintillation bars. Here we use very strict selection cuts to estimate the solid angle carefully. – 5 – To get the coordinates ( 𝜃 𝑔 , 𝜙 𝑔 ) in the geological coordinate frame of the muon track, Euler rotation was performed, 𝑅 𝑥 = [1 0 00 cos𝜃 −sin𝜃0 sin𝜃 cos𝜃 ] 𝑅 𝑦 = [ cos𝜃 0 sin𝜃0 1 0−sin𝜃 0 cos𝜃] 𝑅 𝑧 = [cos𝜃 −sin𝜃 0sin𝜃 cos𝜃 00 0 1] where 𝜃 denotes the corresponding rotation angle. For example, if the orientation of the detector is 20 degrees in the elevation and -10 degrees in the azimuth angle under geological coordinate frame, we first rotate 20 degrees with reference to y-axis and then rotate 10 degrees with reference to the z-axis to transform from the detector frame to the geological coordinate frame. In this way, we get the muon track direction in the geological coordinates. Open-sky data analysis
To better understand the data, a simply Geant4 simulation package was developed. the differential muon flux is defined by the modified Gaisser's Formula (Guan, 2015). Compared to Gaisser's Formula which is applied to the high energy range (100 GeV -100 TeV (Lesparre:2010)) and small zenith angle, modified Gaisser's formula given in Eq. 1 corrects the low energy and large zenith angle flux distributions. The modified Gaisser's formula is given by – 6 –
Φ(𝜃, 𝐸) = 0.14 ( 𝐸 ∗ 𝐺𝑒𝑉) −2.7 ( 11 + 1.1𝐸cos𝜃 ∗ ∗ Where 𝐸 ∗ = 𝐸 (1 + 3.64𝐺𝑒𝑉𝐸(cos𝜃 ∗ ) ) cos𝜃 ∗ = (cos 𝜃 + 𝑝 + 𝑝 (cos𝜃) 𝑝 + 𝑝 cos𝜃 𝑝 + 𝑃 + 𝑝 ) with 𝑝 = 0.102573, 𝑝 = 0.068287, 𝑝 = 0.958633, 𝑝 = 0.0407253, 𝑝 = 0.817285. The predicted number of muons is as follows: 𝑁 𝑖𝑛 (𝜃) = ∫ Φ(𝜃, 𝐸). 𝑆 𝑒𝑓𝑓 (𝜃). ε . Δ𝑇𝑑𝐸 ∞ (2) S eff (θ) denotes the effective area of the detector, ε is the total efficiency in the experiment, Δ𝑇 is the data taking time, Φ(𝜃, 𝐸) is the differential muon flux as a function of the zenith angle 𝜃 and the muon energy 𝐸 as shown above. The integration starts from 0.1 GeV, which is taken from experience. – 7 – Figure 2 Solid angle as a function of the x-direction and y-direction id difference. The solid angle is used in the detection efficiency correction for muons with different incoming directions. Figure 3 Muon flux measurement results in case of open-sky (red dot) and volcano (green triangle), and their comparison with the modified Gaisser formula prediction. The low zenith angle discrepancy is due to noise contamination, such as electrons, protons. In the high zenith angle range, we can clearly see the muon flux decrease due to the attenuation caused by the pass length in the volcano. – 8 –
We did a muon flux dependence measurement with reference to the zenith angle, Figure 3. Various corrections were introduced taking into account the actual detector performance. These include the efficiency corrections of the liquid scintillation bars, the zenith angle dependent solid angle (Figure 2) and acceptance correction, also the detection area and time duration correction. Deviations from the modified Gaisser model can be seen in the lower zenith angle part in Figure 3, it's assumed to be caused by fake muon tracks which are indeed electrons or protons. To avoid this problem, it’s proper to do some shielding, for example, using lead plates, to prevent some low energy protons or electrons from penetrating three layers of liquid scintillator bars. The exact flux number of the protons or electrons is hard to estimate, since our detector has no ability to do particle identification.
Volcano Data analysis
If there is a volcano, we need to modify the expected muon number as: 𝑁 𝑖𝑛 (𝜃) = ∫ Φ(𝜃, 𝐸). 𝑆 𝑒𝑓𝑓 (𝜃). ε . Δ𝑇𝑑𝐸 ∞ 𝐸 𝑚𝑖𝑛 (3) where 𝐸 𝑚𝑖𝑛 reflects the minimal kinetic energy that a muon must possess to hit the detector surface without being absorbed, which depends on the average medium density 𝜌 and the length muon passed through 𝑙 . The transmission power of the muon is taken from (Groom, 1999). In this reference, the CSDA (continuously slowing down ability) for various materials can be found. Here we use liquid water as a reference to account for how muons lost its energy via ionization collision or radiation when passing-by liquid water. – 9 – For data analysis, muon transmission 𝐾 𝑎 is defined as the ratio of the muon event rate recorded by the detector in the tunnel and that in the open air under the investigation, which can be calculated by 𝐾 𝑎 = 𝑁 𝑉 (𝜃)/Δ𝑇 𝑉 𝑁 𝑂 (𝜃)/Δ𝑇 𝑂 (4) where subscripts 𝑉 , 𝑂 denote data selected from the volcano direction and the open-sky, respectively. Substituting Eq. (1) into Eq. (4), we have 𝐾 𝑎 = 𝑁 𝑉 (𝜃)/Δ𝑇 𝑉 𝑁 𝑂 (𝜃)/Δ𝑇 𝑂 = ∫ Φ(𝜃, 𝐸)𝑑𝐸 ∞ 𝐸 𝑚𝑖𝑛 ∫ Φ(𝜃, 𝐸)𝑑𝐸 ∞ 𝐸 (5) since 𝑆 𝑒𝑓𝑓 (𝜃) and ε do not change during the experiment. By using Eq. (5) and muon energy spectrum described in modified Gaisser's formula as in Eq. (1). we can easily obtain 𝐸 𝑚𝑖𝑛 . The relationship between the minimum energy 𝐸 𝑚𝑖𝑛 that muons must possess to penetrate the volcano without being absorbed and its density-length 𝑅(𝐸 𝑚𝑖𝑛 ) have been calculated by Groom, et al. (Groom, 1999), where the density of the standard rock is 𝜌 = 2.65𝑔/𝑐𝑚 and 𝑍𝐴 = 0.50 . In this paper, the density of the volcano rock is taken as . The density-length 𝑅(𝐸 𝑚𝑖𝑛 ) is defined by 𝑅(𝐸 𝑚𝑖𝑛 ) = 𝜌 ⋅ 𝐿 . Above all, it is possible to get the functional relationship between
𝑅(𝐸 𝑚𝑖𝑛 ) and muon transmission 𝐾 𝑎 power, ie the muon survival ratio. The calculated results are as follows and can be repeatedly used in future analysis. – 10 – Figure 4 The survival muon ratio (transmission K a ) and the minimum muon energy under different zenith angle. The calculation is based on modified Gaisser formula. We can see that, the energy spectrum for horizontal muons is much harder than the vertical muons. To get the mountain profile via muon radiography, we can directly use the flux ratio between the volcano and open-sky case. The ratio plot is obtained by the ratio of the elevation and azimuth angle dependent flux measurement results when the detector facing the volcano and facing the open-sky. The decrease in muon ratio reflects the average density of the volcano in the specific direction of the muon passage. For example, if the ratio equals 1.0, that means the density is the same, i.e. muons go merely through air in this specific direction. The raw data-taking time normalized muon flux ratio is shown as in top and bottom plot. In this plot, only data-taking time difference is corrected. The detector related correction, such as efficiency correction of scintillation bars, solid angle and acceptance correction can be cancelled out naturally in the case of pixel by pixel muon flux ratio calculation. In Figure 5, we can clearly see the volcano profile represented by the red part. Since the whole data-taking period is composed of three stages, one is northwest 10 degrees, towards the west part of the volcano. The other is open-sky when – 11 – detector faces no volcano, and the third one is facing northeast 30 degree, towards the east part of the volcano. −60 −40 −20 0 20 40 Azimuth (deg) Figure 5 Raw data-taking time normalized ratio plot from the data. The top one is the ratio when detector facing northwest 10 degrees and open-sky, and the west part of the volcano is clearly represented by the small ratio. The bottom plot is the ratio when detector facing northeast 30 degrees and open-sky, similarly the east part of the volcano is clearly represented by the decreased ratio. The blue part in both the top and bottom plot, where the ratio is close to 1, means there is no decrease in muon rate, which means there is no volcano in the passage of the muons in this specific direction. – 12 – From the ratio plot, we can see the impact from the volcano, due to the density length of the volcano, some of the muon tracks with energy smaller than minimum required energy is terminated inside the volcano, thus we can see ratio less than 1. For the specific elevation and azimuth angle, where there is no volcano, we can see the ratio value is very close to 1.0 with small uncertainty. However, we can see that the ratio is much larger than the expectation value (taking into account the drone scan results, where the volcano is hundreds of meters in length, and assume the average density of the volcano is ~ 1.7 g/cm ). Discussions are made in the following. First of all, the noise from other radioactive sources act as backgrounds, which can also pass through all the three layers of the plastic scintillators. Second, in the analysis, we treat the detector as a point. Then it’s possible that parallel muons in certain direction may or may not pass through the volcano. The parallel muons which don’t pass through the volcano also acts as background. However, calculation indicates that parallel muons effects can hardly affect the ratio. Thirdly, considering the geometrical of the detector, the backward muons can also increase the ratio, for horizontal muons, we cannot differentiate which side the muons actually come from, that’s the possible reason why in low elevation angle where the path length in volcano is the largest, such as 0 degree, the ratio increase again. Last but not least, the low energy muons which are scattered at the surface of the volcano, and the changed their direction can also be taken as the high energy muons passing through volcano but with a wrong direction, therefore gives us a wrong ratio value. However, when overlay the ratio plot with the path-length inside the volcano plot, we can see good boundary agreement between these two plots. The path length inside the volcano and the ratio is in good coincidence taking into account the relatively large size – 13 – of the scintillation bars, i.e. ~5 cm, and the various background. The separation between air and the volcano is obvious, with blue part represents the air and the colourful part represents the volcano. After the ratio data analysis, we obtained the profile imaging of Laoheishan volcano, see Error! Reference source not found. . However, we are more interested in the inner structure of the volcano. Density abnormal is one of the key issues we would like to study. The procedure of density unfolding is as follows. First, we have the ratio value with azimuth and elevation angle set. From the modified Gaisser formula, we can calculated the threshold of the muon energy, serving as the minimum energy to pass through the volcano, under this zenith angle. Then according to the stopping power curve from (Groom, 1999), we can infer the path length inside the certain material. In this manuscript, we use the CSDA curve for the liquid water (1 g/cm ), and the derived density length if using CSDA curve from standard rock (2.65 g/cm ) is very similar. Since we already have the path length from the drone scanning data, we can deduce the density directly by dividing the path length. – 14 – Figure 6 Top: the path length inside the volcano calculated from the drone scan data. The white line clearly visible in the plot is the 1meter equal length contour. The maximum path length is close to 1000 meters, which almost stops all the muons. Bottom: the – 15 – overlaid ratio plot and the path length plot. We can see clear boundary separation between air and volcano. Results and Discussions
Following the procedures introduced in the previous section, we can get the density map. However, as a result of the too high ratio value, the unfolded density is too low and has no physics meaning. Here we publish the relative density results, ignoring the exact numbers. The results are shown as Figure 7. Figure 7 The unfolded density plot of the volcano. Due to noise, the ratio value is not correct. And as a result, the unfolded density result is too low to have physical meaning. Here we only publish the relative density. In the plot, we can see some density pattern is presented: the top of the volcano has relatively large density and the bottom of the volcano – 16 – has relatively small density. This is similar to the observation by (Tanaka,2007). They also observed high density region beneath the dome. Since the measured ratio between volcano scenario and open-sky scenario is one order larger than expected, we can only have the relative density map inside the volcano. We raised up the following factors that can influence the measurement. The electronic noise from the SiPM or readout board. The electrons or protons that pass through all the 3 layers and form a straight line, and behaves as a muon track. There are also some back-ward muons, which comes from the side without volcano, especially the incoming direction of the horizontal muons can not be distinguished. The scattered muons is another very important background, it’s muons that passing through the edge of the volcano and scattered by the volcano, and changed its original direction. As a result, the scattered muons smear the angular distribution of the muons after passing through the volcano. Another reason which lead to the unphysical density resutls is that the violent temperature change during the data-taking period. During the data-taking period, the temperature in Laoheishan region goes from 20 Celsius degree to minus 10 Celsius degree. The daily event rate fluctuation is within 10% range at three different experiment stage. SiPM with automatic temperature compensation function is more desirable. To suppress the fake muon tracks, such as electrons and protons, lead plate between scintillator layers is suggested. The lead plate can terminate some low energy electrons or protons. Besides, it can reflect the direction of high energy electrons and protons. Check the direction before and after the lead plate can be used as a method to do electron and muon discrimination. – 17 –
For the back-ward muons, considering the practical condition of the observation field, muons from backward relative to the volcano direction can dilute the effective signal as mentioned before. To evaluate the backward muon noises, we used a Genat4 simulation package to estimate the noise level. The simulation results suggest that muons from frontward versus muons from backward is 100:1.494 in the whole zenith angle range, and contributed hugely in large zenith angle range. The scattered muons can be very annoying and ruined the measurement. (
Nishiyama, 2016 ) has put up a method to do evaluate scattered muon noise with some Geant 4 simulation. Due to the backgrounds and noise mentioned above, using muon radiography in actual volcano internal density measurement can be still challenging and should be considered in R&D process of the muon detector. Conclusion
In this study, we conducted the first muon radiography of volcano in China. In case of volcano, we can clearly see muon flux defecit in the zenith angle range covered by the volcano, ie 70-90. We did the relative ratio measurement and successfully get the profile image of the volcano. The boundary of the volcano is in well coincidence with the drone scaned image. The seperation of air and volcano is clear. We also had the relative density to inspect the inner structures of the volcano. The relative density map reflects that the top of the volcano is high-density region and the bottom of the volcano is low density region. (Tanaka, 2007) has similar observation. – 18 –
In summary, the volcano works as a filter to modulate the muon flux and angular distribution. It requires good background control and evalution, and can be used as a method to derive the inner structure of the kilo-meter scale object.
Acknowledgements
This work was supported by the National Laboratory Foundation (Grant No. 6142004180203).
References
Adloff, C., et al, (2013), Construction and test of a 1×1 m2 micromegas chamber for sampling hadron calorimetry at future lepton colliders,
Nucl. Instrum. Methods Phys. Res., Sect. A, , 90–101. Ambrosino, F., Anastasio, A., Bross, A., Béné, S., Boivin, P., Bonechi, L., Cârloganu, C., Ciaranfi, R., Cimmino, L., Combaret, Ch., et al. (2015), Joint measurement of the atmospheric muon flux through the Puy de Dôme volcano with plastic scintillators and Resistive Plate Chambers detectors, J. Geophys. Res. Solid Earth, , 7290– 7307, doi:10.1002/2015JB011969. Anne Barnoud, et al, Bayesian joint muographic and gravimetric inversion applied to volcanoes,
Geophysical Journal International , Volume , Issue 3, September 2019, Pages 2179–2194, https://doi.org/10.1093/gji/ggz300. Baxter, P. J., R. Boyle, P. Cole, A. Neri, R. Spence, and G. Zuccaro, The impacts of pyroclastic surges on buildings at the eruption of the Soufri`ere Hills volcano,
Montserrat, Bull. Volcanol. , , 293–313, doi:10.1007/s00445-004-0365-7, 2005. – 19 – Daniele Carbone, Dominique Gibert, Jacques Marteau, Michel Diament, Luciano Zuccarello, Emmanuelle Galichet, An experiment of muon radiography at Mt Etna (Italy),
Geophysical Journal International , Volume 196, Issue 2, February, 2014, Pages 633–643, https://doi.org/10.1093/gji/ggt403. D, Z, Zhou, et al. (1987), Physical and Engineering applications for researches of underground cosmic rays,
ACTA GEOPHYSICA SINICA , Volume , Issue 5, Sept., ELENA GUARDINCERRI, et al, 3D Cosmic Ray Muon Tomography from an Underground Tunnel, Pure and Applied Geophysics , (2017), 2133–214, DOI 10.1007/s00024-017-1526-x. Groom, D E, Mokhov, N V, and Striganov, S I. Muon Stopping-Power and Range Tables: 10MeV - 100 TeV. United States: N. p., 1999. Web. Guan, Mengyun, et al, 2015, arxiv: 1509.06176. H. K. M. Tanaka, et al , Cosmic-ray muon imaging of magma in a conduit: Degassing process of Satsuma-Iwojima Volcano, Japan GEOPHYSICAL RESEARCH LETTERS , VOL. , L01304, doi:10.1029/2008GL036451, 2009. H. K. M. Tanaka, et al, Imaging the conduit size of the dome with cosmic‐ray muons: The structure beneath Showa‐Shinzan Lava Dome, Japan, GEOPHYSICAL RESEARCH LETTERS , Volume , Issue 22, 2007, DOI:10.1029/2007GL031389. Katherine Cosburn, et al, Joint inversion of gravity with cosmic ray muon data at a well-characterized site for shallow subsurface density prediction, Geophysical Journal International , Volume , Issue 3, 2019, Pages 1988–2002, DOI:10.1093/gji/ggz127. Kunihiro Morishima, et al, Discovery of a big void in Khufu’s Pyramid by observation of cosmic-ray muons,
Nature , , pages386–390(2017) , doi:10.1038/nature24647. Luis W. Alvarez, et al. (1970), Search for Hidden Chambers in the Pyramids, Science , – 20 – Vol. , Issue 3919, pp. 832-839, DOI: 10.1126/science.167.3919.832. Li, Zhiwei, et al, Shallow magma chamber under the Wudalianchi Volcanic Field unveiled by seismic imaging with dense array,
Geophysical Research Letters , Volume , Issue 10, pp. 4954-4961, doi:10.1002/2016GL068895 Morishima, K., Kuno, M., Nishio, A. et al. Discovery of a big void in Khufu’s Pyramid by observation of cosmic-ray muons. Nature , 386–390 (2017). DOI:10.1038/nature24647. Peter G Lelièvre, et al, Joint inversion methods with relative density offset correction for muon tomography and gravity data, with application to volcano imaging,
Geophysical Journal International , Volume , Issue 3,(2019), Pages 1685–1701, https://doi.org/10.1093/gji/ggz251. Ran, H., Qian, Y. , et al., Cosmic muon flux measurement and tunnel overburden structure imaging,
Journal of Instrumentation , Volume , June 2020, DOI: Journal of Geophysical Research : Solid Earth, (2014), , 699– 710, doi:10.1002/2013JB010234. Ryuichi Nishiyama, Akimichi Taketa, Seigo Miyamoto, Katsuaki Kasahara, Monte Carlo simulation for background study of geophysical inspection with cosmic-ray muons,
Geophysical Journal International , Volume206