Simulation of muon radiography for monitoring CO 2 stored in a geological reservoir
J. Klinger, S.J. Clark, M. Coleman, J.G. Gluyas, V.A. Kudryavtsev, D.L. Lincoln, S. Pal, S.M. Paling, N.J.C. Spooner, S. Telfer, L.F. Thompson, D. Woodward
SSimulation of muon radiography for monitoring CO stored in a geological reservoir J. Klinger a, ∗ , S.J. Clark b , M. Coleman c , J.G. Gluyas d , V.A. Kudryavtsev a , D.L. Lincoln e , S. Pal a , S.M. Paling f , N.J.C. Spooner a ,S. Telfer a , L.F. Thompson a , D. Woodward a a Department of Physics and Astronomy, University of She ffi eld, She ffi eld, S3 7RH, UK b Department of Earth Sciences, Durham University, Durham, DH1 3LE, UK c NASA Jet Propulsion Laboratory, California Institute of Technology, Pasadena, California, USA d Durham Energy Institute, Durham University, Durham, DH1 3LE, UK e Department of Civil and Structural Engineering, University of She ffi eld, She ffi eld, S1 3JD, UK f STFC Boulby Underground Science Facility, Boulby Mine, Cleveland, TS13 4UZ, UK
Abstract
Current methods of monitoring subsurface CO , such as repeat seismic surveys, are episodic and require highly skilled personnelto acquire the data. Simulations based on simplified models have previously shown that muon radiography could be automated tocontinuously monitor CO injection and migration, in addition to reducing the overall cost of monitoring. In this paper, we presenta simulation of the monitoring of CO plume evolution in a geological reservoir using muon radiography. The stratigraphy in thevicinity of a nominal test facility is modelled using geological data, and a numerical fluid flow model is used to describe the timeevolution of the CO plume. A planar detection region with a surface area of 1000 m is considered, at a vertical depth of 776 mbelow the seabed. We find that one year of constant CO injection leads to changes in the column density of (cid:46) plume is already resolvable with an exposure time of less than 50 days. Keywords:
Muon radiography, CCS, Carbon capture monitoring, Carbon capture, Cosmic-ray muons
1. Introduction
The regulation of atmospheric greenhouse gas concentra-tions is required to stabilise the e ff ects of anthropogenic cli-mate change. Most scenarios of economic development pre-dict a steep rise in greenhouse gas concentrations due to an in-crease in demand for energy, and hence fossil fuel usage [1].The option of CO Capture and Storage (CCS) is regarded asone technological solution for atmospheric CO mitigation [2].CCS technologies have the capacity to reduce CO emissionsfrom power stations by up to 90% [3], by the injection of lique-fied CO into on-shore or o ff -shore geological repositories forpermanent storage.The principle of transporting and storing CO in depleted oilfields or saline formations is already well understood, althoughthere are technological and commercial issues which need to beaddressed [4]. In particular, there is a demand for a continuousand economically viable method for monitoring the injectionand migration of subsurface CO . The migration of subsur-face fluids can be highly unpredictable, even in developed oilfields [5]. Current monitoring methods are episodic and requirehighly skilled personnel to acquire the data. Such technolo-gies include repeat seismic surveys, measuring the subsurfaceresponse to electromagnetic waves, measuring specific gravityand monitoring pressure [6]. Any new, and preferably low-cost,technology which is able to provide automated and continuous ∗ Corresponding author
Email address: [email protected] (J. Klinger) monitoring will greatly enhance the long-term viability of CCStechnologies. One such option that has been proposed is muonradiography [7].High-energy muons, which are produced due to the interac-tions of cosmic-rays with the Earth’s upper atmosphere, aredeeply penetrating charged particles. Muons are point-likecharged particles which are very similar to electrons, but havea rest mass approximately 200 times that of the electron. Allcharged particles interact with other charged particles via theelectromagnetic interaction. In this paper we consider muonsinteractions within matter, through which they are travelling.Although the materials in the Earth are electrically neutral, theatoms within the materials are composed of positively chargednuclei and negatively charged electrons. When charged parti-cles, such as cosmic-ray muons, interact with nuclei and elec-trons they scatter and lose energy by ionising and excitingatoms, as well as occasionally producing secondary particlesthat will be absorbed in matter quickly. As with any stochas-tic energy transfer in a bulk material, the energy is eventuallydispersed as heat.There are other cosmic-ray particles which one can considerat the surface of the Earth. Other point-like cosmic-ray par-ticles such as electrons and photons lose energy at a muchfaster rate compared to muons underground, due having a sig-nificantly lower rest mass (or zero rest mass in the case of the Although photons are electrically neutral particles, they are the ‘force-carriers’ of the electromagnetic interaction, and therefore by definition interactwith electrically charged particles.
Preprint submitted to Elsevier October 22, 2015 a r X i v : . [ phy s i c s . g e o - ph ] O c t hoton) compared to muons. Cosmic-ray particles which havea more complex internal structure (and are therefore not point-like) such as protons, neutrons, mesons and nuclei can addition-ally interact in the Earth via hadronic interactions. This addi-tional interaction significantly reduces their range undergroundcompared to muons. Therefore muons are the only charged par-ticles produced in the atmosphere which are observed in experi-ments based in underground observatories. There is a relativelylarge flux of neutrinos underground, although their extremelylow interaction probability means that they can be neglected.Muons lose energy as they traverse dense materials and, oncethey have a su ffi ciently low energy, will stop and be absorbedinto the material or decay to produce an electron which willalso be absorbed. Muon energy loss is fundamentally relatedto the density of the medium through which they travel. It isalso related to the mean atomic number and weight of the bulkmaterial through which they pass. It is therefore possible to mapthe density profile of a large object by observing local variationsin the flux of muons emerging from the target object.Muon radiography is already used for a wide variety of mon-itoring purposes, from relatively small-scale applications suchas the monitoring of cargo containers for high- Z materials [8–11], to the mapping of density variations in large-scale struc-tures such as voids in the ancient pyramids [12] and magmachambers in volcanoes [13, 14]. Recent research has shownthat muon radiography may be able to provide a low-cost, con-tinuous subsurface CO monitoring system for the injection andsubsequent movement of CO , as an alternative or complemen-tary system to seismic technology [7].In this paper, we will show the results of the first detailed sim-ulation of muon radiography for this application. Our model in-corporates geological data to describe a subsurface CO storagesite, with a realistic model of subsurface CO plume evolution.Such models will be an essential part of any future use of thistechnology. The experience of fluid injection processes in theoil industry suggests that natural systems do not necessarily be-have as anticipated. An example of this is the use of waterfloodto enhance oil production, which can lead to the rapid break-through of water without oil [5]. Thus, in the context of CCS, itwill be essential to have a detailed model of expected behaviourso that deviations from those expectations can be identified atthe earliest stage.It is clear that muon radiography will not be able to imagesubsurface CO with a higher resolution than technologies suchas seismic surveying. We propose, however, that muon radio-graphy will be su ffi ciently accurate such that it may reduce thetotal cost associated with CCS, if used in combination with ex-isting technology. Therefore the utility of muon radiography asa tool to monitor subsurface CO will ultimately depend uponthe practicality of deployment of the technology at CO injec-tion sites and the cost of deployment relative to other moni-toring technology which could be used. We will discuss theseaspects in more detail in Section 8. The monitoring of high- Z materials typically utilises the large scatteringangle of muons incident on high mass nuclei, rather than changes in tshe muonflux due to muon attenuation, although Ref. [8] also uses muon disappearance. We will show that this technology is conceptually viable, andas such this paper represents a crucial milestone in the realisa-tion of muon radiography for the purpose of CCS.
2. Conceptual overview
Projected implementations of CO capture technologiesutilise pre-combustion, post-combustion [15] or oxyfuel pro-cesses [16, 17]. Independent of the setup, liquefied or supercrit-ical CO can be subsequently piped to an on-shore or o ff -shoregeological storage site, of which there are already numerous ex-amples [4]. In the context of this paper we consider o ff -shoresites, which may be more commercially attractive due to theprospect of enhanced oil recovery [18, 19].In the current examples of this technology, liquefied or su-percritical CO is injected underground via a vertical well to ageological reservoir. It is required that the reservoir is sealedby an impermeable stratum (‘cap-rock’), in order to contain thesubsurface CO volume, which will rise under buoyancy forces.In this work, the geological storage reservoir considered isa saline aquifer, although the simulation procedures would besimilar for a depleted oil or gas field with residual hydrocarbonsand brine. As a consequence of constant CO injection, the den-sity distribution of the reservoir will evolve. To model this, weuse a numerical multiphase fluid flow model, which we describein detail in Section 4, following the approach of Refs. [20, 21].The model describes the vertical and horizontal distribution ofCO at discrete time steps, which evolves spatially from theCO injection point. Crucially, the model computes the changein the local density due to CO and brine distributions withinthe reservoir. As a result of the evolving density profile, whichalters the survival probability of passing muons, one should ex-pect localised modifications to the muon flux at depths belowthe storage reservoir.We propose that muon detectors are placed in containers thatmeet standard dimensions used for oil-well boreholes. Full de-tails of the detector dimensions are given in Section 6. The de-tectors are assumed to be deployed in horizontal arrays belowthe reservoir, such that each array of detectors is fully containedwithin a surface area of approximately 1000 m . We assume forthis study that each borehole container will contain a numberof bars of plastic scintillator. Plastic scintillators, which emitlight in response to the energy depositions of charged particlessuch as muons, are a low-cost and durable technology which iswidely used in particle physics research.Given the relatively low flux of muons at the depths whichare considered in this study, we employ a simple imaging strat-egy. We propose to resolve changes in the density of the to-tal geological section via the change of the underground muoncount due to CO injection. Any background to the muon sig-nals underground can easily be suppressed by requiring a mini-mum number of scintillator bar hits. As the e ff ect of CO injec-tion underground is to reduce the overall density of the aquifer,one can measure an enhancement of the number of muons be-low the reservoir, with respect to the number of muons that onewould expect in the absence of subsurface CO .2 urface area at sea level = 40.8 km
776 m
1: Lias shale, ρ = -3
2: Mercia mudstone, ρ = -3
3: Bunter sandstone, ρ = -3
4: Potash, ρ = -3
5: Evaporites, ρ = -3 Figure 1: A diagram of geological stratigraphy model used in this simulation. The layer number, simplified composition andaverage bulk density is indicated for each layer. There is additionally a layer of sea water above layer 1, with a depth of 32 m.Layer ‘3’ corresponds to the saline aquifer, which is a suitable storage medium for CO . The approximate location of the proposeddetector site is indicated by the ‘X’ symbol. In this simulation, CO is piped underground via a vertical well, with the point ofinjection vertically above the detector site.
3. Geological modelling
An unfaulted geocellular model of the geology in the regionof Boulby Mine, situated in Cleveland, United Kingdom, hasbeen produced for this study. There are a number of reasonswhy we have performed the simulation at this location, despitethe site not being an option for o ff -shore CCS. Firstly, the geol-ogy of the site at Boulby is well understood. The stratigraphyfeatures a permian evaporite layer at a depth of 0.75–1.1 km,which is typical of geological repository sites. Furthermore,there is safe, supported and versatile access provided by Is-rael Chemicals Ltd UK, who operate a commercial facility inBoulby Mine, and the STFC Boulby Underground Science Fa-cility [22], which allows for comprehensive testing of prototypedetectors in the underground environment.The geocellular model describes five strata below sea level,with a detector situated 776 m below the seabed, which is theapproximate depth of a proposed subsurface muon detector test-site. The five strata are referred to sequentially as ‘layer 1’,with the shallowest depth, to ‘layer 5’, with the greatest depth.There is additionally a layer of seawater above layer 1, with adepth of 32 m vertically above the test-site. The data samplewhich is used to describe the geological stratigraphy of this siteis described in Section 3.1, and the parameters which are usedto describe the saline aquifer, as a hypothetical CO storage site,are described in Section 3.2.Figure 1 depicts the model of the stratigraphy used in thissimulation. The approximate location of the detector site isindicated in the figure. The implementation of the geologicalcomposition and strata boundary (‘horizon’) data in the simula-tion framework is described in Section 5.1. Data describing four seismic horizons have been provided byIsrael Chemicals Ltd UK. These horizons correspond to the topof layer 1 (the seabed), the top horizon of layer 4, the commonhorizon of layers 4 and 5, and the bottom horizon of layer 5.In order to represent the stratigraphy of layers 1–3, two addi-tional layers were added to the model, based on the estimatedthickness taken from well data. The horizon data are converted into a triangular-based meshin order to reduce CPU and memory overheads. The averagedistance between points on the mesh is approximately 100 m,which represents the level of accuracy of this model. The hori-zon data are combined with information describing the compo-sition of each layer. The average bulk density and simplifiedcomposition of each layer are indicated in Figure 1. This in-formation is used to populate each layer in the simulation de-scribed in Section 5.
A numerical model is employed to describe the injection ofCO and the subsequent density evolution of the geologicalreservoir. The typical properties of the storage reservoir, whichare required to parameterise the model, are described in thissection.The numerical model, which is described in detail in Sec-tion 4, is based on a system of equations which define a uni-form axisymmetric region of the geological reservoir whichfully contains the volume of subsurface CO . The fraction ofthe reservoir containing void spaces, initially occupied by the in situ brine, is given by the porosity φ . The injected CO isentirely confined above and below by impermeable horizontalboundaries, separated by a distance H . This also defines theheight of the CO injection well boundary, which is centred onthe axis with a radius r . A constant mass-rate of CO injection, M n , is imposed which outwardly displaces the in situ brine.Accurate equations of state [23, 24] are used with referencevalues of pressure and temperature for the reservoir system, inorder to determine the corresponding fluid properties as listedin Table 1. The rock intrinsic density and porosity values aredetermined from a log of the test storage sandstone layer at thesite of Boulby Mine, Cleveland, United Kingdom. The parame-ters which follow are typical experimental values characterisinga sandstone-brine-CO system following the study of Ref. [25],where further descriptions of their meaning are also given. Typ-ical values for well depth below an impermeable boundary andwell radius are given along with a feasible mass-rate of injec-tion, which is at a lower limit of valuation deemed practical(3–120 kg / s) for commercial CCS purposes [26].3able 1: The parameters used to model the CO injection andmultiphase interactions in the geological reservoir presented inthis analysis. All parameters are assumed to have typical values,except for those indicated by ‘ * ’, which have been acquiredfrom data provided by Israel Chemicals Ltd UK. Parameter Sym. Value UnitsBrine density ρ w / m Brine viscosity µ w µ Pa sBrine bulk modulus K w density ρ n
720 kg / m CO viscosity µ n µ Pa sCO bulk modulus K n * ) ρ s / m Porosity ( * ) φ k . × − m Brine residual saturation S rw end-point relative permeability k rn m k relative permeability exponent n k m v p v / reservoir height H
170 mWell radius r mass-rate of injection M n
20 kg / s
4. Numerical subsurface CO simulation Numerical models are widely used in order to solve large-scale fluid flow problems in porous media. In this section, wewill present the formulation and implementation of a numeri-cal model for the simulation of a geological storage reservoirover discrete time steps after the start of a period of constantCO injection. This model will then be used in Section 5.3 toassess to what extent bulk density changes in such a system canbe tracked over time using muon radiography.We present a simple numerical model based on the coupledtheoretical approach of Ref. [21]. The model is then embeddedinto our G eant [27], described in Section 5,in order to account for the coupled interaction of CO and brinewithin the reservoir, and hence the corresponding macroscopicchanges in composition and density.In Section 4.1, we formulate the numerical model which wewill use to generate subsurface CO plume formations, usingthe input parameters described in Section 3.2. The bulk densitydistribution of the reservoir after CO injection, which representthe solutions of this model, are discussed further in Section 4.2.One should note that we only intend for the model to be usedin this simulation of muon radiography for the application ofmapping density changes in a large overburden. It is not in-tended that this model could be used for predictive purposes G eant at this site. The site that we model is anyway unsuitable forCO storage, rather it represents a site for which a completedataset was available and, as discussed in Section 3, contains aserviced laboratory for testing prototype detectors. The storage of CO in a deep saline aquifer, generally ofmoderate to high temperatures and pressures, presents a two-phase fluid system within the host rock; namely a ‘non-wetting’and ‘wetting’ phase. The non-wetting phase, modelled as a su-percritical or liquid CO -rich fluid, displaces the wetting-phase,which is modelled as a liquid H O-rich fluid, within a poroussolid (rock) phase.We assume that bulk composition and density changes areprimarily due to the multiphase displacing and drainage be-haviour within the reservoir system. We therefore concentrateon modelling the saturation e ff ects, whilst neglecting the mul-tiphase miscibility, thermal, inertial and solid deformation ef-fects, which is typical for a first-order reservoir analysis.Considering the continuity of multiple fluid phases in aporous medium, the following general macroscopic mass bal-ance equation is given for each phase π , ∂ ( φ S π ρ π ) ∂ t + ∇ · ( φ S π ρ π v π ) = , (1)where φ is porosity (the fraction of void volume to total vol-ume), S π is fluid saturation (the fraction of fluid phase volumeto void volume, such that (cid:80) π S π = ρ π is intrinsic phase den-sity and v π is velocity. The relevant constitutive relationshipsfor bulk fluid compressibility and flow respectively are:1 ρ π ∂ρ π ∂ t = K π ∂ p π ∂ t and φ S π v π = k r π k µ π ( −∇ p π + ρ π g ) , (2)where K π is bulk modulus, p π is fluid pressure, g is the gravityvector, and k , k r π and µ π are the extended Darcy’s law terms,for multiphase flow of intrinsic rock permeability, relative per-meability and dynamic viscosity respectively. The constitutiverelationships are substituted into Equation (1) by assuming φ constant and dividing through by ρ π , whilst neglecting the gra-dient of ρ π , giving φ S π K π ∂ p π ∂ t + φ ∂ S π ∂ t + ∇ · (cid:34) k r π k µ π ( −∇ p π + ρ π g ) (cid:35) = . (3)The fluid saturation capacity relationship is introduced: ∂ S w ∂ t = ∂ S w ∂ p c ∂ p c ∂ t = ∂ S w ∂ p c (cid:32) ∂ p n ∂ t − ∂ p w ∂ t (cid:33) , (4)where S w ( p c ) is the saturation of the wetting phase, which iscontrolled by the capillary pressure p c . The capillary pressureis defined as the di ff erence in pressure between the non-wettingand wetting phases, p c = p n − p w . In this formulation π = n , w, which subscript the terms belonging to thenon-wetting and wetting phases respectively. Ω and on its boundary Γ , for whichthe usual initial and boundary conditions are incorporated [21]: (cid:34) C w QQ T C n (cid:35) dd t (cid:40) p w p n (cid:41) + (cid:34) H w
00 H n (cid:35) (cid:40) p w p n (cid:41) = (cid:40) f w f n (cid:41) (5)where C π = (cid:90) Ω N T (cid:32) φ S π K π − φ ∂ S w ∂ p c (cid:33) N d Ω , Q = (cid:90) Ω N T (cid:32) φ ∂ S w ∂ p c (cid:33) N d Ω , H π = (cid:90) Ω ( ∇ N ) T k r π k µ π ∇ N d Ω , f π = (cid:90) Ω ( ∇ N ) T k r π k µ π ρ π g d Ω − (cid:90) Γ N T q π ρ π d Γ , (6)are compressibility, coupling, permeability and supply matri-ces respectively. The N terms are a set of linear shape functionsinterpolating over the discretised domain. The mass flux, q π ,imposed normal to the boundary, is also introduced as a conse-quence of the boundary conditions.To model the injection scenario, the system of equations (5)are solved with the following specific boundary conditions. Theinner (inflow) vertical injection boundary is prescribed with amass-flux boundary condition, q n , which is related to the mass-rate of CO injection by M n = q n π rH . At the outer (out-flow) radial extent of the domain, a constant hydrostatic pres-sure far-field boundary condition is prescribed. Otherwise, no-flux boundary conditions are assumed.This system of equations is non-linear due to the dependenceof the coe ffi cient matrices on the primary variables; namely, thesaturation and relative permeability terms, S π ( p c ) and k r π ( S π ).These relationships are determined experimentally for a givensystem and are parametrised in this work by the van Genuchten S π - p c model [29] and by a power law k r π - S π relationship [25],respectively (see Table 1). The model parameters, φ , k , ρ π , µ π and K π are assumed constant and uniform, which is typical forbasic hydro-geological investigations.Once the primary pressure variables are determined for agiven point in time, the corresponding secondary saturationvariables S π ( p c ) are determined. From this, the followingmacroscopic bulk density can be given by summing the intrin-sic phase densities multiplied by their corresponding volumefractions, ρ b = (1 − φ ) ρ s + φ S n ρ n + φ S w ρ w , (7)which introduces the intrinsic density of the solidrock / minerals, ρ s . The spatial integrations of (6) are carried out using Gaus-sian quadrature. The temporal integration of (5) is carried outusing finite di ff erencing; the set of equations are solved mono-lithically using a fully implicit (unconditionally stable) embed-ded scheme producing solutions of adjacent order in accuracyto allow for error control and adaptive time-stepping. This isconsidered appropriate given that a coupled injection scenariois being modelled. Each time-step is solved iteratively due tothe non-linearities, which is done via an accelerated fixed-pointtype procedure until convergence is met. The linearised sys-tem is solved at each iteration by a standard direct multi-frontalsolver. Further details on the implementation of a system ofequations of this type are given in Refs. [21, 28, 30]. The solutions to the model presented in Section 4.1 describethe density distribution of the reservoir at time t after the startof continuous CO injection. The radial extent r E ( t ) of the ex-panding volume of CO has the relationship r E ∝ √ t [31].It is found that density changes of (cid:46)
4% are expected due tothe presence of injected CO , with respect to the bulk density ofthe reservoir pre-injection. For the purposes of muon radiogra-phy, which is sensitive to the density of the total geological sec-tion, this corresponds to a change in the total column density of (cid:46) eant
5. Simulation framework
The simulation framework that is used for muon transport, inwhich muons descend through the geological model from Sec-tion 3, is described in this section. All stages of the simulationare performed within the G eant eant ρ b , of thereservoir due to the presence of CO . In Section 5.3, we presentthe modelling of the muon flux distribution that we use at sealevel, and the subsequent transport to the detector site. The model of the stratigraphy of the test-site presented inSection 3 is interfaced with G eant
4. The implementation ofthe geological model in G eant eant
4, which allows formulti-threading to be used in order to reduce total computationtime. We construct each stratum of rock using two parallel (andvery approximately horizontal) meshes, bound by four verticalfaces. Each mesh represents one of the horizons described inSection 3. We opt to construct the less detailed horizons usingquadrangular-based meshes and the highly detailed horizons us-ing triangular-based meshes.5he muon transport simulation, which is described in detailin Section 5.3, is performed in two stages: firstly only trans-porting muons to the lower horizon of layer 2 (or, equivalently,the top horizon of layer 3); and in a subsequent stage of simu-lation, transporting muons through the remaining three layers.One reason for this, in addition to a general improvement inCPU performance, is that fewer constructed strata are used inG eant
4, which therefore reduces the memory overhead. Fur-thermore, as the first stage of the simulation only considersmuons above the subsurface CO formation, it is su ffi cient toperform this stage of the simulation only once and then per-form the second stage of the simulation once for each time stepof subsurface CO migration. The numerical model describing the variation in bulk densityof the reservoir, presented in Section 4, is implemented in thesimulation using the voxelisation capabilities of G eant
4. Inthis framework, a regular grid of voxels (each of size { × × } m ) populates the domain of interest, which is rendered andparameterised by the numerical model.We take steps to preserve the accuracy of the numericalmodel in the voxelisation procedure, particularly in regions ofhigh numerical gradients caused by the fluid interface, whichare due to the displacing behaviour of the CO . The first stageof the procedure is to overlay the finite element mesh, whichdescribes density profile and gradient with high precision, ontothe voxel grid. The voxels are then parameterised with the com-position and density outputs, which are interpolated from thefinite element mesh at all of the overlapping voxel centroids.The voxelisation process is carried out at selected points intime. As described in Section 4.2 and Ref. [31], the radialextent r E ( t ) of the CO distribution expands with time t as r E ∝ √ t . Accordingly, the voxelisation is performed at inter-vals (2 n ) days which are therefore assumed to represent thereservoir state for the period between (2 n − and (2 n + days, where n = →
10. For this simulation, using the pa-rameters presented in Table 1, the radial extent is increased byapproximately 36 m at each time-step until reaching a distanceof approximately 360 m from the well.A quarter-symmetric rendering of the system bulk densityas predicted at 64 days, alongside its equivalent voxelisationof density bins for implementation with G eant
4, is shown inFigure 2.
For the simulation of muon energy loss in solid material, weuse the G eant E at sea level, as a function of theangle θ of the muon trajectory from the zenith. The spectrumis based on the Gaisser parameterisation [32], which accountsfor correlations between the θ and E due to di ff erent muon-production mechanisms in the atmosphere. The parameterisa-tion accounts for muon energy losses in the atmosphere, which (a) Before voxelisation (b)
After voxelisation
Figure 2: Quarter-symmetric rendering of the sandstone-brine-CO system bulk density as predicted after 64 days (a), along-side its equivalent voxelisation of density bins for implementa-tion with G eant E >
400 GeV so we do not include recent measurementsof the muon spectrum such as Ref. [33] in which the spectrumis only measured up to 100 GeV. Muons which have been gen-erated according the Gaisser parameterisation and subsequentlypropagated to these depths have been previously shown to agreewith data within <
10% [34]. We apply additional correctionsto the parameterisation to account for the muon lifetime and theEarth’s curvature.In order to reduce computation time, a number of optimisa-tions are implemented. There are no charged particles, otherthan muons, which can survive through large depths of mat-ter. Therefore any secondary particles that are produced due tothe interactions of muons with matter are immediately removedfrom the simulation.We only consider muons with E >
400 GeV and θ < ◦ ,based on a previous study of the muon survival depth [34, 35].Given the detector depth of ∼
800 m below sea level, and thatthe detector array is confined to within an area of 1000 m ,we require that muons originate from an approximately squareregion of { . × . } km at sea-level. This surface area is suf-ficiently large to account for any displacement of the muon’smeasured position r real due to scattering during the muon trans-port, with respect to the position that is linearly projected fromthe muon’s trajectory at sea level, r proj . The distribution of | r real − r proj | of muons in this analysis are shown in Figure 3. Wefind that | r real − r proj | is a steeply falling distribution, with the av-erage displacement of muons being < . d max , using a look-up table. The look- We define ‘primary’ particles as those which are present at sea level, and‘secondary’ particles as those which are produced at any depth below this. [m] proj r real |r Fraction of surviving muons
10 1
Figure 3: The distribution of | r real − r proj | between the horizontalmuon position at the detector site, r real , and the position that islinearly projected from the muon’s trajectory at sea level, r proj .up table, which maps E to d max , is generated with the MU-SIC muon simulation code [34, 35]. We then redefine d max ,as a function of E , as the distance beyond which a muon hasa survival probability of less than 10 − . The MUSIC codeonly considers materials with a uniform density, so we con-servatively choose to evaluate d max for a material with density ρ = .
17 g / cm , which corresponds to the lowest density ofthe five strata (layer 4). We therefore remove muons at anypoint in the simulation which have a remaining transport dis-tance greater than d max . One should note that the uniform den-sity assumed in the MUSIC simulation code is only used toprovide a conservative optimisation to the simulation and doesnot a ff ect the detail of our model.The energy spectra of muons, at sea level and undergroundrespectively, for muons which survive to the detector are shownin Figure 4. The average energy of muons which reach the de-tector site is found to be approximately 1500 GeV at sea level,and 220 GeV underground. The flux of muons at the detectorsite is found to be approximately 2 . × − cm − s − .
6. Muon detectors
In this section, we discuss the general configuration of muondetectors that we propose for mapping changes in the density ofthe total geological section. In the analysis of the muon simu-lation data, we present the maximum possible signal; thereforeneglecting gaps between the detectors and assuming 100% ac-ceptance for detecting muons whilst rejecting background fromradioactivity and detector noise. The discussion presented inthis section is therefore intended to provide a degree of quan-tification of the total e ffi ciency of these e ff ects, whilst still al-lowing us to present our results in a manner that is independentof the detector configuration. In order to ensure the long-term stability of the monitoringsystem, we assume that plastic scintillators, coupled to appro-priate photo-sensors, will be employed as muon detectors for
Muon energy [GeV] Surviving muons / GeV
10 110 Sea levelUnderground
Figure 4: The energy spectra of muons, at sea level and under-ground respectively, for muons which survive to the detector. Intotal, approximately 140 days of simulation data is presented,which corresponds to approximately 2.8 × muons in eachspectrum. The same muons are used in both distributions, there-fore only muons at sea level with energy greater than 400 GeVare observed underground.this application. We assume that individual detectors will beconfined to containers capable of being deployed in standardoil-well boreholes. The borehole containers would typicallyhave an inner diameter of 15 cm and su ffi cient length (lessthan a few metres) for commercially available plastic scintilla-tor bars. Plastic scintillator bars can be acquired commerciallywith a cross-sectional area of { × } cm and a length of typi-cally 100 cm. It is envisaged that multiple scintillator bars canbe packed into borehole containers. The borehole containerscan be deployed in a horizontal array at some depth below thereservoir.Plastic scintillators emit light in response to charged parti-cles, such as muons. The light signal is converted into an elec-trical signal with photo-sensors, such as photomultipliers, andcan be subsequently digitised for analysis. The position of amuon ‘hit’ along a scintillator bar can be inferred from the timedi ff erence between signals at either end of the bar. A linear fitto multiple muon hits can be used to reconstruct a muon track.The relative position of each of the bars gives the muon anglein the plane perpendicular to the long sides of the bars. Themuon’s angle in the plane parallel to the long sides of the barscan be calculated using the relative positions of the respectivehits along each bar, inferred from the timing information.From geometric arguments, the angular resolution in theplane of circular borehole faces is approximately 3 ◦ for fivescintillator bar hits. Position resolution based on timing andamplitude information from optical sensors leads to an angularresolution down to 5 ◦ in the plane along the long side of thebars. We consider angular cells of 7 ◦ × ◦ in Section 7.For this study we assume 100% trigger e ffi ciency. The resultsthat we present in Section 7 may be scaled accordingly, oncevalues of trigger e ffi ciency are acquired for a particular detectorconfiguration. We consider a single plane of muon detectorsconfined within a surface area of approximately 1000 m . Thelimiting factor for this choice of surface area is the computation7ime associated with muon transport. ffi ciency As there will be multiple scintillator bars per borehole con-tainer, a muon traversing a borehole container will give riseto multiple bar hits. Background signals, for example due tophoto-sensor noise or radioactivity in the surrounding rock,which are assumed to be short-ranged or highly localised, canbe suppressed by requiring a minimum number of bar hits N hitbar .The e ffi ciency for muons to satisfy the condition on N hitbar isclearly dependent on the number and configuration of scintil-lator bars in the borehole containers. Furthermore, the totalacceptance A for all muons in the detector region must also ac-count for the fraction of the total surface area in the 1000 m detector region that actually contains borehole containers. Wetherefore calculate the detector acceptance by considering twoconfiguration parameters; P det , the fractional surface area oc-cupied by detectors contained within the 1000 m region, and N bar , the number of scintillator bars packed into each boreholecontainer. We consider two values of each parameter which, forboth parameters, we refer to as ‘loose’ and ‘tight’: • P det = { loose ∼ , tight ∼ } , • N bar = { loose : 16 , tight : 24 } .It is assumed that the detectors are arranged homogeneouslythroughout the detector region. In reality, as discussed in Sec-tion 8, the muon detectors will be deployed in a number ofboreholes (approximately 18), which extend radially from amother borehole. Practically speaking therefore, the detectorswill not be arranged in a square arrangement of approximately { × } ≈ , as in our simulation. Instead we en-visage the total instrumented area may occupy up to 1000 m .The choices of P det are simply chosen to show the e ff ect of thereducing / scaling the total number of observed muons from themaximum anticipated yield.The arrangements of scintillator bars in the radial plane of theborehole containers, that we consider for the ‘loose’ and ‘tight’values of N bar , are shown in Figure 5. The acceptance of muonsfor a single plane of detectors, A , which are recorded in thedetector region is shown in Figure 6, for all four combinationsof P det and N bar . The acceptance A is shown as a function of the N hitbar , for surviving muons at the detector site crossing a singleplane of detectors. Realistically, ≥
7. Analysis of muon flux distribution
In this section, an analysis of the muon simulation data isperformed, using muons which have been transported throughthe geological strata in the G eant . It (a) Loose bar packing (16 bars) (b)
Tight bar packing (24 bars)
Figure 5: The arrangement of bars that we consider in this sim-ulation, for the ‘loose’ (a) and ‘tight’ (b) values of N bar . hitbar N ≥ ≥ ≥ ≥ ≥ ≥ ≥ ≥ ≥ A bar , N det Ρ ( (tight, tight)(tight, loose)(loose, tight)(loose, loose) Figure 6: The acceptance of muons which are recorded in thedetector region in one plane of detectors, for all four combina-tions of P det and N bar . The acceptance is shown as a function ofthe number of bars traversed by a muon, N hitbar .is assumed that background noise and signals will be negligibleonce a condition is chosen for N hitbar . We consider increasinglylarge time steps in which muon data is collected, in sequen-tial time periods since first CO injection, as described in Sec-tion 5.2. This utilises the relationship between the radial extentof the CO plume ( r E ∝ √ t ), such that the time intervals cor-respond to periods under which the radius of the CO plume isapproximately constant. We parameterise the significance S ofthe change in the number of muons observed in a given timeinterval as: S = N after − N before √ N after + N before (8)where N before and N after are the number of muon events observedover equal lengths of time, before any subsurface CO is in-jected and after subsurface CO is injected respectively. Thenumber N before is calculated using a statistically independentsample of muons from those used to calculate N after . It is as-sumed that systematic uncertainties cancel in Equation 8. The There is no known e ff ect from electrical noise, radioactivity, or othersources that can lead to an appreciable background over muons at these depths,once the requirement for N hitbar is introduced. N before and N after refer to the number ofmuons detected at any point in the 1000 m detector region. Inreality, the number of muons entering this calculation is reducedprimarily due to the geometric factors discussed in Section 6.The e ff ect of the true values of the P det and N bar is to scale theobservables N before and N after by an e ffi ciency, A . It is thereforeclear that S will be reduced by a factor of √ A . We do notapply factors of √ A to the values of S presented in this section,as realistic values of P det and N bar need to be investigated infurther study. A degree of quantification of the acceptance A can be inferred from Figure 6.The muon flux Φ µ at each time step T , and the associatedsignificance S (measured in standard deviations) of the changein muon count are shown in Table 2. The presented uncertain-ties are taken from the statistical uncertainty in the muon count.A change in the global muon flux after 49 days corresponds toapproximately 24 Gaussian standard deviations. This is clearlysignificant even for an acceptance of A = S is shown in Figures 7 and 8 forthe periods 0–169 days and 169–441 days after first CO injec-tion, respectively. The shape of the plume distributions is vis-ible above the background of statistical fluctuations, even after25 days. Figure 9 shows the distribution of S as a function of θ only. The value of θ corresponding to the maximum value of S is shown to shift as a function of time, which is related to thedynamics of the CO plume and also the sensitivity to the muonpath length.Table 2: The muon flux Φ µ as a function of time step T (orequivalently the radial extent of the CO plume, R ) and the as-sociated significance S relative to the muon flux prior to injec-tion. The presented uncertainties are taken from the statisticaluncertainty in the muon count. T [days] R [arb. units] Φ µ [10 − cm − s − ] S . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± .
8. Discussion
In order to perform this simulation, we have chosen param-eters which are applicable to the geology in the vicinity ofBoulby Mine, Cleveland, United Kingdom. In realistic o ff -shore sites several parameters will change, and di ff erent anal-ysis strategies may be adopted. This section presents a short (a) (b) (c) (d) (e) (f) Figure 7: The significance S of the change in muon count inthe detector plane due to the constant injection of subsurfaceCO , using the model described in Section 4, as a function ofmuon zenith and azimuthal angle θ and ψ . The caption of eachfigure indicates the time interval after first CO injection. Eachdistribution contains 100 values of S , calculated in bins corre-sponding to ∆ θ = ◦ radially and ∆ ψ = ◦ about the radialaxis. A linear interpolation is used between adjacent bins.discussion in order to provide a platform for future studies. Fur-thermore, we present a brief discussion on the practical viabil-ity of this technology, in particular with respect to the relativecosting compared to existing technologies.For this study we have chosen a conservative value for themass-rate of CO injection, of 20 kg / s ( ∼ .
63 Mtpa). Thereforeit is clear that for realistic implementations of this technology,the sensitivity to both the plume extent and density variationswill increase.Clearly the depth of the CO storage reservoir will vary be-tween sites, and therefore the depth of the muon detector ar-ray will also vary between sites. Realistic o ff -shore sites, suchas the Hewett Fields complex in the North Sea, have reservoirdepths (cid:38) . ∼ . S by a factor of ∼ .
39. This does not significantly alter9 a) (b) (c) (d) Figure 8: The significance S of the change in muon count inthe detector plane due to the constant injection of subsurfaceCO , using the model described in Section 4, as a function ofmuon zenith and azimuthal angle θ and ψ . The caption of eachfigure indicates the time interval after first CO injection. Eachdistribution contains 100 values of S , calculated in bins corre-sponding to ∆ θ = ◦ radially and ∆ ψ = ◦ about the radialaxis. A linear interpolation is used between adjacent bins.the conclusions of this study.After significantly longer storage times the lateral plume sizewill extend well beyond what is discussed in this paper. Expe-rience from the Sleipner CCS project [4], which operates in ano ff shore deep saline aquifer (depth ∼ . plumes may have a lateral range of approximately 15 kmafter 14 years. In order to image a plume with such a lateralrange, one needs several arrays of muon detectors positioned atdi ff erent sites relative to the injection point. Several arrays ofdetectors are also required to reconstruct a proper tomographicimage of the CO migration (in addition to the ‘line-of-sight re-construction). This will allow for a triangulation of the muondata, in order to reconstruct a spatially three-dimensional imageof the CO plume. With the detection strategy that is presentedin this paper, which only uses a single detector array, one maystill interpret vertical movements from unexpected changes inthe angular distribution. The manner in which this process isoptimised will vary from site to site.One limiting factor, in terms of the sensitivity of this tech-nique, is the choice of packing fraction P det of the boreholecontainers across the total surface area of the detector array. Re-alistic values for this parameter are currently unclear, althoughthis may not be important if a smaller packing fraction can becompensated by a larger area for detector deployment. [deg] θ Muon
Significance [std. dev.] (a) [deg] θ Muon
Significance [std. dev.]
121 169 days169 225 days225 289 days289 361 days361 441 days (b)
121 – 441 days
Figure 9: The significance S of the change in muon count in thedetector plane due to the constant injection of subsurface CO ,using the model described in Section 4, as a function of muonzenith angle θ . The caption of each of the two figure indicatesthe time interval after first CO injection.As the e ff ectiveness of muon radiography for the mapping oflarge scale objects is already proven (Ref. [12–14]), the utilityof this technology will ultimately depend upon the practical-ity of deployment of the technology at CO injection sites andthe cost of deployment relative to other monitoring technologywhich could be used. Clearly the practicality of deploymentwill require further investigation. We can however provide arough estimate of the relative cost of muon radiography, withrespect to seismic technology, based on a review of freely avail-able information and without reference to particular vendors.We choose to compare to seismic surveying, as it is clearly themost widely used monitoring technology for o ff -shore injectionsites. The relative evaluations are based on a 2015 cost base.Before we present our cost evaluation, we emphasise a num-bers of points: • Costs in this industry are volatile although relative costs,as we will discuss, are approximately fixed. Precise costestimates are impossible on this basis. • Multilateral sidetrack technology, which is discussed inthis section, may be used for deploying the detectors, al-10hough further research may yield more favourable or cost-e ff ective methods. • In the early stages of this technology we suggest to usemodified injection wells, rather than drilling new wells formonitoring. • Muon detector technology is expected to drop signifi-cantly in price if large-scale deployment is implemented,although the scale of this reduction is di ffi cult to quantify.In the few sites which have been developed for o ff -shoreCO storage [4], seismic surveys have been performed aheadof any injection and subsequently on an annual basis. The costof such surveys varies widely and is responsive to the oil price,since most seismic data will be acquired to enable mapping ofpetroleum accumulations and their response to production. Atypical single, modestly sized survey of around 100 km willtoday cost in the region of $5 million. For a storage site last-ing 40 years this equates to $200 million over the duration ofthe injection phase. This will increase further if one considerspost-injection monitoring, processing and reprocessing of theannual data.The cost to deploy muon sensors is similarly di ffi cult to as-sess due to the same volatility issues on oil price. Using thesame cost base we can estimate the incremental cost for deep-ening a well and deploying muon detectors beneath the injec-tion site. A mid-price semi-submersible or jack-up rig for theNorth Sea will today cost in the region of $300,000 per day.For the initial deployment of detectors the mobilisation and de-mobilisation costs will fall to the injection well itself and canbe discounted for these comparison purposes. Detectors wouldbe deployed in multilateral sidetracks drilled from the motherborehole beneath the injection site, in a process known as coiledtubing drilling [see, for instance, 36, and references therein].Coiled tubing drilling is used extensively in o ff -shore settings toincrease the drainage surface of well-bores into an oil reservoir.A mother borehole may have several tens of daughter bore-holes. For the purpose of muon radiography, it is likely thatshort radii sidetracks would be drilled. We could consider 18sidetracks (which allows for a total instrumented surface area of1000 m ) which are each constructed over two days, plus twoweeks to deploy and complete the instrumented part of the wellarray. This equates to about 50 days at $300,000 or $15 millionin total.One must then consider the cost for developing a set of de-tectors that would instrument the required 1000 m area. Weassume that a single detector occupies a surface area of 0.2 m ,and contains approximately 25 scintillator bars, each 1 m long.With current (2015) market costs within a laboratory, one couldacquire all scintillator bars (in one detector) for $4,000, photo-sensors for $6,000, a specially designed data acquisition systemfor $2,000, other components for $1,000 and a further $7,000for technical labour costs. One could envisage a cost reductionfactor of at least two due to mass production - which would im-ply a cost of about $10,000 per detector. Therefore in order toinstrument the required 1000 m area, this implies a cost of $50million. Again, the uncertainties in the precise costing at this stage are clearly very large so this figure may rise or fall, butan estimate of the order of tens of millions of dollars is applica-ble. We are confident that significant cost savings will be possi-ble by using novel and cheaper alternatives to various detectorcomponents, and a further cost reduction by mass productionof detectors, when this technology becomes widely used on anindustrial scale for geological repositories and CCS.Based on experience in experimental particle physics exper-iments, one can expect scintillator detectors to operate for atleast ten years, necessitating approximately four sets of detec-tors over the 40 year injection period. In practice though, seis-mic would still be deployed to detect plume extent in the eventthat the plume can be seen to be moving using muon radiogra-phy. One should note that the cost of both seismic acquisitionand drilling are highly volatile and both are currently falling inline with a low oil price, however the ratio between drilling andseismic costs is e ff ectively stable.
9. Conclusion
We have shown the results of the first detailed simulation ofmuon radiography for the purpose of the monitoring of sub-surface CO stored in geological CCS repositories. Our modelis the first to incorporate geological data to describe a subsur-face CO storage site, with a simplified but nevertheless realis-tic modelling of CO plume evolution and muon detector con-figuration. Non-ideal behaviour in real situations is to be ex-pected because of geological heterogeneity or other factors, forexample, variable wettability of mineral surfaces. Simulationslike this will become an integral part of any future use of thisCCS technology, to compare actual outcomes with the initialreservoir model.We use a detailed numerical model of the fluid flow of sub-surface CO to show that muon radiography is sensitive to thetime-evolution of a CO plume formation. We have shown thatafter approximately 49 days of constant CO injection, muonradiography has a strong sensitivity to changes in column den-sity of (cid:46)
10. Acknowledgements
This work was supported by the Department of Energy andClimate Change and Premier Oil plc. We would like to thankthe Science and Technology Facilities Council (STFC, UK) forsupporting this work (grant ST / K001841 / / L502492 / References [1] N. Nakicenovic and N. Swart,
Emissions Scenarios , IPCC (2000), 570.[2] O. Bert Metz et al.,
Carbon Dioxide Capture and Storage , IPCC (2005),431.[3] K. Stephenne,
Start-up of World’s First Commercial Post-combustionCoal Fired CCS Project: Contribution of Shell Cansolv to SaskPowerBoundary Dam ICCS Project , Energy Procedia (2014), 6106–6110,12th International Conference on Greenhouse Gas Control Technologies,GHGT-12.[4] O. Eiken et al., Lessons learned from 14 years of CCS operations:Sleipner, In Salah and Snøhvit , Energy Procedia (2011), 5541–5548,10th International Conference on Greenhouse Gas Control Technologies.[5] B. Bailey et al., Water control , Oilfield Review (2000), 30–51.[6] J. Senior et al., Geological storage of carbon dioxide: an emergingopportunity , Petroleum Geology Conference series (2010), 1165–1169.[7] V. A. Kudryavtsev et al., Monitoring subsurface CO2 emplacement andsecurity of storage using muon tomography , International Journal ofGreenhouse Gas Control (2012), 21–24.[8] T. Blackwell and V. Kudryavtsev, Development of a 3D muondisappearance algorithm for muon scattering tomography , JINST(2015).[9] K. N. Borozdin et al.,
Surveillance: Radiographic imaging withcosmic-ray muons , Nature (2003), 277–277.[10] W. C. Priedhorsky et al.,
Detection of high-Z objects using multiplescattering of cosmic ray muons , Review of Scientific Instruments (2003), 4294–4297.[11] M. Furlan et al., Application of muon tomography to detect radioactivesources hidden in scrap metal containers , Advancements in NuclearInstrumentation Measurement Methods and their Applications(ANIMMA) (2013), 1–7.[12] L. W. Alvarez et al.,
Search for Hidden Chambers in the Pyramids ,Science (1970), 832–839.[13] J. Marteau et al.,
Muons tomography applied to geosciences andvolcanology , Nuclear Instruments and Methods in Physics ResearchSection A: Accelerators, Spectrometers, Detectors and AssociatedEquipment (2012), 23–28, New Developments in PhotodetectionNDIP11.[14] N. Kanetada,
Geo-tomographic Observation of Inner-structure ofVolcano with Cosmic-ray Muons , Journal of Geography (ChigakuZasshi) (1995), 998–1007.[15] U. Desideri
Advanced absorption processes and technology for carbondioxide (CO2) capture in power plants , vol. 1 of
Woodhead PublishingSeries in Energy . Woodhead Publishing. 2010.[16] P. Mathieu
Oxyfuel combustion systems and technology for carbondioxide (CO2) capture in power plants , vol. 1 of
Woodhead PublishingSeries in Energy . Woodhead Publishing. 2010.[17] S. Kluiters et al.
Advanced oxygen production systems for power plantswith integrated carbon dioxide (CO2) capture , vol. 1 of
WoodheadPublishing Series in Energy . Woodhead Publishing. 2010.[18] M. Blunt et al.,
Carbon dioxide in enhanced oil recovery , EnergyConversion and Management (1993), 1197–1204, Proceedings of theInternational Energy Agency Carbon Dioxide Disposal Symposium.[19] A. R. Awan et al., A Survey of North Sea Enhanced-Oil-RecoveryProjects Initiated During the Years 1975 to 2005 , Society of PetroleumEngineers (2015).[20] D. Lincoln et al.,
Coupled two-component flow in deformable fracturedporous media with application to modelling geological carbon storage ,Numerical Methods in Geotechnical Engineering (2014), 989–994.[21] R. W. Lewis and B. A. Schrefler , The Finite Element Method in theStatic and Dynamic Deformation and Consolidation of Porous Media .John Wiley & Sons, second edition. 1998.[22] A. Murphy and S. Paling,
The Boulby Mine Underground ScienceFacility: The Search for Dark Matter, and Beyond , Nuclear PhysicsNews (2012), 19–24, http://dx.doi.org/10.1080/10619127.2011.629920 . [23] I. H. Bell et al., Pure and Pseudo-pure Fluid Thermophysical PropertyEvaluation and the Open-Source Thermophysical Property LibraryCoolProp. , Industrial & engineering chemistry research (2014),2498–2508.[24] A. M. Rowe and J. C. S. Chou, Pressure-Volume-Temperature-Concentration Relation of Aqueous NaClSolutions , Journal of Chemical and Engineering Data (1970), 61–66.[25] S. A. Mathias et al., On relative permeability data uncertainty and CO injectivity estimation for brine aquifers , International Journal ofGreenhouse Gas Control (2013), 200–212.[26] S. A. Mathias et al., Role of partial miscibility on pressure buildup due toconstant rate injection of CO into closed and open brine aquifers ,Water Resources Research (2011).[27] S. Agostinelli et al., Geant4 - a simulation toolkit , Nucl.Instrum.Meth (2003), 250–303.[28] O. C. Zienkiewicz et al.
The Finite Element Method Fifth edition Volume1: The Basis , vol. 1. Butterworth-Heinemann, Oxford, 6 edition. 2005.[29] M. T. van Genuchten,
A Closed-form Equation for Predicting theHydraulic Conductivity of Unsaturated Soils , Soil Science Society ofAmerica Journal (1980), 892–898.[30] K.-J. Bathe Finite Element Procedures . Prentice Hall. 1995.[31] J. M. Nordbotten et al.,
Injection and Storage of CO2 in Deep SalineAquifers: Analytical Solution for CO2 Plume Evolution DuringInjection , Transport in Porous Media (2005), 339–360.[32] T. Gaisser , Cosmic Rays and Particle Physics . Cambridge UniversityPress. 1990.[33] L. Bonechi et al., Development of the ADAMO detector: test with cosmicrays at di ff erent zenith angles , 29th International Cosmic RayConference (2005), 101–104.[34] V. Kudryavtsev, Muon simulation codes MUSIC and MUSUN forunderground physics , Computer Physics Communications (2009),339–346.[35] P. Antonioli et al.,
A three-dimensional code for muon propagationthrough the rock: MUSIC , Astroparticle Physics (1997), 357–368.[36] A. Afghoul et al., Coiled tubing: the next generation. , Oilfield Review(2004), 38–57., Oilfield Review(2004), 38–57.