Monitoring temporal opacity fluctuations of large structures with muon tomography : a calibration experiment using a water tower tank
Kevin Jourde, Dominique Gibert, Jacques Marteau, Jean de Bremond d'Ars, Serge Gardien, Claude Girerd, Jean-Christophe Ianigro
MMonitoring temporal opacity fluctuations of large structures with muon radiography: acalibration experiment using a water tower
Kevin Jourde, Dominique Gibert,
2, 3
Jacques Marteau, Jean de Bremondd’Ars, Serge Gardien, Claude Girerd, and Jean-Christophe Ianigro Institut de Physique du Globe de Paris (CNRS UMR 7154), Sorbonne Paris Cit´e, Paris, France. OSUR - G´eosciences Rennes (CNRS UMR 6118), Universit´e Rennes 1, Rennes, France. National Volcano Observatories Service, Institut de Physique du Globe de Paris (CNRS UMR 7154), Paris, France. Institut de Physique Nucl´eaire de Lyon, Univ Claude Bernard (UMR 5822 CNRS), Lyon, France.
Usage of secondary cosmic muons to image the geological structures density distribution signifi-cantly developed during the past ten years. Recent applications demonstrate the method interest tomonitor magma ascent and volcanic gas movements inside volcanoes. Muon radiography could beused to monitor density variations in aquifers and the critical zone in the near surface. However, thetime resolution achievable by muon radiography monitoring remains poorly studied. It is biased byfluctuation sources exterior to the target, and statistically affected by the limited number of particlesdetected during the experiment. The present study documents these two issues within a simple andwell constrained experimental context : a water tower.. We use the data to discuss the influence ofatmospheric variability that perturbs the signal, and propose correction formulas to extract the muonflux variations related to the water level changes. Statistical developments establish the feasibilitydomain of muon radiography monitoring as a function of target thickness (i.e. opacity). Objects witha thickness comprised between ≈ ±
30 m water equivalent correspond to the best time resolution.Thinner objects have a degraded time resolution that strongly depends on the zenith angle, whereasthicker objects (like volcanoes) time resolution does not.
Keywords: Cosmic muons, Tomography, Monitoring
INTRODUCTION
Using the secondary cosmic rays muon component to image geological bodies like volcano lava domes is thesubject of increasing interest over the past ten years. Much like medical X-ray radiography, muon radiography aimsat recovering the density distribution, ρ , inside the targets by measuring their screening effect on the cosmic muonsnatural flux. This approach was first tested by George[1] to measure the thickness of the geological overburden of atunnel in Australia, and later by Alvarez et al.[2] who imaged the Egyptian Pyramid of Chephren to eventually finda hidden chamber. The method then stayed long dormant until recent years when, thanks to progress in electronicsand particle detectors, field instruments were designed and constructed by several research teams worldwide [3, 4].Muon radiography experiments have successfully been performed on volcanoes where the hard muon componentis able to cross several kilometres of rock [3, 5–12]. Applications to archaeology [13], civil engineering (tunnels,dams) and environmental studies (near surface geophysics) are subject to active research, and monitoring of densitychanges in the near surface constitutes an important objective in hydrology and soil sciences.The material property that can be recovered with muon radiography is the opacity, (cid:36) , which quantifies the amountof matter encountered by the muons along their travel path, L , across the volume to image, (cid:36) = (cid:90) L ρ ( l ) × d l . (0.1)Generally, the opacity is expressed in [ g.cm − ] or, equivalently, in centimetres water equivalent [ cm.w.e. ] . Muonslose their energy through matter by ionisation processes [14] at a typical rate of 2.5 MeV per opacity increment of1 g.cm − . They are relativistic leptons produced in the upper atmosphere at an altitude of about 16 km [15, 16], andreach the ground after losing about 2.5 GeV to cross the opacity of 10 m.w.e. represented by the atmosphere. Muonstravel along straight trajectories across low-density materials, including water, concrete and rocks, and scatteringis significant only in high-density materials like lead and uranium [14]. However, low-energy muons ( E ≤ a r X i v : . [ phy s i c s . i n s - d e t ] M a r − , i.e. several percent variations during perturbedmeteorological conditions. Consequently, monitoring subtle density changes in low-opacity targets necessitates aprecise correction of the muon flux time-variations induced by both atmospheric pressure variations and eventualintense geomagnetic events.The present study aims at contributing to the procedures development to monitor density changes in low-opacitybodies. We apply and discuss a simple way to suppress atmospheric pressure effects from muon counting data.We present a controlled experiment performed on a water tank tower whose opacity fluctuates in a significantrange (3 m.w.e. < (cid:36) < THE SHADOW EXPERIMENT
FIG. 1. Sketch of the SHADOW experiment. The three yellow rectangles are the detection matrices (each with 16 ×
16 pixels of 5 × ), the red dotted lines encompass the detection solid angle and the blue surface represents the watervolume. The SHADOW experiment measured the muonic component time-variations while the water level varied in awater tower. For this purpose, we placed a muon telescope (its description is given in the Methods Section below)along the tower symmetry axis and below the tank. We oriented the instrument vertically (i.e. central zenith angle= 0) as shown in Figure 1 so that the apparent opacity is only zenith angle (azimuthal invariance) and time (whenthe water level h ( t ) is changing) dependent. The water tower is located in Tignieu-Jameyzieu, France, a villagelocated 20 kilometres East from Lyon (altitude 230 m above sea level, X UTM =
31 669490, Y UTM = ◦ ≤ θ ≤ ◦ such that all the telescope 961 lines of sight pass through the water. The solid angle spanned by the telescope equals Ω int = T int =
630 cm sr.The data acquisition started on November 21 th , 2014 and stopped on January 22 nd , 2015. While measuring themuon flux under the tank, the water level was monitored with a several cm accuracy every 5 minutes by thecompany in charge of the tower (Syndicat Intercommunal des Eaux de Pont-de-Ch´eruy – SIEPC). These data arecompleted by atmospheric pressure hourly measurements at the nearby Saint Exupery airport located 7.45 km Westof the water tower at an altitude of 248 m above sea level. FIG. 2. Hourly averages of the data acquired during the calibration period (November 22 nd to December 13 th , 2014): (a) Atmospheric pressure time-variations relatively to p = (b) Muon flux relative raw time-variations (black curve),and corrected from the atmospheric pressure influence (red curve). (c)
Water level variations.
During the measurement period, the geomagnetic activity was monitored using the Kp geomagnetic index pub-lished by the International Service of Geomagnetic Indices (ISGI) of the International Association of Geomagnetismand Aeronomy (IAGA). From November 21 st , 2014 to January 22 nd , 2015, a noticeable geomagnetic activity is re-ported for December 2014 7 th , 12 th , 22 nd , 26 th , 29 th and January 2015 4 th , 5 th , 7 th , 8 th where the geomagnetic activityreached the ”minor” G1 level[23] excepted in polar regions where the G3 level was observed. Sudden storm com-mencements (SSC) are reported on December 21 st (19:11 UTC), 22 nd (15:11 UTC), 23 rd (11:15 UTC), and January 7 th (06:14 UTC). It cannot be excluded that the geomagnetic activity at these dates produced small variations, at thefraction of percent level [24], of the muon flux measured during the SHADOW experiment. However, these varia-tions are expected to occur only a few times in the data time-series and this sparsity prevents a detailed quantitativestudy to identify the corresponding signals.In the next section, we use hourly averages of these data series to document the relationship between the muonflux time-variations and those of both the atmospheric pressure and the water level in the tank. Constant water level: Atmospheric effects contribution
We first consider the data acquired during the measurement period first three weeks, from November 22 nd toDecember 13 th h =
496 cm (Figure 2c). Meanwhile, the atmospheric pressure varied by ±
15 hPa with respect to a referencepressure p = ≤ h ≤
496 cm are retained. A least-squares fit to these points gives a negative slope β p = − ( ) hPa − where the value in parenthesis is the half-width of the 95% confidence interval. Weperformed the fit by assigning to the relative flux averages a standard deviation σ Φ = σ p = σ r = σ Φ and indicates that no higher-order fit is required.Consequently, in the remaining, we shall represent the atmospheric influence on the relative muon flux with a linearrelationship, Φ − Φ Φ = β p × ( p − p ) . (0.2) FIG. 3. Muon flux relative variations versus atmospheric pressure deviation. The red line represents the best least-squares fitsolution. Only the data points corresponding to a water level greater than 495 cm have been kept to compute the line fit.
Dayananda [25] uses the same kind of linear relation and finds β p = − ( ) hPa − from muon countsat the Earth’s surface. Other authors [26–28] also find linear relationships with coefficients falling near β p = − − . We do not expect the barometric coefficient derived in the present study to be strictly equal to thoseobtained for other experiments since this coefficient is sensitive to the site location and especially to the telescopealtitude[29, 30]. However they should be in the same order of magnitude.The correction formula (0.2) says that a ∆ p =
10 hPa increase of the atmospheric pressure induces a relativemuon flux decreases of 1.2%. Applying this correction to the muon flux data (black curve of Figure 2b) efficientlyreduces the long-period variations (light red curve of Figure 2b).
Time-varying water level
We now consider the data from the second measurement period, where we observe large water level variations(Figure 4). It begins on December 13 th nd h (Figure 4c). During the same period, the muon flux variations appearclearly anti-correlated with the water level (Figure 4b), and the highest relative flux deviation reaches 15% whenthe water level is minimum ( ≈
320 cm). Meanwhile, the atmospheric pressure variations (Figure 4a) also produceconspicuous effects on the muon flux like, for instance, the flux bump that occurs around December 28 th during alow-pressure event.The circles in Figure 5 represent the muon flux data versus the water level. Applying the atmospheric correction(0.2) to the muon flux reduces the data points scattering and enhances the correlation between the flux and thewater level (black dots in Figure 5). The uncorrected points standard deviation σ = σ = σ = σ = β h = − ( ) cm − and an intercept ∆Φ w = ( ) . The standard deviations assigned to the water leveland muon flux are respectively σ h = σ Φ = FIG. 4. Hourly averages of the data acquired during the non-stationary period (December 13 th nd (a) Atmospheric pressure time-variations relatively to p = (b) Muon flux relative raw time-variations (black curve),and corrected from the atmospheric pressure influence (red curve). (c)
Water level variations.FIG. 5. Normalized and centered muon flux as a function of water level. The open circles correspond to the measured muonflux (i.e. black curve in Figure 4b) and the black dots correspond to the muon flux corrected from the atmospheric pressureinfluence (equation 0.2). The red straight line is a linear fit to the pressure-corrected points. muon flux, Φ − Φ Φ = β h × h + ∆Φ w . (0.3) DISCUSSION OF SHADOW DATA ANALYSIS
The data analyzed in the previous Section show that linear relationships (equations 0.2 and 0.3) may safely beused to represent the relative muon flux dependence with respect to the atmospheric pressure (Figure 3) and waterlevel (Figure 5) variations. Owing to the fact that the opacity fluctuations produced by ∆ p = ∆ h = β coefficients in equations(0.2) and (0.3) should be the same. This hypothesis is not supported by our experimental results which indicate that β p is significantly larger than β h .We explain the discrepancy between the experimental values for β p and β h by the fact that the atmosphere is notonly, like water in the tank, a screen of matter for the muon flux but it is also the place where muons originate[15, 16, 31]. Consequently, the muon flux at ground level depends on both the pressure and the temperatureprofiles in the atmosphere. For instance, if the atmosphere is warmer, the muon production altitude is higher(roughly at the isobaric level p =
100 hPa) and the muons transit times increase. Then, muons are more likely todecay before reaching the ground and thus the relative muon flux decreases [15, 31]. This is the so-called negativetemperature effect. However, an increase in temperature at the production level decreases the air density, thusreducing the likelihood of pion interactions before their decay into muons. Muon production then increases, andthis phenomena is known as the positive temperature effect. Both the pressure and temperature effect upon theflux of muons at ground level may be summarized by [32, 33], Φ − Φ Φ = β ∗ p × ( p − p ) + β ∗ T × ( T − T ) , (0.4)where T is the temperature at the production level and β ∗ p and β ∗ T are adjustable coefficients for the pressureand temperature effects respectively. The coefficient β ∗ p is always negative while β ∗ T may be either positive ornegative depending on the prevailing temperature effect. For the soft muon component ( ≤
10 GeV) which composesthe main part of the particles detected by our telescope in this experimental context, the negative temperatureeffect dominates and β ∗ T is expected to be negative. The correlation analysis recently performed by Zazyan etal. [34] shows that the pressure and temperature effects are positively correlated. Consequently, for the presentmeasurement conditions, both β ∗ p and β ∗ T are negative and time-correlated. When considering the atmosphericeffect alone like in equation (0.2), the β p coefficient actually accounts for both the pressure and temperature effects.This explains why the experimental value found for β p (0.2) is larger than the value of β h (0.3). STATISTICAL FEASIBILITY AND LIMITS OF OPACITY MONITORING
We now address some statistical issues concerning the monitoring of opacity variations like those produced bywater-level variations measured during the SHADOW experiment.Let us assume that N = N + N particles are detected by the telescope during a time period T , and where N and N are the number of particles respectively counted during the first and second half of T . We want to determineunder which conditions N and N may be considered different at the confidence level α . The particle flux difference ∆ N = N − N obeys a Skellam distribution defined as the difference between two Poisson processes with means µ and µ [35], S ( ∆ N , µ , µ ) = e − ( µ + µ ) (cid:18) µ µ (cid:19) k /2 I ∆ N ( √ µ µ ) , (0.5)where I ∆ N is the modified Bessel function of the first kind.In the case where N > N , the hypothesis ∆ N (cid:54) = α if − ∑ i = − ∞ S ( i , N , N ) + × S ( N , N ) ≤ − α (0.6) N = T /2 × φ = T /2 × φ × ( − ε /2 ) (0.7) N = T /2 × φ = T /2 × φ × ( + ε /2 ) . (0.8)with ε the flux variation percentage.When the inequality (0.6) becomes an equality we get T = T min , the minimum acquisition time necessary toresolve a flux difference given by the following set of parameters ( φ , ε , α ) . When ε is fixed, T min is the best timeresolution achievable to observe temporal relative flux variations larger than ε . When T min is fixed, we derive thebest relative flux variation, ε , detectable on a time-scale larger than T min . FIG. 6. Minimum acquisition time T min versus the average measured flux φ necessary to detect a ε flux variation with α = α curves respectively computed with equation (0.6)and the approximation (0.10) . The curved arrows delimit the resolution domains for the SHADOW experiment and typicalvolcano applications. The horizontal limit marked by the arrows is the measurement whole duration and the vertical limit isthe maximum flux measured. The crosses represent likely sources of muon flux variations, their coordinates depend on the fluxfluctuations amplitude α and their typical period T min . Note that if ( N , N ) (cid:38)
10 the Poisson laws can be approximated with Gaussians and equation (0.6) is simplifiedto, ( N − N ) − ˜ α × (cid:115) N × N N + N ≥ T ≥ T min = ˜ α (cid:0) − ε /4 (cid:1) ε × φ (0.10)where α = erf ( ˜ α ) .We numerically compute T min from equation (0.6) with a confidence level α = ε =
1, 0.1, 0.01, and 0.001 (i.e. 100%, 10%, 1%, and 0.1%).Observe that the approximation (0.10) is suitable for our range of applications, T min will be underestimated startingfrom ε (cid:38) N ≈ ε = φ > − must be measured. This solution is represented by the black crosslabelled “water tank” on Figure 6. The lower-left domain delimited by the curved black arrow in Figure 6 representsthe solution-domain for time scales and opacity variations of the SHADOW experiment category. This solution-domain is the region where flux variations can be resolved at a high confidence level. The arrow horizontal branchis limited by the experiment duration, and the vertical branch is placed at a level corresponding to the maximumflux that can be measured by the telescope. This latter quantity increases with the telescope acceptance, e.g. byincreasing the angular aperture (i.e. by reducing the distance between the detection matrices), or by groupingseveral lines of sight, or by using instruments with a larger detection surface.The feasibility domain for a typical volcano experiment is also represented on Figure 6 and delimited by the bluecurved arrow. Note that for this kind of experiments we have a longer acquisition time and a tiny measured fluxas the total opacity of the geological body facing the telescope is much bigger than for the SHADOW experiment:about 1000 m.w.e. for a volcanic lava dome versus 5 m.w.e. for the water tank.We can rewrite equation (0.6) into a form more suitable for radiography applications by replacing the flux fluctu-ations by opacity fluctuations, N = T × φ ( T , (cid:36) , θ ) (0.11) N = T /2 × φ ( T , (cid:36) × ( − ε /2 ) , θ ) (0.12) N = T /2 × φ ( T , (cid:36) × ( + ε /2 ) , θ ) , (0.13)where the muon flux φ , is explicitly written to depend on telescope acceptance T , opacity (cid:36) , and zenith angle θ .As before, ε represents the variation of opacity relative to the average opacity (cid:36) . We warn the reader that a given ε -variation of (cid:36) corresponds to a much larger ε -variation of φ . Putting equations (0.12) and (0.13) in equation (0.9)we obtain the following feasibility condition, T ≥ T min ( (cid:36) , ε , θ , T ) = × ˜ α × φ × φ ( φ − φ ) × ( φ + φ ) (0.14)where T min is the measurement period minimum duration necessary to resolve the sought opacity variation. Notethat the feasibility formula from Lesparre et al. [7] is the first order development of equation (0.14).A subset of T min solutions of equation (0.14) is represented on Figure 7 for the confidence level α = θ = ◦ , 30 ◦ , and 60 ◦ , and opacity variations ε = T =
10 cm .sr,typical of our telescopes has been used in the computation. Roughly, a one order of magnitude ε variation inducesa T min change by two orders of magnitude.Observe that there is an optimal opacity range where the measurement time, i.e. the time resolution that isachievable, is minimum to resolve a given opacity variation. The optimal opacity range depends on the zenith angleand goes roughly from 40 to 100 m.w.e for θ = ◦ , and from 20 to 40 m.w.e for θ = ◦ . For low-opacity conditions,measurements at high zenith angles are wise to optimize the time resolution. This is particularly conspicuous forthe SHADOW experiment where the average opacity (cid:36) ≈ ε ≈ T min > θ = ◦ to resolve the fluctuations while T min > θ = ◦ . Thetime resolution strong dependence with respect to the zenith angle disappears at larger opacities (cid:36) >
500 m.w.elike those encountered in volcano muon radiography.
DISCUSSION
Muon radiography is a powerful method to monitor opacity/density variations inside geological bodies. Notice-able advantages of the method are the possibility to remotely radiography unapproachable dangerous volcanoesand to image the density distribution of large volumes from a single view-point [8, 10, 36]. Muon radiography isentering an era of precision measurements not only for structural imaging but also for dynamical monitoring pur-poses. Some monitoring experiments have been performed on active volcanoes that demonstrate the usefulness ofsuch measurements to constrain the evolution of eruption crisis [12]. However, as shown above, monitoring opacityvariations is subject to external sources of bias, and statistical and experimental constraints that limit the achievableresolution. Understanding these limits is of primary importance to improve the method and to assess the muonradiography monitoring feasibility and validity.Experimental constraints are partly dictated by statistical considerations, and mainly come from the telescopeacceptance that limits the maximum flux which fixes the resolution domain right boundary in Figure 6. Thisboundary may be moved rightward by increasing the acceptance T of the instrument. Recalling that T is expressedin [ cm sr ] , the acceptance may be augmented by several means: 1) increasing the solid angle encompassed by theinstrument by reducing the distance between the detection matrices; 2) increasing the detection surface by couplingseveral telescopes (actually our telescopes may be merged into a single one); 3) grouping lines of sight to increaseboth the detection surface and the solid angle at the price of reducing the angular resolution of the radiographies. Inthe present study, the latter solution was retained and all lines of sight were merged to obtain an effective acceptanceof 630 cm sr.Statistical constraints bound the resolution domain of a given experiment (Figure 6), and the main concern whendoing measurements is to ensure that the monitored phenomena fall inside the boundaries. As will be discussed inthe next paragraph, the telescope configuration may be adapted to comply with the ongoing experiment objectives.As shown in the preceding sections, the statistical constraints are quite different whether the opacity is high or low.This is conspicuous in Figure (7) where the feasibility solutions for T min strongly differ in the low- and high-opacitydomains. It is remarkable that high-opacities variations are equally resolved whatever the zenith angle while, FIG. 7. Minimum acquisition time T min as a function of the average opacity (cid:36) to detect an ε fluctuation at the α = ◦ , 30 ◦ ,60 ◦ ) and are computed for T =
10 cm .sr using the modified Gaisser model from Tang et al (2006). instead, the resolution for low-opacities strongly depends on this angle. Another conspicuous feature present inFigure (7) is the existence of an optimal medium-opacity range (cid:36) ≈ ±
30 m.w.e. where telescopes offer their bestperformance. These two effects are due to the cosmic muon energy spectrum nature [15, 17, 37], and changing thetelescope acceptance has no effect on the optimal opacity values but only changes T min by translating the solutioncurves of Figure (7) either upward (decrease of acceptance) or downward (increase of acceptance). METHODS
The muon count series analysed in the present study were acquired with one of our standard telescopes shownin Figure 8 [4, 10, 38]. The picture was taken during an open-sky calibration phase where the muon count servesto determine the efficiency of the scintillator bars forming the detection matrices. Each matrix is formed by anassemblage of two sets of 16 bars arranged perpendicularly to obtain a 16 ×
16 square 5 × pixels array. Thetelescope upper and lower matrices allow 31 ×
31 pixels combinations, i.e. 961 distinct lines of sight. The distancebetween the matrices may be changed to adapt the solid angle spanned by the trajectories. In the present study, thedistance was tuned to encompass the entire water tank (Figure 1).Once geometrically configured, the telescope is totally characterised by its acceptance function T i [ cm sr ] whichrelates the muon count, N i , to the muon flux, ∂φ [ s − cm − sr − ] received by the telescope in its i th line of sight, N i = T × (cid:90) π P i ( ϕ , θ ) × ∂φ ( (cid:36) , ϕ , θ ) × d Ω , (0.15) = T × T i × ∂φ i , (0.16)where T is the acquisition duration, P i [ cm ] is the line of sight detection surface function, T i is the integrated ac-ceptance, and ∂φ i is the muon flux in the line of sight central direction. It must be understood that ∂φ [ s − cm − sr − ] is the differential muon flux that reaches the instrument after crossing the target. Consequently, ∂φ depends bothon the open sky differential flux ∂φ ( (cid:36) = ϕ , θ ) and on the muon absorption law inside matter. These are deter-mined through experiments [26, 27, 37, 39–42], theoretical works [17] or thanks to Monte-Carlo simulations [43, 44]depending on the precision expected and the available information.Figure (9) shows the telescope acceptances T i for i = · · · , 961 used in the SHADOW experiment. This accep-tance function is determined experimentally to account for the detection matrices deffects, mainly imperfect optical0 FIG. 8. Picture of the muon telescope used for the SHADOW experiment, here during the open-sky calibration phase. Thethree detection matrices are horizontal. The calibration gives access to the effective acceptance. The control box embedding amini-PC, a common clock distribution system, a network switch is visible on the middle matrix.FIG. 9. Telescope experimental acceptance for the configuration shown in Figure (1). The acceptance maximum value, T max = sr is obtained for the line of sight perpendicular to the detector planes and corresponding to ( x , y ) = (
0, 0 ) . The x and y coordinates represent the horizontal offsets between the pixels defining a given line of sight of the telescope (one pixel in theupper detection matrix, and the other one in the lower matrix). The acceptance integrated over the instrument entire detectionsurface equals T int =
630 cm sr for a solid angle aperture Ω int = couplings at the scintillator bars outputs and on the multichannel photomultiplier front. The latter causes the dis-tortions visible in the Figure (9) 3D plot. In practice, the acceptance computation is performed by measuring the“open-sky” muons flux coming from the zenith.The detected particles number N may be increased by grouping several adjacent lines of sight belonging to a1subset E , N E = T × ∑ i ∈E T i × ∂φ i = ∑ i ∈E N i . (0.17)It results in an acceptance increase and thus a better time resolution. The counterpart is an angular resolutiondegradation induced by the merging of the small solid angles spanned by the trajectories. In the present study, theentire solid angle spanned by the telescope trajectories were grouped to obtain a total acceptance T total =
630 cm sr.Such a large acceptance dramatically improves the time resolution which falls to the order of tens of minutes in thecase of the SHADOW experiment. ACKNOWLEDGEMENTS
This study is part of the DIAPHANE project ANR-14-CE 04-0001. We acknowledge the financial support from theUnivEarthS Labex program of Sorbonne Paris Cit´e ( anr -10- labx -0023 and anr -11- idex -0005-02). We would like tothank also the members of the SIEPC, Tignieu-Jameyzieu, France, for their help, support and the access to the data.In particular we thank Mr Gilbert Pommet and Mr R´emi Cachet. The results presented in this paper rely on rely ongeomagnetic Kp indices calculated and made available by ISGI Collaborating Institutes and on the data collected atChambon-la-Foret magnetic observatory. We thank the Institut de Physique du Globe de Paris, for supporting itsoperation and the INTERMAGNET network and ISGI (isgi.unistra.fr). This is IPGP contribution ****.
AUTHOR CONTRIBUTIONS STATEMENT
K.J., J.M. and D.G. conceived the experiment, all authors designed and constructed the apparatus, J.M. conductedthe experiment, K.J., J.M. and D.G. analysed the data. All authors reviewed the manuscript.
ADDITIONAL INFORMATION
The authors declare no competing financial interests.
The corresponding author is responsible for submitting a competing financial interests statement on behalf of allauthors of the paper. This statement must be included in the submitted article file. [1] George, E. P. Cosmic rays measure overburden of tunnel.
Commonwealth Engineer et al.
Search for hidden chambers in the pyramids.
Science , 832–839 (1970).[3] Tanaka, H. et al.
Development of the cosmic-ray muon detection system for probing internal-structure of a volcano.
HyperfineInteract. , 521–526 (DOI:10.1023/A:1020843100008) (2001).[4] Marteau, J. et al.
Implementation of sub-nanosecond time-to-digital convertor in field-programmable gate array: applicationsto time-of-flight analysis in muon radiography.
Meas. Sci. Technol. , 035101 (DOI:10.1088/0957-0233/25/3/035101) (2014).[5] Nagamine, K. Geo-tomographic observation of inner-structure of volcano with cosmic-ray muons. J. Geogr. , 998–1007(1995).[6] Tanaka, H., Nagamine, K., Nakamura, S. & Ishida, K. Radiographic measurements of the internal structure of mt.west iwate with near-horizontal cosmic-ray muons and future developments.
Nucl. Instrum. Methods
A555 , 164–172(DOI:10.1016/j.nima.2005.08.099) (2005).[7] Lesparre, N. et al.
Geophysical muon imaging: feasibility and limits.
Geophys. J. Int. , 1348–1361 (DOI:10.1111/j.1365-246X.2010.04790.x) (2010).[8] Lesparre, N. et al.
Density muon radiography of la soufriere of guadeloupe volcano: comparison with geological, electricalresistivity and gravity data.
Geophys. J. Int. , 1008–1019 (DOI:10.1111/j.1365-246X.2012.05546.x) (2012).[9] Marteau, J. et al.
Muons tomography applied to geosciences and volcanology.
Nucl. Instrum. Methods
A695 , 23–28(DOI:10.1016/j.nima.2011.11.061) (2012).[10] Jourde, K. et al.
Experimental detection of upward going cosmic particles and consequences for correction of densityradiography of volcanoes.
Geophys. Res. Lett. , 6334–6339 (DOI:10.1002/2013GL058357) (2013).[11] Portal, A. et al. Inner structure of the puy de dˆome volcano: cross-comparison of geophysical models (ert, gravimetry, muonimaging).
Geosci. Instrum. Method Data Syst. , 47–54 (DOI:10.5194/gi-2-47-2013) (2013).[12] Tanaka, H. K., Kusagaya, T. & Shinohara, H. Radiographic visualization of magma dynamics in an erupting volcano. Nature (DOI:10.1038/ncomms4381) (2014). [13] Menichelli, M. et al. A scintillating fibres tracker detector for archaeological applications.
Nuclear Instrumentsand Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment , 262–265(DOI:10.1016/j.nima.2006.10.317) (2007).[14] Nagamine, K.
Introductory muon science (Cambridge University Press, 2003).[15] Gaisser, T. K.
Cosmic rays and particle physics (Cambridge University Press, 1990).[16] Grieder, P. K.
Cosmic rays at Earth (Gulf Professional Publishing, 2001).[17] Tang, A., Horton-Smith, G., Kudryavtsev, V. A. & Tonazzo, A. Muon simulations for super-kamiokande, kamland, andchooz.
Phys. Rev.
D74 , 053007 (DOI:10.1103/PhysRevD.74.053007) (2006).[18] Yanchukovsky, V., Filimonov, G. Y. & Hisamov, R. Atmospheric variations in muon intensity for different zenith angles.
Bull. Russ. Acad. Sci. Phys. , 1038–1040 (DOI:10.3103/S106287380707043X) (2007).[19] Kremer, J. et al. Measurements of ground-level muons at two geomagnetic locations.
Phys. Rev. Lett. , 4241(DOI:10.1103/PhysRevLett.83.4241) (1999).[20] Munakata, K. et al. Precursors of geomagnetic storms observed by the muon detector network.
J. Geophys. Res. , 27457–27468 (DOI:10.1029/2000JA000064) (2000).[21] Chilingarian, A. et al.
Correlated measurements of secondary cosmic ray fluxes by the aragats space-environmental centermonitors.
Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and AssociatedEquipment , 483–496 (DOI:10.1016/j.nima.2004.12.021) (2005).[22] Maghrabi, A., Al Harbi, H., Al-Mostafa, Z., Kordi, M. & Al-Shehri, S. The kacst muon detector and its application tocosmic-ray variations studies.
Advances in Space Research , 700–711 (DOI:10.1016/j.asr.2011.10.011) (2012).[23] Poppe, B. B. New scales help public, technicians understand space weather. Eos, Transactions American Geophysical Union ,322–328 (DOI:10.1029/00EO00247) (2000).[24] Ishibashi, S. et al. Observation of solar-terrestrial variation of the cosmic muon flux. In
Nuclear Science Symposium ConferenceRecord, 2005 IEEE , vol. 3, 1441–1446 (DOI:10.1109/NSSMIC.2005.1596591) (IEEE, 2005).[25] Dayananda, M., Zhang, X., Butler, C. & He, X. Understanding the effect of atmospheric density on the cosmic ray fluxvariations at the earth surface. arXiv:1303.7191 (2013).[26] Sagisaka, S. Atmospheric effects on cosmic-ray muon intensities at deep underground depths.
Nuovo Cimento C9 , 809–828(DOI:10.1007/BF02558081) (1986).[27] Motoki, M. et al. Precise measurements of atmospheric muon fluxes with the bess spectrometer.
Astropart. Phys. , 113–126(DOI:10.1016/S0927-6505(02)00195-0) (2003).[28] Rigozo, N. R. & Petry, A. The atmospheric pressure effect on muon data normalization by spectral analysis studies. RevistaBrasileira de Geof´ısica , 507–514 (DOI:10.1590/rbgf.v31i3.324) (2013).[29] Dorman, L. I. Cosmic Rays in the Earths Atmosphere and Underground (Springer Netherlands, 2004).[30] Chilingarian, A. & Karapetyan, T. Calculation of the barometric coefficients at the start of the 24th solar activity cycle for par-ticle detectors of aragats space environmental center.
Advances in Space Research , 1140–1146 (DOI:10.1016/j.asr.2010.12.001)(2011).[31] Anchordoqui, L. et al. High energy physics in the atmosphere: Phenomenology of cosmic ray air showers.
Ann. Phys. ,145–207 (DOI:10.1016/j.aop.2004.07.003) (2004).[32] Blackett, P. M. On the instability of the barytron and the temperature effect of cosmic rays.
Phys. Rev. , 973(DOI:10.1103/PhysRev.54.973) (1938).[33] Duperier, A. A new cosmic-ray recorder and the air-absorption and decay of particles. Terrestr. Magnet. Atm. Electr. , 1–7(DOI:10.1029/TE049i001p00001) (1944).[34] Zazyan, M., Ganeva, M., Berkova, M., Yanke, V. & Hippler, R. Atmospheric effect corrections of mustang data. J. SpaceWeather Space Clim. , A6 (DOI:10.1051/swsc/2015007) (2015).[35] Skellam, J. G. The frequency distribution of the difference between two poisson variates belonging to different populations. J. Roy. Stat. Soc.
A109 , 296–296 (1945).[36] Nagamine, K., Iwasaki, M., Shimomura, K. & Ishida, K. Method of probing inner-structure of geophysical substance withthe horizontal cosmic-ray muons and possible application to volcanic eruption prediction.
Nucl. Instrum. Methods
A356 ,585–595 (DOI:10.1016/0168-9002(94)01169-9) (1995).[37] Hebbeker, T. & Timmermans, C. A compilation of high energy atmospheric muon data at sea level.
Astropart. Phys. ,107–127 (DOI:10.1016/S0927-6505(01)00180-3) (2002).[38] Lesparre, N. et al. Design and operation of a field telescope for cosmic ray geophysical tomography.
Geosci. Instrum. MethodData Syst. , 33–42 (DOI:10.5194/gi-1-33-2012) (2012).[39] Ambrosio, M. et al. Seasonal variations in the underground muon intensity as seen by macro.
Astropart. Phys. , 109–124(DOI:10.1016/S0927-6505(97)00011-X) (1997).[40] Adamson, P. et al. Observation of muon intensity variations by season with the minos far detector.
Phys. Rev. D , 012001(DOI:10.1103/PhysRevD.81.012001) (2010).[41] Tilav, S. et al. Atmospheric variations as observed by icecube. arXiv:1001.0776 (2010).[42] Poirier, J. & Catanach, T. Periodic variations in muon flux at project grand. In
Proc. Int. Cosmic Ray Conf. , vol. 11, 173–176(2011).[43] Heck, D., Schatz, G., Knapp, J., Thouw, T. & Capdevielle, J. CORSIKA : A Monte Carlo Code to Simulate Extensive AirShower , FZKA Report 6019, Forschungszentrum Karlsruhe (1998).[44] Wentz, J. et al.
Simulation of atmospheric muon and neutrino fluxes with corsika.
Phys. Rev.