Monitoring temporal opacity fluctuations of large structures with muon radiography: a calibration experiment using a water tower

Usage of secondary cosmic muons to image the geological structures density distribution significantly developed during the past ten years. Recent applications demonstrate the method interest to monitor magma ascent and volcanic gas movements inside volcanoes. Muon radiography could be used to monitor density variations in aquifers and the critical zone in the near surface. However, the time resolution achievable by muon radiography monitoring remains poorly studied. It is biased by fluctuation sources exterior to the target, and statistically affected by the limited number of particles detected during the experiment. The present study documents these two issues within a simple and well constrained experimental context: a water tower. We use the data to discuss the influence of atmospheric variability that perturbs the signal, and propose correction formulas to extract the muon flux variations related to the water level changes. Statistical developments establish the feasibility domain of muon radiography monitoring as a function of target thickness (i.e. opacity). Objects with a thickness comprised between ≈50 ± 30 m water equivalent correspond to the best time resolution. Thinner objects have a degraded time resolution that strongly depends on the zenith angle, whereas thicker objects (like volcanoes) time resolution does not.

Usage of secondary cosmic muons to image the geological structures density distribution significantly developed during the past ten years. Recent applications demonstrate the method interest to monitor magma ascent and volcanic gas movements inside volcanoes. Muon radiography could be used to monitor density variations in aquifers and the critical zone in the near surface. However, the time resolution achievable by muon radiography monitoring remains poorly studied. It is biased by fluctuation sources exterior to the target, and statistically affected by the limited number of particles detected during the experiment. The present study documents these two issues within a simple and well constrained experimental context: a water tower. We use the data to discuss the influence of atmospheric variability that perturbs the signal, and propose correction formulas to extract the muon flux variations related to the water level changes. Statistical developments establish the feasibility domain of muon radiography monitoring as a function of target thickness (i.e. opacity). Objects with a thickness comprised between ≈50 ± 30 m water equivalent correspond to the best time resolution. Thinner objects have a degraded time resolution that strongly depends on the zenith angle, whereas thicker objects (like volcanoes) time resolution does not.
Using the secondary cosmic rays muon component to image geological bodies like volcano lava domes is the subject of increasing interest over the past ten years. Much like medical X-ray radiography, muon radiography aims at recovering the density distribution, ρ, inside the targets by measuring their screening effect on the cosmic muons natural flux. This approach was first tested by George 1 to measure the thickness of the geological overburden of a tunnel in Australia, and later by Alvarez et al. 2 who imaged the Egyptian Pyramid of Chephren to eventually find a hidden chamber. The method then stayed long dormant until recent years when, thanks to progress in electronics and 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 component is 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 density changes 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,  which quantifies the amount of matter encountered by the muons along their travel path, L, across the volume to image, Generally, the opacity is expressed in [g cm −2 ] or, equivalently, in centimetres water equivalent [cm.w.e.]. Muons lose their energy through matter by ionisation processes 14 at a typical rate of 2.5 MeV per opacity increment of 1 g cm −2 . They are relativistic leptons produced in the upper atmosphere at an altitude of about 16 km 15,16 , and reach the ground after losing about 2.5 GeV to cross the opacity of 10 m.w.e. represented by the atmosphere. Muons travel along straight trajectories across low-density materials, including water, concrete and rocks, and scattering is significant only in high-density materials like lead and uranium 14 . However, low-energy muons (E ≤ 1 GeV) have strong scattering in almost all materials.
Muon radiography of kilometre-size objects like volcanoes involves the hard muonic component with energy above several hundredths of GeV. In such cases, the muons incident flux may reasonably be considered stationary, azimuthally isotropic and to only depend on the zenith angle 7,15 , and simple flux models can be used to determine the screening effects produced by the target to image 7,17 . The situation is different in environmental and civil engineering applications where the bodies have low opacities (i.e. several tens of m.w.e.) that can be crossed by the soft muonic component (i.e. several GeV) which can no more be considered stationary and isotropic.
The soft muon component main causes of non-stationarity and anisotropic characteristics 15,16 are pressure variations at the ground level 18 , and geomagnetic storms together with solar coronal mass ejections (CME) 19,20 . Both time-variations of atmospheric pressure and CME involve time constants of several hours or days 21 which are of a critical importance when monitoring fast density changes in low-opacity bodies. During CME and associated magnetic storms, one generally observes Forbush decreases that correspond to a deficit of cosmic rays of 1% or 2% at ground level, but variations up to 10% have been reported 22 . The atmospheric pressure variations typically produce muon flux relative changes of 0.1% hPa −1 , i.e. several percent variations during perturbed meteorological conditions. Consequently, monitoring subtle density changes in low-opacity targets necessitates a precise correction of the muon flux time-variations induced by both atmospheric pressure variations and eventual intense geomagnetic events.
The present study aims at contributing to the procedures development to monitor density changes in low-opacity bodies. 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 significant range (3 m.w.e. <  < 5 m.w.e.) where atmospheric effects are expected to significantly perturb the incident cosmic muons flux. We made measurements during a several weeks period while the opacity remains steady at its maximum level before fluctuating. Meanwhile, water level in the tank, atmospheric pressure, and geomagnetic activity are monitored in order to evaluate their relative importance to produce muon flux variations across the water volume. Finally, a discussion about the time resolution in muon radiography monitoring is presented with a particular emphasis for the low-opacity targets case.

The SHADOW experiment
The SHADOW experiment measured the muonic component time-variations while the water level varied in a water 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 Fig. 1 so that the apparent opacity is only zenith angle (azimuthal invariance) and time (when the water level h(t) is changing) dependent. The water tower is located in Tignieu-Jameyzieu, France, a village located 20 kilometres East from Lyon (altitude 230 m above sea level, X UTM = 31 669490, Y UTM = 5067355). The distance between the upper and the lower matrices is set at 195 cm to cover a zenith angle range 0° ≤ θ ≤ 22.3° such that all the telescope 961 lines of sight pass through the water. The solid angle spanned by the telescope equals Ω int = 0.161 sr, and the total effective acceptance int  = 630 cm 2 sr. The data acquisition started on November 21 th , 2014 and stopped on January 22 nd , 2015. While measuring the muon flux under the tank, the water level was monitored with a several cm accuracy every 5 minutes by the company in charge of the tower (Syndicat Intercommunal des Eaux de Pont-de-Chéruy-SIEPC). These data are . It cannot be excluded that the geomagnetic activity at these dates produced small variations, at the fraction of percent level 24 , of the muon flux measured during the SHADOW experiment. However, these variations are expected to occur only a few times in the data time-series and this sparsity prevents a detailed quantitative study to identify the corresponding signals.
In the next section, we use hourly averages of these data series to document the relationship between the muon flux 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 to December 13 th 2014, when the water level in the water tower remained almost constant at its maximum level h 0 = 496 cm (Fig. 2c). Meanwhile, the atmospheric pressure varied by ± 15 hPa with respect to a reference pressure p 0 = 1016.8 hPa (Fig. 2a). The muon flux shown on Fig. (2b) not only randomly fluctuates as expected for a Poissonian process but also contains long-period variations with an amplitude of less than 3%. These long-period variations are clearly anti-correlated with those of the atmospheric pressure (Fig. 2a).
Since the water level is mainly constant during the considered period, we expect the muon flux time variations to be principally caused by atmospheric effects. The graph in Fig. 3 represents the muon flux hourly averages with respect to the atmospheric pressure from Saint Exupery airport. In this graph, only the data points such that the water level 495 cm ≤ h ≤ 496 cm are retained. A least-squares fit to these points gives a negative slope β p = − 0.0012 (0.0001) hPa −1 where the value in parenthesis is the half-width of the 95% confidence interval. We performed the fit by assigning to the relative flux averages a standard deviation σ Φ = 0.0081 derived from the events arrival times statistics. A standard deviation σ p = 1 hPa is assigned to the atmospheric pressure data. The linear fit residuals standard deviation, σ r = 0.0093, falls near σ Φ 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 linear relationship,  [26][27][28] also find linear relationships with coefficients falling near β p = − 0.001 hPa −1 . We do not expect the barometric coefficient derived in the present study to be strictly equal to those obtained for other experiments since this coefficient is sensitive to the site location and especially to the telescope altitude 29,30 . However they should be in the same order of magnitude. The correction formula (2) says that a Δp = 10 hPa increase of the atmospheric pressure induces a relative muon flux decreases of 1.2%. Applying this correction to the muon flux data (black curve of Fig. 2b) efficiently reduces the long-period variations (light red curve of Fig. 2b).

Time-varying water level
We now consider the data from the second measurement period, where we observe large water level variations (Fig. 4). It begins on December 13 th 2014 and ends on January 22 nd 2015. The largest decrease in the water level is up to nearly 200 cm with respect to h 0 (Fig. 4c). During the same period, the muon flux variations appear clearly anti-correlated with the water level (Fig. 4b), and the highest relative flux deviation reaches 15% when the water level is minimum (≈ 320 cm). Meanwhile, the atmospheric pressure variations (Fig. 4a) also produce conspicuous effects on the muon flux like, for instance, the flux bump that occurs around December 28 th during a low-pressure event.
The circles in Fig. 5 represent the muon flux data versus the water level. Applying the atmospheric correction (2) to the muon flux reduces the data points scattering and enhances the correlation between the flux and the  water level (black dots in Fig. 5). The uncorrected points standard deviation σ = 0.019 is reduced to σ = 0.011 for the pressure-corrected data. The uncorrected points standard deviation regularly increases from σ = 0.011 to σ = 0.021 when the water level increases from 3-5 meters while the pressure-corrected points standard deviation remains constant in the whole range of water levels. We explain this feature by the fact that high water levels are more frequent than low levels. Consequently the atmospheric pressure fluctuates in a wider range during the time period of high water level, causing larger muon flux variations.
A linear fit to the pressure-corrected points is displayed by the red line in Fig. 5 with a negative slope β h = − 0.0009 (0.00002) cm −1 and an intercept ΔΦ w = 0.444 (0.007). The standard deviations assigned to the water level and muon flux are respectively σ h = 1.7 cm and σ Φ = 0.0082. The residuals standard deviation is not reduced when fitting a second-order polynomial, and we adopt a linear relationship to represent the water level influence on the muon flux,

Discussion of SHADOW data analysis
The data analyzed in the previous Section show that linear relationships (equations 2 and 3) may safely be used to represent the relative muon flux dependence with respect to the atmospheric pressure (Fig. 3) and water level (Fig. 5) variations. Owing to the fact that the opacity fluctuations produced by Δp = 1 hPa and Δh = 1 cm are identical (i.e. they represent the same mass of matter), it may be deduced that the β coefficients in equations (2) and (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 not only, 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 temperature profiles 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 to decay before reaching the ground and thus the relative muon flux decreases 15,31 . This is the so-called negative temperature effect. However, an increase in temperature at the production level decreases the air density, thus reducing the likelihood of pion interactions before their decay into muons. Muon production then increases, and this phenomena is known as the positive temperature effect. Both the pressure and temperature effect upon the flux of muons at ground level may be summarized by 32,33 , where T is the temperature at the production level and β ⁎ p and β ⁎ T are adjustable coefficients for the pressure and temperature effects respectively. The coefficient β ⁎ p is always negative while β ⁎ T may be either positive or negative depending on the prevailing temperature effect. For the soft muon component (≤ 10 GeV) which composes the main part of the particles detected by our telescope in this experimental context, the negative temperature effect dominates and β ⁎ T is expected to be negative. The correlation analysis recently performed by Zazyan et al. 34 shows that the pressure and temperature effects are positively correlated. Consequently, for the present measurement conditions, both β ⁎ p and β ⁎ T are negative and time-correlated. When considering the atmospheric effect alone like in equation (2), the β p coefficient actually accounts for both the pressure and temperature effects. This explains why the experimental value found for β p (2) is larger than the value of β h (3).

Statistical feasibility and limits of opacity monitoring
We now address some statistical issues concerning the monitoring of opacity variations like those produced by water-level variations measured during the SHADOW experiment.
Let us assume that N = N 1 + N 2 particles are detected by the telescope during a time period T, and where N 1 and N 2 are the number of particles respectively counted during the first and second half of T. We want to determine under which conditions N 1 and N 2 may be considered different at the confidence level α. The particle flux difference ΔN = N 2 − N 1 obeys a Skellam distribution defined as the difference between two Poisson processes with means μ 1 and μ 2 with ε the flux variation percentage. When the inequality (6) becomes an equality we get T = T min , the minimum acquisition time necessary to resolve a flux difference given by the following set of parameters (φ 0 , ε, α). When ε is fixed, T min is the best time resolution achievable to observe temporal relative flux variations larger than ε. When T min is fixed, we derive the best relative flux variation, ε, detectable on a time-scale larger than T min . Note that if (N 1 , N 2 )  10 the Poisson laws can be approximated with Gaussians and equation (6) is simplified to, where α α =  erf( ). We numerically compute T min from equation (6) with a confidence level α = 0.05 and represent it on Fig. 6 for a range of measured muon flux and variation threshold ε = 1, 0.1, 0.01 and 0.001 (i.e. 100%, 10%, 1% and 0.1%). Observe that the approximation (10) is suitable for our range of applications, T min will be underestimated starting from ε  0.5 which implies N ≈ 20. Figure 6 shows that to detect a daily variation in the muons count of 2% (ε = 0.02), as is typically observed in the SHADOW experiment, an average flux φ 0 > 2 s −1 must be measured. This solution is represented by the black cross labelled "water tank" on Fig. 6. The lower-left domain delimited by the curved black arrow in Fig. 6 represents the 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 branch is limited by the experiment duration, and the vertical branch is placed at a level corresponding to the maximum flux that can be measured by the telescope. This latter quantity increases with the telescope acceptance, e.g. by increasing the angular aperture (i.e. by reducing the distance between the detection matrices), or by grouping several lines of sight, or by using instruments with a larger detection surface.
The feasibility domain for a typical volcano experiment is also represented on Fig. 6 and delimited by the blue curved arrow. Note that for this kind of experiments we have a longer acquisition time and a tiny measured flux as 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 (6) into a form more suitable for radiography applications by replacing the flux fluctuations by opacity fluctuations, where the muon flux φ, is explicitly written to depend on telescope acceptance  , opacity  0 , and zenith angle θ. As before, ε represents the variation of opacity relative to the average opacity  0 . We warn the reader that a given ε-variation of  corresponds to a much larger ε-variation of φ. Putting equations (12) and (13) in equation (9) we obtain the following feasibility condition, min 0 where T min is the measurement period minimum duration necessary to resolve the sought opacity variation. Note that the feasibility formula from Lesparre et al. 7 is the first order development of equation (14). A subset of T min solutions of equation (14) is represented on Fig. 7 for the confidence level α = 0.05, for zenith angles θ = 0°, 30° and 60° and opacity variations ε = 100%, 10%, 1%. An acceptance,  = 10 cm 2 sr, typical of our telescopes has been used in the computation. Roughly, a one order of magnitude ε variation induces a 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 is achievable, is minimum to resolve a given opacity variation. The optimal opacity range depends on the zenith angle and goes roughly from 40-100 m.w.e for θ = 0°, and from 20-40 m.w.e for θ = 60°. For low-opacity conditions, measurements at high zenith angles are wise to optimize the time resolution. This is particularly conspicuous for the SHADOW experiment where the average opacity  0 ≈ 5 m.w.e and ε ≈ 10%. For these parameters, Fig. 7 gives T min > 1 day is necessary at θ = 0° to resolve the fluctuations while T min > 0.2 day is sufficient at θ = 60°. The time resolution strong dependence with respect to the zenith angle disappears at larger opacities  0 > 500 m.w.e like those encountered in volcano muon radiography.

Discussion
Muon radiography is a powerful method to monitor opacity/density variations inside geological bodies. Noticeable advantages of the method are the possibility to remotely radiography unapproachable dangerous volcanoes and to image the density distribution of large volumes from a single view-point 8,10,36 . Muon radiography is entering an era of precision measurements not only for structural imaging but also for dynamical monitoring purposes. Some monitoring experiments have been performed on active volcanoes that demonstrate the usefulness of such measurements to constrain the evolution of eruption crisis 12 . However, as shown above, monitoring opacity variations is subject to external sources of bias, and statistical and experimental constraints that limit the achievable resolution. Understanding these limits is of primary importance to improve the method and to assess the muon radiography monitoring feasibility and validity.
Experimental constraints are partly dictated by statistical considerations, and mainly come from the telescope acceptance that limits the maximum flux which fixes the resolution domain right boundary in Fig. 6. This boundary may be moved rightward by increasing the acceptance  of the instrument. Recalling that  is expressed in [cm 2 sr], the acceptance may be augmented by several means: 1) increasing the solid angle encompassed by the instrument by reducing the distance between the detection matrices; 2) increasing the detection surface by coupling several telescopes (actually our telescopes may be merged into a single one); 3) grouping lines of sight to increase both the detection surface and the solid angle at the price of reducing the angular resolution of the radiographies. In the present study, the latter solution was retained and all lines of sight were merged to obtain an effective acceptance of 630 cm 2 sr.
Statistical constraints bound the resolution domain of a given experiment (Fig. 6), and the main concern when doing measurements is to ensure that the monitored phenomena fall inside the boundaries. As will be discussed in the 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 Fig. (7) where the feasibility solutions for T min strongly differ in the low-and high-opacity domains. It is remarkable that high-opacities variations are equally resolved whatever the zenith angle while, instead, the resolution for low-opacities strongly depends on this angle. Another conspicuous feature present in Fig. (7) is the existence of an optimal medium-opacity range  ≈ 50 ± 30 m.w.e. where telescopes offer their best performance. These two effects are due to the cosmic muon energy spectrum nature 15,17,37 , and changing the telescope acceptance has no effect on the optimal opacity values but only changes T min by translating the solution curves of Fig. (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 shown in Fig. 8 4,10,38 . The picture was taken during an open-sky calibration phase where the muon count serves to determine the efficiency of the scintillator bars forming the detection matrices. Each matrix is formed by an assemblage of two sets of 16 bars arranged perpendicularly to obtain a 16 × 16 square 5 × 5 cm 2 pixels array. The telescope upper and lower matrices allow 31 × 31 pixels combinations, i.e. 961 distinct lines of sight. The distance between the matrices may be changed to adapt the solid angle spanned by the trajectories. In the present study, the distance was tuned to encompass the entire water tank (Fig. 1).
Once geometrically configured, the telescope is totally characterised by its acceptance function i  [cm 2 sr] which relates the muon count, N i , to the muon flux, ∂ φ [s −1 cm −2 sr −1 ] received by the telescope in its i th line of sight, where T is the acquisition duration, i  [cm 2 ] is the line of sight detection surface function, i  is the integrated acceptance, and ∂ φ i is the muon flux in the line of sight central direction. It must be understood that ∂ φ [s −1 cm −2 sr −1 ] is the differential muon flux that reaches the instrument after crossing the target. Consequently, ∂ φ depends both on the open sky differential flux ∂ φ( = 0, ϕ, θ) and on the muon absorption law inside matter. These are determined through experiments 26,27,37,[39][40][41][42] , theoretical works 17 or thanks to Monte-Carlo simulations 43,44 depending on the precision expected and the available information. The three curves (resp. blue, red, green) correspond to three different observation zenith angles (resp. 0°, 30°, 60°) and are computed for  = 10 cm 2 sr using the modified Gaisser model from Tang et al. 17 .
Scientific RepoRts | 6:23054 | DOI: 10.1038/srep23054 Figure (9) shows the telescope acceptances i  for i = 1, … , 961 used in the SHADOW experiment. This acceptance function is determined experimentally to account for the detection matrices deffects, mainly imperfect optical couplings at the scintillator bars outputs and on the multichannel photomultiplier front. The latter   Fig. (1). The acceptance maximum value, max  = 2.80 cm 2 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 the upper detection matrix, and the other one in the lower matrix). The acceptance integrated over the instrument entire detection surface equals int  = 630 cm 2 sr for a solid angle aperture Ω int = 0.161 sr. causes the distortions visible in the Fig. (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 a subset ε , It results in an acceptance increase and thus a better time resolution. The counterpart is an angular resolution degradation induced by the merging of the small solid angles spanned by the trajectories. In the present study, the entire solid angle spanned by the telescope trajectories were grouped to obtain a total acceptance total  = 630 cm 2 sr. Such a large acceptance dramatically improves the time resolution which falls to the order of tens of minutes in the case of the SHADOW experiment.