Muographic data analysis method for medium-sized rock overburden inspections

Muographic measurements of rock overburdens are of particular interests because they can be applied to natural resources and undiscovered cave explorations, and even to searching for hidden chambers in historic architectural structures. In order to derive the absolute density distribution of the overburden, we conventionally needed to know the accurate information about the measurement conditions, e.g., the detector’s geometrical acceptance, detection efficiency, measurement time, etc. in order to derive the absolute value of the transmitted muon flux. However, in many cases, it is not a simple task to accurately gauge such conditions. Open-sky muon data taken with the same detector are useful as reference data to cancel these factors, however, if the detector is not transportable, this data taking method is not feasible. In this work, we found that the transmitted muon flux will follow a linear function of the areal density along the muon path as long as the incident muon energies are below a few hundred GeV. Based on this finding, we proposed a simple analysis method that does not require detailed knowledge of the detector’s conditions by combining the independently measured density information for the partial volume of the target. We anticipate that this simple method is applicable to future muographic measurements of rock overburdens. Geosci. Instrum. Method. Data Syst. Discuss., doi:10.5194/gi-2016-6, 2016 Manuscript under review for journal Geosci. Instrum. Method. Data Syst. Published: 15 March 2016 c © Author(s) 2016. CC-BY 3.0 License.


Introduction
Earth's subsurface density structures have been extensively measured with muography in the last decade. High energy muons originating from cosmic rays have strong penetrative power to resolve the density distribution of gigantic objects. The muography technique was first applied by Alvarez et al. (1970) as a method to search for hidden chambers inside the second pyramid of Chephren. Muon detectors installed inside the Belzoni chamber recorded the arrival angles of comic ray muons after they had penetrated the Pyramid and reached the chamber. Muography utilizes the natural tendency of highenergy muons to penetrate gigantic materials in the same way that the x-ray does; however, while the x-ray is applicable only to objects up to 1 m in thickness, muography visualizes density distributions inside much larger objects. Since George(1955) first proposed and implemented densimetric measurements of the rock overburden above a gallery tunnel of the Snowy Mountains Hydro Electric Authority in Australia with cosmic ray muons, particle detection techniques had been developed to the extent that 12 years later, Alvarez's group could successfully apply spark chambers with digital readout units to the technique of muography. After analyzing data collected during several months of operation, they concluded that areal density of the pyramid was measured with a precision of 2% along the muon paths that penetrated over 100 meters through limestone in the Pyramid. This pioneering experiment was a crucial step that eventually led to recent muographic experiments that explored inside volcanoes (Tanaka et al., 2007;Tanaka et al., 2008;Tanaka et al., 2009;Lesparre et al., 2012;Carloganu et al., 2012;Carbone et al., 2013;Kusagaya and Tanaka, 2015a;Kusagaya and Tanaka, 2015b), industrial plants (Tanaka, 2013;Ambrosino et al., 2015), seismic faults (Tanaka et al. 2011;Tanaka et al., 2015), and caves (Caffau et al., 1997;Barnafoldi et al., 2012;Oláh et al., 2012).
Cosmic ray muons are generated in the Earth's atmosphere as secondary cosmic rays and the integrated vertical open-sky muon flux ISKY is known to be ~70/m 2 ▪sr▪s. Since uncertainties in these theories (Groom et al., 2001). Once muons are irradiated to the surface of the rock over burden, they are detected and recorded by the detector underneath the overburden to generate the histogram of the number of muon events as a function of the muon's arriving angle.
Muography has the capability to derive an areal density along the muon path. George (1955) compared the muon flux inside and outside of the Guthega-Munyang tunnel, Australia to calculate the muon's transmission rate, and measured the areal density of the rock overburden of 163±8 m water equivalent (m.w.e.) which was consistent with the result of the drilling and sampling at the same site: 175±6 m.w.e. However, the open-sky muon counts are not always available as reference data particularly when the observation system is not transportable. For example, Alvarez et al. (1970) constructed their observation system inside the Belzoni chamber of the Chephren Pyramid. Their entire system weighed more than 10 tons, and thus it was not realistic to transport the system to the outside of the Pyramid to take the reference data without changing the configuration.
In this work, we proposed a new method of muographic data analysis for deriving the density distribution of the middle scale rock overburden (up to a few hundred meters in thickness). With our proposed technique, by combining the value of the density independently measured for a partial volume of the target, we can derive the areal density along the muon path without detailed knowledge about the muon detector. In this paper, we formulated this method and re-analyzed the data presented by Alvarez et al. (1970), Caffau et al. (1995), andLiu et al. (2012) as examples in order to test the method.

Cosmic-ray muons
High-energy muons are produced when primary cosmic rays interact with nuclei in the Earth's atmosphere. The energy spectrum of the primary cosmic ray is expressed with the following power law that has an index of -2.7 (= 1.7) and is almost constant with energies up to 10 6 GeV.
particles of the primaries. Up to this date, several authors have derived the analytical expressions for the differential atmospheric muon spectra in the reaction between primaries and atmospheric nuclei (Bull et al. 1965;Matsuno et al. 1984;Bugaev et al. 1998;Gaisser and Stanev 2008). The parameters in their models have been adjusted by comparing these with the observed spectrum of muons taken from experiments made at sea level (Jokisch et al. 1979;Matsuno et al. 1984;Allkofer et al. 1985;Haino et al. 2004;Achard et al. 2004). The integrated vertical muon flux is proportional to E -2.2 in the energy region between 50 and 200 GeV.
The geomagnetic deflection depends on the geomagnetic latitude, and has a significant effect on the muon flux with energies up to 10 GeV (Haeshim and Bludman 1988) while Hansen et al. (2005) reported that the east-west effect is negligible in the vertical cosmic ray muon flux. Kamiya et al. (1976) reported that the muon charge ratio measured by the MUTRON spectrometer had to be modified with correction factors of 1.35 and 0.75 for muons in the N-W and S-E direction, respectively, indicating a strong geomagnetic effect even on high-energy horizontal muons.

Muon range
The processes through which muons interact with matter can be divided into two types: continuous and stochastic. The rate of the muon's energy loss through matter is expressed where E is the muon's energy and X is an areal density along the muon path. a(E) is the energy loss caused by the continuous process, and b(E)E is the energy loss caused by the stochastic process.
Via the continuous process, the muon makes frequent encounters with atoms, each losing a very small fraction of its energy via the ionization process. Fluctuations in range arise from stochastic processes: bremsstrahlung, direct pair production, and photonuclear interactions. In these processes, the muon loses a large but random fraction of its energy.
However, if the muon's energy is much lower than the critical energy, 708 GeV in SiO2, the continuous process is the main process in muon's energy loss, and thus muons of a given energy would have almost a unique range. For example, an energy loss rate of 100-GeV muons via the continuous and stochastic processes is 2.46 MeV cm 2 g -1 and 0.11 MeV cm 2 g -1 , respectively. 100 GeV muons have a range of 406 m.w.e.
The cutoff energy (Ec) is defined as the minimum energy that a muon can penetrate for a given areal density of rock (X) and is derived by finding the value at which the muon's continuous slowing down approximation (CSDA) range matches X (Fig. 3). The Ec solely depends on an areal density along the muon path in matter, and its material dependence is small. In order to derive the muon flux after passing through matter, the zenith-angular dependence of the muon flux and the continuous slowing down approximation (CSDA) range (Groom et al., 2001) are utilized. The CSDA range is derived by integrating Eq.
(2) over a muon's energy E within a range between 0 and the incident energy E0. Inversely, once an areal density is determined along the muon paths, the minimum energy that a muon can penetrate through the path is uniquely determined if the muon's energy is much lower than the critical energy.

Results
In this work, we modeled muographic observations of the rock overburden. We divided the overburden into two horizontal layers (Layer 0 and Layer 1). The muon detector was assumed to be located underneath this overburden. The areal density was then calculated by assuming a uniform density along the muon paths, and multiplying it with the geometrically exploited muon's path lengths as follows: where r0 and r1 are assumed densities for the Layer 0 and Layer 1, respectively. Both  0 and  1 are the path lengths averaged over the elevation angle () and azimuth angle (),according to an angular resolution of the muogram when it is generated. Most muons traverse matter in a linear trajectory. Therefore, the path length of muons can be precisely determined by reading the outer geometry of the target volume. To calculate the CSDA range, the muon's ionization, bremsstrahlung, direct pair production and photonuclear processes are considered.
In order to derive the transmitted muon intensity (I) (the muon flux after passing through the overburden), the number of muon events counted in each histogram bin (N) is divided by the detector's active area (S) detection efficiency (Aeff) solid angle () and measurement time (t). I is compared with the theoretical muon intensity (i) to derive the areal density along the muon path (X), where i is calculated by integrating the muon energy spectrum over a range between the muon's cutoff energy (Ec) and infinity. In this work, we proposed the flux ratio I0/I1 to cancel S, Aeff, , and t, where I0 and I1 is the measured intensity after passing through Layer 0 and Layer 1 respectively.
When the muon's energy is less than 100 GeV , corresponding to a muon range of 400 m.w.e., the contribution of the stochastic energy loss (b(E)E) is less than 10% to the total energy loss. Furthermore, the energy loss rate due to the ionization process (a(E)) does not have a strong energy dependence. For example, a(E) = 2.3 MeV cm 2 g -1 when the muon energy is 20 GeV (that corresponds to the muon range of 90 m.w.e.), and it becomes 2.5 MeV cm 2 g -1 when muon energy is 200 GeV (that corresponds to the muon range of 740 m.w.e.). Therefore, we can expect that the CSDA range can be fit by a near-linear function within this energy range.
In order to confirm this assumption, the CSDA range reported by Groom (2004) was fit by the following equation.
where X is the CSDA range. A and e are parameters. The results of the fitting 5 data points (40,80,100,140,and 200 GeV) are shown in Eq. (5).
The fitting accuracy was about 0.7%. Within this energy range the muon range varies from 180 m.w.e. to 740 m.w.e. As can be seen in Eq. (5-2), the cutoff energy (Ec) has an almost linear relationship with the areal density. Therefore, we approximated the muon flux ratio (I0/I1) as follows: where I0 and I1 is the transmitted muon flux after passing through Layer 0 and Layer 1, respectively. If the muon counts (N) contain background events (n), i.e., N0=n0+n and N1=n1+n, the muon count ratio N0/N1 will not be equal to the actual flux ratio I0/I1, where n0 and n1 is the number of events without background events, respectively. However, we can reasonably assume that n>> n when the rock overburden above the detector is thicker than 20 m. This thickness is equivalent to 100 times longer than the electron's radiation length in matter. Also, the hadron's interaction length is up to 100 gcm -2 which is 50 times shorter than this thickness. Therefore, we can expect that electromagnetic and hadronic components would be effectively removed by the time they reach the detector. This assumption is supported by a result obtained in the similar experiment in 1955. George (1955) installed a Geiger counter into the Guthega-Munyang tunnel, where the rock overburden was previously measured to be an areal density of 175±6 m.w.e. with the drilling and sampling metrhod. He simply compared the counting rate measured inside and outside the tunnel without any background treatments to obtain an areal density of 163±8 m.w.e. that is in agreement with the muographically derived areal density.
Therefore, it is reasonable for us to approximate Eq. (6) as as long as the overburden thickness exceeds 20 m.

Discussion
In this section, we examined our technique by applying it to three different kinds of targets: the Pyramid of Chephren in Egypt, the limestone cave in Italy called Grotta Gigante, and the Price 5 deposit at the Myra Falls mine, Canada. For all of the targets, the muography data were collected in the past, and had already been published. These data were used for testing our technique.

Case study 1: the Pyramid of Chephren
The Pyramids of Giza are not geological products. However, their sizes are remarkably large for man-made architecture, and the observational configuration of muography performed by Alvarez et al. (1970) was essentially the same as other types of muography performed inside tunnels to measure the rock overburdens (Olah, 2012;Caffau, 1995;George, 1955). The Pyramid of Chephren offers us a unique target volume to test our technique for the following reasons: (A) the geometrical shape of the Pyramid is much simpler than regular geological targets and its topographic features have been studied with aerial surveys (Alvarez et al. 1970 made of Tura-limestone remain on the upper region of the Pyramid (Fig. 2). Below the lower border of the existing casing stones, all casing stones as well as several back stones (those stones located behind the casing stones) have been lost, and as a result the subsurface structure of the Pyramid is exposed. This subsurface exposure allows us to see a partial cross section of the interior structure and reveals that both casing stones and back stones are tightly packed without spaces and one layer of casing stones originally covered the surface of the entire Pyramid (Fig. 2b). The Pyramid was probably intact until at least until 100 B.C. since Diodorus Siculus wrote in that year that, "the stones remain to this day still preserving their original position and the entire structure undecayed" (Siculus and Oldfather, 1933). Therefore, trusting this description, we can assert that in the intervening 2100 years, the subsurface materials (casing and back stones) of the Chephren Pyramid were lost. Croci and Biritognolo (2000) reported, based on their close observations, their conclusion that the missing materials had been systematically removed in a manner similar to quarrying techniques.
Two different types of limestone were used to build Chephren's Pyramid: Tura and Mokattam. Tura limestone (high quality stone similar to marble) was used for the decorative casing stones. Arnold (1991) reported that the density of limestone used for the casing stones ranges from 2.65 to 2.85 gcm -3 . Tura limestone is thought to have  covered the outside layer of the Pyramid and accounted for about 5% of its entire volume.
In comparison to the Tura limestone, the Mokattam limestone is more porous and less dense, and was used for the core of the pyramid. Although the blocks of Mokattam limestone were exploited near the Chephren Pyramid, some of them were transported from a quarry located in southeastern Cairo, Egypt. Geology of this quarry is attributed to the Middle Eocene series. According to Arnold (1991), the average density of limestone used for building the core of the pyramid ranges from 1.7 to 2.6 gcm -3 . Gerald Lynch did a direct measurement of the rock piece exploited from the surface of the Chephron Pyramid in 1968, and derived a density of 1.8 gcm -3 (Alvarez, 1987).   Alvarez et al. (1970) collected 650,000 muons during their muography experiment in 1968. 100,000 of these muons had passed through the upper zone that consists of casing, back stone, and core layers, as defined in Fig. 2. However, there is no description about the measurement time, t in Alvarez et al. (1970). This unknown parameter is cancelled if we take the ratio of the number of muon events accumulated in the different histogram bins. The muon counts are statistically sufficient to apply Eq. (7) to compare the subsurface density (the density of the deviation near the apex that consists of casing and back stones) with the Pyramid's core density (the average density of the Pyramid without this deviation part). For the calculation of  0 and 1  , the virtual detector was located 13.5 m east and 4 m north of the center on the ground level as described in the report written by Alvarez et al. (1970) (Fig. 2) In Table 1, we show the number of muon events (N) that Alvarez et al. (1970)  In Table 2 and Fig. 3, the ratio r0/r1 is shown as a function of elevation angle. The ratio does not vary with the muon's arriving angles beyond the error bars associated to these

Weight of the Pyramid
As an outlet of deriving the density ratio between the subsurface density and the core density of the Pyramid, we attempted to derive the total weight of the Pyramid. Lehner Therefore, we assumed that a quarter of the volume of the deviation consists of Tura limestone within a density range of 2.75±0.10 gcm -3 . If we assume the back stone density is 1.8 g/cm 3 , the Pyramid's core density will be 1.89±0.20 gcm -3 that is in agreement with the back stone density within the statistical error. Since both casing stones and back stones are tightly packed without spaces (at least near the surface) this assumption is reasonable.
By considering the total volume of the Chephren Pyramid to be 2,211,096 m 3 , the total weight of the Pyramid is derived to be 3.98 × 10 6 tons.

Case study 2: Grotta Gigante
The area around the natural limestone cave called "Grotta Gigante" consists of 400 m thick limestone bedrock including a known aquifer located 250 m below the ground surface, and a number of branch caves are undiscovered. Karst topography is a unique geological landscape formed by the natural dissolution of soluble bedrock, mostly made up of a thick layer of limestone that is eroded by being dissolved gradually into the underground aquifer. This dissolution is likely to occur in the crack of the ground, and thus, the specific part of the ground tends to be eroded and forms a doline (sink hole) on the ground surface and a limestone cave underground. As a result, a complicated cave system that consists of a number of small caves will spread around the region where the large main cave is developed. Many of the caves are not opened to the outside, and therefore it is expected that the unique ecosystem is developed there without being disturbed by exterior biological activities (e.g. Rohwerder et al., 2003). Caffau et al. (1997)  (7), X /X0 is then calculated to be 2.6 and here we assumed that the vertical muon flux closely matches with the flux in the direction of 50 o from zenith (Haino et al. 2004;Achard et al. 2004). Consequently, the areal density of the overburden in the direction of 50 o is derived to be 140 m.w.e. If we assumed the uniform density of 2.7 gcm -3 , the overburden thickness in this direction will be 52 m that is slightly shorter than the expected 1  .
Once the muographic anomalies are detected, we can compare them with the gravimetric data. In the area of the Grotta Gigante, the Osservatorio Geofisico Sperimentale (OGS) conducted a gravimetric survey using a LaCoste-Romberg microgravity meter. More than 200 data points were acquired in an area of ~800 x 650 m 2 to map out the gravity anomaly in this area. The region where the muographic anomaly was observed showed a gravity deviation of 0.1 mgal from the value expected from the terrain topography. By combining these data with the gravimetric data, it was revealed that this anomaly came from a red-soil deposit laying beneath the doline, and its volume and total weight were estimated to be 5 x 10 3 m 3 and 8.5 x 10 6 kg, respectively (Caffau et al. 1997).

Case study 3: Ore body explorations
This muographic data analysis may be applied also to localization of ore bodies because typically the ore density is 1.4-1.8 times higher than the surrounding medium. Liu et al. (2012) conducted mugraphic observations in the Price 5 deposit at the Myra Falls mine, Canada in order to estimate the total weight of the zinc ore body. This target was suitable for the proof of concept trial for muographic mineral explorations because the mine gallery exists at relatively shallow depths from the ground surface and a density map is available based on the diamond drilling results. The averaged ore and the surrounding medium density were respectively measured to be 3.2 gcm -3 and 2.7 gcm -3 based on the drilling and sampling method prior to their muographic measurements. A  (7), (X 0/X1 ) -1 is then calculated to be 1.6. Since the areal density of Overburden A is calculated to be ~270 m.w.e., the areal density of Overburden B is therefore derived to be ~430 m.w.e., therefore the average density is 3.1 gcm -3 . By inputting the ore density of 3.2 gcm -3 , the thickness of the ore can be calculated to be ~100 m which is in agreement with the result from the prior drilling and sampling work.
We anticipate that this method is also applicable to exploring the pyrite-polymetallic and the wolfram deposits of the Greater Caucasus (Eppelbaum and Khesin, 2012) This technique is applicable to future underground muography observations by utilizing, e.g., a borehole (Fig 4). Once we determine the near surface density (r0) with drilling and sampling methods, the average density of the depth region (d1) between the core sampling depth (d0) and the detector location (d) will be measured without requiring knowledge of the active area (S), detection efficiency (Aeff), and solid angle () of the detector. If the measurement time is fixed for each measurement, the ratio of the number of muon events (N0 / N1) counted at depths of Depth 0 and Depth 1 will give the density ratio r0/r1 by considering the muon's path length in each layer, which is proportional to the depths, d0 and d1. Therefore, if we compare the number of muon events measured at different depths, the average density above the detector locations are derived one after another, and thus the vertical density distribution of the soil will be obtained. If two or more boreholes are available for this kind of the measurement, three dimensional information on the density distribution will be obtained.

Conclusion
In this work, we evaluated the relationship between the transmitted muon flux and the areal density along the muon path, and found that it has a near-linear relationship as long as the overburden thickness is thinner than a few hundred meters. Based on this finding, we proposed a simple analysis method to cancel the contributions from the active area (S), detection efficiency (Aeff), and solid angle () of the detector in the generated muogram by combining the independently measured density information for the partial  volume of the target.
We showed two examples as possible applications of this analysis method. By reanalyzing the muographic data collected by Alvarez et al. (1970) and combining the surface sampling results with them, we derived the bulk density of the core of the Pyramid of Chephren, Egypt; hence its total weight. By combining the drilling and sampling results taken at the region near Grotta Gigante, Italy with the muographic data collected by Caffau et al. (1997), we calculated the size of the void associated with the doline caused by karst process. The derived size with our method was consistent with the result obtained with the conventional method, which utilized the absolute value of the transmitted muon flux.
Combining non-muographic densimetric techniques with muography has been shown to be useful for improving the spatial resolution of the density image. We showed in the case study of Groitta Gigante, that this simple analysis method could be even more powerful by combining gravimetric measurements to provide useful geological information. Recently, Jourde et al. (2015) showed that gravimetric and muographic joint measurements enhance the resolving power of the technique, in particular, when the muographic measurement is unidirectional. Furthermore, since muographic and gravimetric measurements respectively derive horizontally and vertically integrated density, a resolution of the deeper region where muography is not applicable becomes greatly improved by this kind of joint measurement.
The near-surface density is relatively easier to measure, e.g., with a drilling sampling method, in comparison to deeper region densities. Therefore, we anticipate that this technique is useful for future muographic measurements of the rock overburden for the purpose of exploring geological structures that require an accurate density value; for example, surveys of natural resources and mechanical fracture region of the fault zones.