Interpreting muon radiographic data in a fault zone : possible application to geothermal reservoir detection and monitoring

Rainfall-triggered fluid flow in a mechanical fracture zone associated with a seismic fault has been estimated (Tanaka et al., 2011) using muon radiography by measuring the water position over time in response to rainfall events. In this report, the data taken by Tanaka et al. (2011) are reanalyzed to estimate the porosity distribution as a function of a distance from the fault gouge. The result shows a similar pattern of the porosity distribution as measured by borehole sampling at Nojima fault. There is a low porosity shear zone axis surrounded by porous damaged areas with density increasing with the distance from the fault gouge. The dynamic muon radiography (Tanaka et al., 2011) provides a new method to delineate both the recharge and discharge zones along the fault segment, an entire hydrothermal circulation system. This might dramatically raise the success rate for drilling of geothermal exploration wells, and it might open a new horizon in the geothermal exploration and monitoring.


Introduction
Fluid flow through geological-scale porous media is of great importance in many areas ranging from civil engineering to geothermal technologies.Such large-scale porous media can be generated in response to lithostatic, tectonic, and thermal stresses and high fluid pressures.The high crack density in the damage zone results in significant pore volume increase, hence improving the water transport properties.Water influx is therefore often dominated by such high-conductivity fracture zones, and thus economically significant geothermal and water supply reservoirs often form in fractured rocks.For example, the fault-hosted geothermal water is captured in fractured rocks within seismic faults.Tunneling through a fault zone can produce a sudden water inflow if it intersects a high-pressure fracture zone.Fracture zones also affect stability, deformation, and fluid flow into underground structures.The penetration of water into fractures can be also hazardous since water-induced pore pressures may increase in response to rainfalls if the fractures are blocked, and as a result failure may follow.Understanding of the geometrical, porous and permeability structures of fracture zones would improve the design and construction of tunnels, dams, and geothermal plants.
A typical seismic fault structure includes a narrow finegrained core surrounded by a damage zone of highly fractured rock (Caine and Forster, 1999).A zone that includes the fine grain size in the fault core is called a fault gouge, and expected to have a low porosity.The surrounding damage zone contains a high density of microcracks, and thus is likely to have higher porosity.Conventionally, such a contrast between damage zones and the surrounding rock has been indirectly measured by magnetotelluric imagery (Lockner and Byerlee, 1986), by seismic tomography (Li and Leary, 1990), and more directly by drilling boreholes (Lockner et al., 2009).Although bore hole sampling method provides density contrast along the shaft with a great accuracy, it only provides local information around the drilling hole.
Muon radiography is becoming a commonly used imaging technique for geological or plant evaluation and monitoring (e.g., Tanaka et al., 2009).Radiographs are reproducible and can identify changes in density inside the target volume with fairly high precision and accuracy.However, there are several disadvantages to muon radiography.Firstly, technical limitations include the need for a topographically prominent feature of interest, and there will only be results for the volume located above the detector.This is why a seismic fault on a hill slope made a good study target in the work by Tanaka et al. (2011) (Fig. 1).However, this problem could be solved with borehole/tunnel observations from underground muon detectors (Tanaka and Sannomiya, 2013).Secondly, this technique only resolves the average density distribution along individual muon paths.Therefore, the user must end up making assumptions or interpretations about more localized structure along those muon paths, or must use more than one detector to resolve the three-dimensional density structure.
The muon attenuation is directly proportional to the density of the material with minimal uncertainty due to chemical composition.Muon radiography therefore provides varying densities inside the geological formations.Among various physical quantities measured in the earth's crust, density plays a special role because it is readily interpreted in terms of porosity of fracture-filling rocks; hence the permeability of the fault zones can be estimated from muon radiography.Tanaka et al. (2011) measured rainfall-triggered fluid flow in a fracture zone associated with a seismic fault located in Itoigawa-Shizuoka tectonic line with muon radiography.It has been shown that the water front was clearly imaged in response to rainfall events.There, correlations were observed between rainfall events and muon intensities for different depths in the fault zone.The decrease in the muon flux has been interpreted as a result of active fluids trapped inside cracks of the damage zone.After the rainfall event, the trapped fluid is gradually released from the cracks, and as a result the increase in the muon flux was measured.In this work, the data taken by Tanaka et al. (2011) were reanalyzed to estimate the density distribution as a function of a distance from the fault gouge.The method directly measures an average density along the muon path, covering a larger spatial extent.In this report, possible geothermal exploration and monitoring utilizing muon radiography is also discussed.

Reanalysis
Itoigawa-Shizuoka tectonic line (ISTL) is a major fault that runs from the city of Itoigawa, Niigata Prefecture, through Lake Suwa to the city of Shizuoka in Shizuoka Prefecture.The ISTL is a major geological structure dividing the Japanese mainland into northeastern and southwestern parts.A part of this fault line was discovered by trenching a fault scarp at Fossa Magna Park, the first UNESCO Geopark in Japan, located at the north end of the ISTL.The discovered fault cuts across 16-million-year-old basaltic rock.The fault outcrop on a steep hill slope provides an ideal location to study the properties of the fracture (Fig. 1).The data discussed here were obtained using a detection system with an active area of 3969 cm 2 that was installed 6 m south of the fault outcrop, pointing towards the geologically estimated direction of the fault line (Fig. 1) (Tanaka et al., 2011).
In order to identify the path of cosmic-ray muons, two segmented detectors were employed as position-sensitive detectors (PSDs) in this experiment.Each segmented detector consisted of two plane arrays of scintillator strips, one each in the x and y direction.Each scintillator strip used a plastic scintillator 70 cm long by 7 cm wide and 2 cm thick, and a 2-inch photomultiplier tube.The two arrays each had nine counters.The path of a muon can be determined by the combination of two signals from an x and a y plane detector, which defined a 7 cm square within which the muon passed.The path of a muon can be determined by the locations of its signals at the two segmented detectors, which were separated by a distance of 70 cm.The triggering requirement for the soft-component rejection analysis is a coincidence within 20 ns of at least two out of nine counters within the same plane (multiplicity cut analysis method) (Tanaka et al., 2001).Therefore within this framework, the selection criteria for a candidate muon path are as follows: there is one signal from each scintillator array plane; the two signals of the counter are detected simultaneously; and the four signals are in a coincidence window.An angle of the incident muon is calculated using the difference between detection positions of two PSD planes.
The angular resolution of the present segmented detection system was evaluated by utilizing a steel block (Tanaka, 2013).The root mean square (RMS) angular resolution of ±40 mrad defines the ±24 cm accuracy of the position resolution at the outcrop.However, at a distance of 100 m from the detection system, it becomes ±4.0 m.The present detection system cannot resolve the density of the fault gouge whose width is 60 cm.The image acquisition was made at a 27-day exposure time.The distance travelled in 2 × 10 −8 s is equivalent to that between the detection system and the outcrop, and is much shorter than muon's decay constant of 2.2× 10 −6 s.We therefore neglected the effect of muon decays.The time resolution of the data acquisition system is 10 ns, and the width of the coincidence window between two planes was set to be 40 ns.

Data analysis
Since we expect that the geomagnetic effect on cosmic-ray muon flux (east-west effect) is negligible within the detector's effective azimuth angle of −400 mrad < ϕ < +400 mrad, the azimuth distribution of muons from the sky (elevation of 600 mrad) was used to correct the data.Figure 2 shows the histogram of the muon events collected in the present experiment.The counts in negative elevation angles mainly consist of the "backward-directed" muons that can be used to confirm the integral muon intensity, its zenith-angular dependence, and the geometrical acceptance of the detection system.The azimuthally symmetric behavior of the counts confirms that the muons are azimuthally isotropic.The combination of the scintillator strips to determine the arrival direction is different between forward-directed and backward-directed events.It is more useful for us to normalize the raw azimuth distribution with the forward open sky flux.However, in our case, the target was within the entire detector viewing angles, and therefore the backward-directed events had to be used for normalizing data.By utilizing high yield scintillators and photomultiplier tubes (PMTs), and by setting the lowest discrimination level as possible, the efficiency loss of the scintillator strips is suppressed.As shown in Fig. 3, the azimuth distributions at different elevations of sky are almost isotropic for negative angles (backward-directed) as well as that from the thick part of the target (+100 mrad), where the events are dominated by the azimuthally isotropic EM shower-originated background.The raw radiograph was corrected by applying the different horizontal distribution at a different elevation in order to cancel the geometrical acceptance of the system, and compared with the result of Monte Carlo (MC) simulations.
The azimuth distributions at different elevations (θ = −700, −500, −400, −300, +700, +600, +500, and +100 mrad) were normalized to the reference azimuth distribution (θ = −600 mrad).Figure 3   caused by accidental events simultaneously hitting the two segmented planes of the detector, the extended air shower MC simulation compared the muon intensity arriving at θ zenith > 50 • (horizontal muons) in zenith and the intensity of multiple (≥ 2) shower particles/m 2 (> 10 MeV) arriving at θ zenith < 50 • (vertical EM particles).By considering the vertical and horizontal projected area of the detector, roughly 0.6 % of the horizontal muons can be counted as background noise.As shown in Fig. 2, the ratio of the number of events from θ zenith = 100 mrad and −700 mrad is ∼ 0.0057 (2000/350 000), which is consistent with the estimation.Therefore, it is very difficult for us to analyze the data for θ < +300 mrad.As shown in Fig. 3, a statistically significant excess can be seen at ϕ = +50 mrad.From the topographic information, the muon events will simply decrease as the azimuth angle increases if the density of the target is uniform.This can be a systematic error associated with the geometrical efficiency of the segmented detection system.However, such an excess cannot be seen for θ ≤ −300 mrad and θ = +100 mrad.
Figure 4 shows the derived density length (muon path length times average density along the path line) as a function of angle from horizontal.The estimated density length was compared with the actual path length as measured with the topographic map, assuming a uniform average density to find the lower density region.The lines in Fig. 4 show the thickness of the target, as obtained from the position and the resolution of the apparatus, and the local topography around the fault zone.The reading error of the map is ±1.5 m, and the horizontal uncertainty of the topography model is ±1.0 m.This error corresponds to 2.5 % for 100 m rock at maximum.The thickness of the line in Fig. 3   assumed, and density variations relative to the reference were determined.Overall, there is an agreement between the muon transmission rate and the local topography.By comparing the estimated density length with the actual path length, the average density along the path line can be calculated as shown in Fig. 5, where it is plotted in g cm −3 against the horizontal distance from the gouge.Statistically well-determined three data points at ϕ = ±50 and ±150 mrad in Fig. 4 were used to draw Fig. 5.As shown in Fig. 3, since the number of fake track N BG is of the same order of magnitude as true muon tracks for θ < 400 mrad, the statistical error on the values for these angular regions is seriously masked by the noise.
In Fig. 5, the azimuth angle was converted to the horizontal distance at a distance of 100 m from the detector, representing the middle point of the target whose thickness is ∼ 200 m above 300 mrad.It was assumed that the width of the clay gouge (60 cm, which is not resolvable with the present detector) is constant throughout the strata.The actual density of 2.2 g cm −3 was taken as the density of the clay gouge.In this figure, the density values at ±5 m are identical for every elevation angles, because both values were taken from the data at ϕ = 50 mrad.It should be noticed that the density values at ±5 m include the density of the clay gouge, but its contribution is negligible.

Results and discussion
Given that the data point at ϕ = 50 mrad in Fig. 4 must, by the detector orientation, pass through the fault gouge (or the narrow fine-grained core), we see that systematic decrease in the density length (or the average density along the muon path in Fig. 5) from the assumed initial density length (or the average density) around the fault gouge.The center of the fracture zone, due to its high clay content and fine gouge grain size, has high density in general.However, we do not see such a high density in Fig. 5.This is because the resolution of the detector is much lower to capture the gouge width, and thus the density was averaged over a larger spatial extent.However, we can assume that the contribution from the clay gouge to the average density is small.
We also can see that, farther away from the fault gouge, density increases.This is because grain damage decreases as a function of the distance from the fault core.The damage zone half width, as indicated by the density profiles, is ∼ 10 m.For reference, a fault zone half width of 10-25 m for the Nojima fault was suggested by Lockner et al. (2009).The result shows that the technique can be applied to monitor the process of sealing and restrengthening of the seismic fault.It has to be stressed that dynamic muon radiography may provide a good opportunity to monitor how rapidly the enhanced fault zone porous structure will be reduced by sealing and crack healing processes.
Most geothermal reservoirs consist of fractures including faults, where permeability is a key factor for estimating the potentials of the hydrothermal convection.The permeability of fractures is theoretically and empirically a function of the porosity as well as porosity-induced density of fracturefilling rocks.The density recovery pattern as a function of the distance from the fault axis demonstrates that muon radiography can be one of the most powerful exploration techniques that directly delineate the permeability structure of the fault zones in detail (Fig. 5).
Drilling cost of each geothermal well is estimated in the range of 2-20 million US dollars (Tester et al., 2006) and is the main interest in geothermal exploitations.Nevertheless, available major exploration methods to determine a target fracture by imaging subsurface geothermal reservoirs until the present are confined only to the resistivity and seismic surveys, where they use electromagnetic and elastic waves, respectively, both of which do not go straight on.The positioning resolutions of imaging by these methods are not so satisfactory such that the drilling of exploration wells in greenfield areas is reported to have a success rate of typically about 50 to 60 % (Goldstein et al., 2011).About half of the drilling cost is wasted.
Even in a single fault segment, one part acts as a recharge zone and the other part acts as a discharge zone depending on the position of geothermal heat sources and piezometric gradient induced by topography along the fault segment.The dynamic muon radiography by Tanaka et al. (2011) and this paper provide a new method to delineate both the recharge and discharge zones along the fault segment, an entire hydrothermal circulation system.Dynamic muon radiography also provides insights on the flow directions of hydrothermal systems where fluid flow prevails along the faults through both sides of permeable zones of the central shear zone, but fluid flow is impeded across the faults by their fine-grained central shear zones.As a similar concept, it has been shown that dynamic muon radiography has the potential to be a practical technique for monitoring the movement of carbon dioxide during carbon storage in deep geological formations (Kudryavtsev et al., 2012).

Conclusions
In conclusion, when the geology is well characterized, the present method is likely to give useful porosity measurements, as otherwise differences in density may relate to different rock types.Under this condition, the technique we reported in this paper can be used in ordinary soil where geothermal reservoirs are expected.The density recovery pattern as a function of the distance from the fault axis is consistent with the borehole sampling result.However, the present target is limited to a topographically prominent feature of interest, because we can image only the volume located above the detector.
Borehole muon sensors are one of the candidates that can be deployed down drill holes to survey an underground reservoir, as shown in Fig. 6.However, the muon detector should be compact, watertight and water-resistant to large temperatures and pressures and should have a size that can be inserted in a standard drill hole.We anticipate that high-resolution borehole muon sensors have a potential to open a new horizon in the geothermal exploration and monitoring. 17

Fig. 2 .
Fig. 2.Muon events collected in the present experiment as a function of azimuth and elevation angle.The horizontal and vertical bin width is 100 × 100 mrad, larger than the experimental angular resolution, which was estimated to be ±40 mrad (root mean square).

Fig. 4 .
Fig. 4. Elevation dependence of the density length as a function of azimuth angle.The solid lines show the density lengths as measured from the topographic map by assuming a uniform density of 2.0 g cm −3 .The measurements and the related error bars can be referred to the scale and unit on the vertical axis.The bottom axes indicate the azimuth angles (ϕ in units of mrad) with reference to the line perpendicular to the detector plane and path lengths respectively.

Fig. 5 .
Fig. 5. Density distribution as a function of the distance from the fault gouge.The vertical error bar shows the statistical error, and horizontal error bar indicates the positioning resolution of the detector.The data points are fitted to guide the eye for visible convenience.