Experimental investigation and hybrid numerical analytical hydraulic mechanical simulation of supercritical CO2 flowing through a natural fracture in caprock.

7 Recent work on analogue CO 2 storage sites has shown that the single most defining natural 8 factor as to whether the store can be successfully utilised to retain 99% of the injected CO 2 9 for 1000 years is the behaviour of fractures within the low permeability strata acting as 10 caprock. Here we present experimental and numerical investigation of the hydro-mechanical 11 behaviour of a natural fracture in a caprock during the flow of supercritical CO 2 through it. 12 The caprock is a naturally fractured dolomitic limestone sample recovered from a depth of 13 ~1500m, and is the primary seal to the natural CO 2 storage analogue, the Fizzy field, in the 14 Southern North Sea. For the first time the hydro mechanical behaviour of the fracture is 15 examined using unique experimental equipment applying multiple high pressure single 16 phase supercritical CO 2 fluid flow experiments at representative in situ reservoir pressures 17 (10 MPa to 30 MPa, with confining pressures from 35 MPa to 55 MPa) and a temperature of 18 40°C. The fracture surfaces are scanned to provide high resolution images both prior to and 19 after the experimental investigation. The results are modelled through the further 20 development of a hybrid numerical analytical approach to fluid flow through a discrete 21 fracture, implemented in the open source code OpenGeoSys. The work indicates that 22 through the statistical approximation of the fracture surface and combination of the 23 application of standard nonlinear flow models and analytical mechanical solutions, the key 24 features of the hydro-mechanical behaviour of the supercritical fluid flow through the sample 25 can be replicated. The results provide insight into erroneous effective stress assumptions at 26 higher fluid pressures and the importance of understanding the coupled process multi-27 physics behaviour of fractures in a CO 2 storage setting. Over-simplistic approximations 28 using the effective stress law lead to a Biot’s coefficient greater than 1 being predicted under 29 varying fluid and confining pressures.


Introduction
The long term secure storage of CO2 in the subsurface is one of the key elements of Carbon Capture and Storage, the only industrial scale technology currently planned to reduce CO2 emissions to the atmosphere from fossil fuelled power plants and large industrial point sources.In order for this technology to contribute to climate mitigation efforts, it is crucial that there is no significant migration of CO2 out of the storage site for a minimum of 10,000 years (IPCC, 2005).Storage is most likely to take place in sedimentary basins.The importance of fractures and faults for fluid migration within such basins has long been recognized e.g.Bjorlykke, (1993).Significant effort has been focused on identifying risks related to fractures and faults for CO2 storage (e.g.Rutqvist, 2012, and references therein).Work on natural CO2 reservoirs as analogues for anthropogenic CO2 storage sites has shown the importance of faults and fractures for CO2 migration from the subsurface to the surface (e.g.Shipton et al., 2005, Kampmann et al., 2012).Recent work on natural analogues worldwide has confirmed fractures and faults in the low-permeability rocks overlying CO2 reservoirs as the single most defining natural factor as to whether CO2 can be successfully retained or leakage of CO2 out of the reservoir occurs (Miocic et al, 2013), outranked only by abandoned boreholes.
Obtaining samples of the naturally occurring fractures in caprocks of potential storage sites for experimental investigation is challenging.Caprocks are seldom cored, natural fractures are rare and when fractures are encountered the sample will often fall apart during recovery.
Substantial effort was invested into searching through well records and core logs to identify and then sample a natural fracture occurring in dolomite of the Permian Zechstein formation in the Wissey field, southern North Sea (well 53/04a-9).The Wissey field is adjacent to the CO2 rich gas field, the Fizzy field, which has held a ~50% CO2 concentration in the Rotleigend formation over long periods of time without any indication of leakage (Wilkinson et al., 2009, Yielding et al., 2011).
The Permian Zechstein formation forms the regional caprock for the Permian Rotliegend formation, the main reservoir for gas fields in the Southern North Sea (Glennie, 1998).This caprock represents the seal for many potential anthropogenic CO2 storage sites.The Zechstein formation comprises evaporitic cycles in a landlocked depression which was periodically catastrophically flooded by saline sea waters and in which shales, carbonates, anhydrites and salts up to several 100s of meters thickness were deposited (Ziegler, 1990, Legler & Schneider, 2008).
As the natural fractures within the caprock control the efficiency of the storage system, it is important to understand the coupling of the rock stresses confining fractures and the fluid pressures within fractures on the flow properties of the fractures in a CO2 storage situation.
A number of authors have experimentally investigated the relationship between the confining pressure and the fracture permeability, usually with limited (<2 MPa) fluid pressures, e.g.Yasuhara et al. 2006, Walsh et al. 2008 and several references therein.However there is very little experimental work on the flow of supercritical fluids through quasi-tight natural fractures under reservoir conditions of stress, fluid pressure and temperature.
Experimental investigation of such hydro-mechanical behaviour requires state of the art equipment capable of; containing supercritical CO2;, measuring differential flow pressures of a few thousands of Pascal where the downstream pressure is of the order of 10 to 30 MPa, and maintaining constant but adjustable temperature conditions, as demonstrated in Edlmann et al. (2013).
Simulation of the results to improve understanding of the key processes operating requires state of the art multi-physics coupled process simulators.Several authors have published on the coupling of hydraulic and mechanical responses, and appropriate solution techniques, a key work being Lewis and Schrefler, 1998.Ignoring the coupling of these processes can lead to a significant oversimplification of the system, and approaches which do not adequately represent the systems being modelled (e.g.McDermott et al. 2006).However incorporation of coupling processes and associated phenomena at different scales leads to numerical stability issues related to spatial distribution, time and parametrical heterogeneity (e.g.Ouyang and Tamma, 1996;Kolditz, 1997, Wang et al., 2011a;Wang et al., 2011b).
Replicating, understanding and thereby having the ability to predictively model the behaviour of these systems is dependent upon identifying the key processes operating under a given set of conditions.Numerically, and using modern computational codes, it is possible to effectively overkill the solution of coupled process problem, whereby the key processes are missed, a large amount of computational time is required and several fitting variables are included.To simplify computational demands it is possible to integrate analytical and physical models into standard numerical modelling techniques.
In order to simulate the multi-physics coupling of the hydraulic and mechanical behaviour of flow through a porous or fractured media, usually some form of the effective stress law is invoked, whereby the fluid pressure within the porous media or fracture is considered to work against the confining pressure.The coupling between the fluid pressure and the mechanical stress is moderated by the Biot's coefficient, representing the efficiency of the transfer of the fluid pressure into the mechanical skeletal stress.A Biot's coefficient of greater than 1 would indicate that the fluid pressure has an overall influence greater than its pressure can provide which is not physically possible, and therefore of some significant interest.Such apparent behaviour has been observed experimentally in porous media, e.g.Zoback andByerlee, 1975 andMoghadam et al. 2014.In this paper we present the results of new experimental methods to investigate the flow of supercritical CO2 through fractured rock samples under varying fluid pressure and confining pressures representative of natural CO2 storage reservoir conditions.We investigate one of the few naturally fractured caprock samples recovered from a depth of ~1500m from the North Sea and use fracture scanning techniques to provide the geometry of the fracture surface.We then simulate the experimental results using a hybrid numerical analytical approach developed in open source software OpenGeoSys.Coupled physical models for fracture aperture closure due to pressure changes, turbulent flow, elastic and plastic deformation were implemented and are described here.
The main features of the experimental data could be replicated with parameterisation well within the expected literature values.The effect of an apparent Biot's coefficient of greater than 1 can be reconciled to processes of compression of the fracture matrix, highlighted due to the high fluid pressures required for the experimental investigation.

Experimental Investigation
A well core sample containing a natural fracture within the Zechstein dolomite of the Southern North Sea field, Wissey, was sampled from well 53/04a-9.A 38mm diameter cylindrical sample of the natural fracture was cored at the University of Edinburgh, and the end surfaces were smoothed using a lapping plate to ensure they were exactly perpendicular to the sample axis.The fracture surfaces of this prepared sample are shown in Figure 1.In-situ depth 5525' / 1684 m In-situ pore pressure 18.5 MPa (Noy et al, 2012) In-situ temperature ~47°C (Harper, 1971) In-situ vertical stress 38 MPa (Noy et al, 2012)

Experimental rig description
The experimental rig of Edlmann et al. (2013)  The fractured Wissey sample is held vertically within a Hassler-type uniaxial coreholder contained within an oven to provide temperature control and minimise external laboratory temperature fluctuations that could influence the hydraulically controlled confining pressure.
The sample is contained within an elastomer sleeve and an even external radial pressure is applied to the sample by pressurising confining oil on the outside of the elastomer sleeve using a high pressure hydraulic hand pump.Fluid ports are located on both the upstream and downstream end platens which allow CO2 to be pumped into the bottom of the sample and out of the top to minimise any unwanted gravitational effects.Steady state scCO2 flow is obtained through the use of paired high pressure syringe pumps (Teledyne ISCO 100DX in Hastelloy C-276 for resistance to CO2) to ensure accurate steady control of both the CO2 flow rate (from upstream) and the downstream fluid pressure.
To ensure fluid bypass around the outside of the core sample does not occur during experiments, the fluid pressures is kept significantly below the radial confining pressure acting on the sample at all times.In order to facilitate the filling of the syringe pumps with liquid CO2, the syringe pumps must be cooled to 2 o C therefore it is necessary to heat up the CO2 prior to it entering the fractured sample.The CO2 is pumped through a heat exchanger set to the experimental temperature before the CO2 pipework enters the core holder oven, heating tapes on all pipework ensure the steady temperature is maintained.Due to adiabatic expansion of carbon dioxide released from the downstream syringe pumps during their emptying process, and resultant Joule-Thomson cooling effects, there is a danger of the disposal pipework downstream of the pumps freezing and blockages occurring.In order to minimise this risk, an automatic back pressure regulator is fitted downstream of the syringe pumps to step down the pressure in the syringe pumps to 5 MPa (725 psi) before release to the atmosphere.Downstream of the automatic back pressure regulator heating tape is used to prevent the risk of pipe freezing occurring as the CO2 expands on release to atmosphere.For successful operation of the experimental rig, the integrity of component materials in the presence of scCO2 is important and all metal parts are of 316 stainless streel, however there is a corrosion risk and all valves (especially needle) needed to be inspected regularly.CO2 diffusion was observed through the traditional Viton core holder sleeves so HNBR sleeves were used.To further ensure there was no CO2 diffusion through the core holder sleeve the core sample was coated prior to loading into the cell.The coating comprised a layer of PTFE tape to hold the fractured core sample together then a layer of aluminium foil finished with a lead liner jacket.The ends of the sample were left uncovered to allow CO2 flow through the sample fracture within the coreholder.
Continuous (5s frequency) measurement and recording of confining pressure; fluid pressure both upstream and downstream of the sample; and sample temperature is carried out during experiments.In addition, a high accuracy differential pressure transducer is used to continuously measure and record differential fluid pressure across the sample.Table 2 provides details of the measurement instruments and accuracies.

Experimental program
The hydraulic mechanical response of the Wissey natural fracture to flow of supercritical Within the experiments presented here, fluid pressures were varied at each change in confining pressure, then the confining pressure was stepped up to 35,and fluid pressures varied again, this also continued at confining pressures of 45 and 55 MPa.For each static confining pressure, flow rates within the range 1-10 ml/min were tested as shown in Figure 5.
The whole sequence of experiments (A to E, Figure 4) lasted more than 1 month, the experimental results considered here are typical of the whole profile, and represent conditions after several cycles of loading and unloading.This paper demonstrates proof of concept to the experimental work, the numerical modelling and some significant insights into the processes operating.Although different temperatures were investigated in the experimental sequence, the results presented and modelled in this paper are at constant temperature of ~40 C. At this stage only a selection of the several thousand data points have been analysed using this modelling approach.During loading and unloading of the fracture, hysteresis was observed, as is commonly observed in rock mechanical tests.With increasing cycles of loading and unloading the response of the system becomes more elastic, the hysteresis reduces significantly.The area chosen was taken as being a representative selection of the different pressure scenarios during which investigation had taken place, and after initial loading and unloading cycles where the elastic effects could be examined and the plastic deformation effects were limited.They demonstrate proof of concept of the numerical approach and the conceptual understanding of the processes operating during the experimental investigation.An estimate of effective stress during the experiment,   , acting on the sample may be calculated using the effective stress law:   =   −   , where   is the confining pressure and   is the average fluid pressure (Terzaghi, 1943).This estimate is shown for each experimental step in Figure 6.(2) In practice under the experimental conditions we examined there was minimal difference between the compressible and incompressible flow estimates of transmissivity.
Given the extremely low porosity of the sample matrix, flow through the sample can be considered to occur within the fracture only, with negligible transport through the sample matrix.Thus the cubic law can be used to estimate a hydraulic aperture for the fracture.
The hydraulic aperture will not be equivalent to the mechanical aperture of the fracture (Renshaw, 1995) due to flow losses as a result of tortuosity and surface roughness within the rough rock fracture.The mismatch between these apertures is likely to increase as the aperture decreases.However, the hydraulic aperture can be considered to be indicative of the mechanical aperture and is therefore an extremely useful parameter in evaluating fracture closure and deformation during the flow experiments.
According to the cubic law, the intrinsic fracture permeability can be estimated as below (Witherspoon et al, 1980): Thus, using estimates of the parameters above (T, eh, kfrac) derived from experimental observations (∆P), we can analyse the fracture response to the differing flow, confining pressure and back pressure conditions considered during the experiments.

Surface Scan Data
Very seldom is reliable discrete fracture aperture data available for simulation of flow through a fracture surface.However, it is possible to scan fracture surfaces to provide information on the asperities of the fracture surfaces, which can then be recombined numerically to provide an estimate of the fracture aperture.
A Micro-Epsilon ScanCONTROL high resolution laser scanner (model: LLT2700-100), at the University of Strathclyde was used for 3D profiling of the experimental sample fracture surfaces both before and after the experimental investigation.Use of the Micro-Epsilon scanner enabled fast, non-contact, scanning of the experimental samples with elevations recorded to a resolution of 5 µm.The scanner is mounted on a mechanical arm that moves the scanner along the y-axis.This enables continuous scanning along the length of the sample.The full width of the surface can be captured within a single laser scan profile.Thus, the full fracture surface can be captured through multiple scan profiles taken along the y axis.
Using this technique the fracture surface profile was scanned with an xy resolution of 200 µm x 40 µm.This provided per surface circa 250,000 data points comprising an x,y,z coordinate.Scanned fracture surface data of a fracture surface is illustrated in Figure 7

Statistical modelling of the fracture surface
Examination of both Figure 1 (the photographic image of the fracture surface) and Figure 7 (scanned topography of the fracture surface) shows that there are two scales of surface profile roughness, the aperture "waviness" at a scale of 2 mm or greater providing the topography of the surface, and the small scale roughness at a sub-millimetre scale.Previous authors have addressed the issues of different scale of roughness, and aperture matching, e.g. Brown et al. 1986, Brown 1995, Glover et al., 1999, Walsh et al., 2008, Zou et al. 2015.At the scale of the experimental investigation, the flow through the fracture is controlled by the small scale roughness of the fracture surface and the mismatch of the upper and lower surfaces.The larger scale "waviness" is not so relevant to the experimental data where the sample is of the order of 4 cm long.In the field the larger scale "waviness" would have an impact, see for example Zou et al. 2015, but as there is only one wave length at this scale in the sample (see Figure 8), the effect is can not be separated in our experimental set up.The fracture scan data is determined as x,y,z data, therefore it is necessary to remove the influence of the larger scale waviness to be able to assess the characteristics of the smaller scale roughness.To do this a low resolution asperity reference surface (LRARS) was calculated by generating a linearly interpolated grid at a significantly coarser resolution than the detailed surface roughness (see Figure 8).This grid, however had a much higher resolution than the "waviness".The LRARS has a resolution of 2 mm in the x and y direction, the z coordinate was generated as the best linear interpolation fit of the experimental data.
From the LRARS a high resolution asperity reference surface (HRARS) was generated as the resolution of the experimental data, the information of the LRARS being used to predict the surface profile of the HRARS.The HRARS was then used as a reference surface to evaluate the aperture profile from and determine the statistical distribution of the small scale asperity variation.This surface asperity data is presented in Figure 9 To generate the aperture data from the asperity data, two normal distributions are fitted to the asperity data (see Figure 9, Distribution fit 1 & 2), mathematically described by a normal distribution as.

Numerical Simulation
We The fracture plane is discretized into individual elements, and an aperture is mapped to each element, this concept is illustrated in Figure 12.The spatial density of the mesh was chosen at 0.5mm x 0.5mm, creating 3933 nodes and 3807 elements.This mesh density was a compromise between the computational expense, the resources available and based on previous investigation on the optimal resolution needed to capture material behaviour of fracture surface at a similar scale (Bond et al. 2014).
The permeability of each element representing part of the fracture surface is assumed to be given by the cubic law approximation, e.g.(Koyama et al. 2006;Kalbacher et al. 2007) where e is the aperture, and e k  represents the assigned element permeability in coordinate directions  and  .
where S is the storativity coefficient [Pa -1 ], k denotes the permeability tensor [ 2 m ] derivation of which is described above, s rel k is a relative permeability correction described to account for nonlinear flow, P is the fluid pressure in Pa, and Q is source/sink [s -1 ].This equation is valid for a saturated, non-deforming porous medium with heterogeneous hydraulic conductivity.To evaluate the fluid pressure distribution within the fracture a numerical solution is required.The solution of equation ( 8) using the finite element technique is covered in standard works such as Istok, 1989;Lewis andSchrefler, 1998 andZienkiewicz andTaylor, 2005.A steady state solution of ( 8) is selected, with the reason that the fluid volume change as a result of the change in aperture of the fracture with time is minimal compared to the total amount of fluid flow through the fracture during the experiment.The elastic and plastic response of the fracture aperture is simulated using analytical and empirical solutions discussed below.As the contact area in the fracture changes dynamically with time (plastic response) and with the fluid and skeleton stress (elastic), it is necessary to solve the mechanics using an iterative solution to deal with the deal with this nonlinear coupling.
Iteration 1 relates to the closure of the fracture under confining pressure.Once the closure of the fracture has been evaluated this is related nonlinearly to the fluid pressure, which is a function of the permeability in the fracture, related to the fracture aperture.It is therefore necessary to address this second nonlinear coupling using a further iterative approach, Iteration 2. Iteration 1 uses a standard bisecting intervals method approach, more details for using this method in fracture closure are found in McDermott 2015, iteration 2 uses a standard Newton Raphson iterative approach.

Figure 14 Numerical solution procedure
In (8) flow is assumed to be laminar and incompressible.However supercritical CO2 has both liquid and gas properties, i.e. the density of a liquid, but the viscosity of a gas.It is also known that scCO2 is compressible, and has pressure and temperature dependent density and viscosity functions.As a first approximation we use a nonlinear flow correction described in Kolditz, (1997) and simulate the fluid with constant viscosity and density functions.At this stage the complexity of including temperature and density dependent fluid properties, as well as extra iterative functions required to simulate compressible flow seemed beyond the accuracy of the data available for the fracture surface profile given we apply a generic statistical representation of the surface, and beyond the measurement accuracy of the experimental work.However this is a step which could be added at a later date in a further development if the accuracy was required.Zhang and Nemcik (2013), and references therein, discuss fluid flow regimes and nonlinear flow characteristics in deformable rock fractures.They measure in four fractured sandstone samples nonlinear flow parameters, describe the nonlinear flow using the Forchheimer equation (Forchheimer 1901) and Izbash's approach (Izbash 1931), and show a relationship of the nonlinear flow to the amount of confining stress, acting as a proxy for the degree of closure across the fracture surface.
After Forchheimer, 1901 the inertial effect as a function of the kinetic energy term can be included in the relationship between the pressure gradient and the flow velocity as where hydraulic head is used rather than the pressure gradient, and  is some parameter such that 21   .This approach is closely related to that from Izbashes (1931).This concept is integrated in the current model such that we simulate a linear and a nonlinear flow Where the suffix L P  is the value of the pressure gradient causing linear flow and the extra addition pressure t P  is that part considered to be causing turbulent flow.The model described in this paper has been integrated into the scientific OpenGeoSys code, a standard Galerkin finite element solver, (Kolditz et al., 2012).
Simulation of the elastic response of the fracture to confining pressure and fluid pressure changes.
The contact stress, a  , is defined as the normal confining stress across the fracture minus the fluid pressure u in the fracture carried by the fractional contact area (CA) of the fracture.

 
The analytical elastic response is given by a remarkably simple assumption that the contact areas and channel surfaces act as loading plate in semi elastic space.This is illustrated in Fracture closure or opening occurs as a consequence of compression changes across the fracture at the contact areas in response to changes in the contact stress.We assume that the contacts in the fracture represented by quad elements in our mesh can be represented by cylindrical pillars.In reality we do not know exactly what shape the contacts are, and the cylindrical assumption ensures uniform deformation over the element.The amount of closure is given by a standard elastic solution for the deformation of a ridged circular plate loading an elastic half space.The diameter of the sample is 2 orders of magnitude larger than the assumed area of loading allowing the application of the infinite depth approximation associated with this analytical expression.We chose the edge of the ridged plate as most representative of the loading conditions at the contacts.The deformation z  at the edge of the plate for an infinite body is given by Poulos and Davis, 1974, as here L represents the loading length, being the equivalent radius of the loaded circular plate, E is Young's modulus, v is Poisson's ratio.
Within the channels the same basic solution is applied, in this case the deformation at the centre of the plate is taken, This solution is for single plate loading and does not allow for the superposition of concurrent stress fields.Additionally it relies on an approximation of the geometry of the fracture surface to evaluate an average value for L, both for the channels and for the contact areas.However these uncertainties are within the overall uncertainty of the data, and in this case are most likely accommodated for within the range of possible realistic E modulus values.
Stress corrosion is considered to be the process occurring leading to plastic deformation.
Stress corrosion occurs when during closure of the fracture surface mechanical across the asperities leads to the development of micro fractures in the rock surface and the destruction of asperities.Stress corrosion is discussed in more detail in Zhang et al., 2012.They describe stress corrosion as velocity of crack propagation as a function of the geochemical interaction between the fluid and the rock grains and the impact of the stress across the grain contacts.The impact of supercritical CO2 on the rate of stress corrosion development is not well known.Zhang et al. 2012, derive an expression for the contact stress which can be shown to be identical to equation ( 16).They then derive the stress following Yasuhara and Elsworth 2008 at the edge of the asperities, such that Following Atkinson (1982) and references therein tensile crack velocity can be described as Here v is the crack velocity, 0 v is an experimentally determined value, H  is an activation enthalpy, R is the gas constant and T is the absolute temperature.
1 n K is the stress intensity factor for tensile (mode 1) failure and n is a material constant known as the stress corrosion index.This model is linked with (19) using the approximation that where r is a characteristic crack length and Y is a geometrical constant, as per Atkinson (1982).
Several authors express more detailed chemical dependent formulas to describe the impact of chemical attack and stress intensification on the mineral surfaces and fracture tips enhancing stress corrosion.These approaches include experimentally measured fitting constants, e.g.Dove (1995), and Zhang et al. (2012) and references therein.In the current paper we note that there are so many experimental uncertainties that it is impossible to reproduce exactly the formulation of crack velocity for supercritical CO2 from the literature.
Although a number of experimentally determined constants have been measured for the development of stress corrosion in carbonate systems (e.g.Atkinson,1982, Henry, 1977) we are not aware of the values of experimentally determined factors being measured for supercritical CO2 both in the carbonate or silicate system in such a setting.Therefore to model the plastic deformation caused by stress corrosion at this stage we adopt the overall approach illustrated above, but simplify it such that where f is some factor relating the crack propagation velocity to the geometry of the system and the closure of the fracture.We do not try to break it down further into its individual components as the data is not available to validate this.The parameter s  is chosen to fit the experimental results, and controls the non-elastic behaviour of the fracture surface during repeated loading cycles.

Results and discussion
Flow of supercritical CO2 through a natural fracture in a major caprock lithology found extensively throughout the southern North Sea was investigated under different conditions of flow rate, confining stress and fluid (pore) pressure.The confining pressure was increased and decreased from 35 MPa to 55MPa, supercritical CO2 fluid flow rate was increased and decreased from 1 to 13 ml/min, and supercritical CO2 fluid pressures altered from 10 MPa to 30 MPa simulating a wide range of conditions and enabling the hydro-mechanical response of the fractured sample to be investigated.The experimental results were simulated using a hybrid numerical approach implemented in the open source code OpenGeoSys.
The results of the comparison between the simulated results and the experimental results are presented in Figure 16.The upper section of this graphic illustrates the experimental conditions under which the sample was placed.The lower section of this graphic indicates the hydraulic aperture calculated from the experimental results, and the hydraulic aperture predicted by the modelling approach.The experimental conditions can be seen to comprise three repeated cycles at different frequencies ensuring a wide range of conditions were investigated.The confining pressure is stabilised at three levels, 35, 45 and 55 MPa, the upper black line labelled "Confining stress" in the diagram.Whilst at a constant confining pressure the downstream fluid pressure is altered from 10 to 20 to 30 to 20 to 10 MPa.This cycle is illustrated for 35 and 45 MPa confining pressure and partially for 55 MPa confining pressure, and is seen in the lower black line, labelled "Downstream fluid pressure".For a constant downstream fluid flow pressure, a series of different flow through rates are applied.
Illustrated in green and labelled as "Flow rate" on the right hand axis.
The hydraulic aperture of the fracture under changing conditions of flow, confining pressure and downstream pressure is calculated (3), and illustrated in Figure 16 as "Experimental results".Simulation of the experimental results is presented in this figure labelled as "Model results" and the parameterisation of the model presented in Table 3. Better calibration of the model addressing a number of model sensitivities could potentially improve the fitting of the results, but would not necessarily add to the process understanding, or be more valid than the current set of results.
In the literature there is some discussion about the validity of the cubic law approximation (Witherspoon 1980) which is applied in this paper for both the analysis of the experimental It can be seen that despite the approximations inherent in the application of the cubic law to derive both the experimental and numerical permeability, the elastic model approximations, the use of statistically generated fracture aperture surfaces based on the surface scans, and the approximations inherent in the use of combined analytical and numerical solutions, the results follow the experimental trends remarkably well.Where external data is available to validate the results, the parameterisation in the model are well within the bounds usually expected for the materials and fluid modelled.
The implication of the reasonable fitting of the experimental results is that the key processes and the controls operating during the experimental work have been well understood in the modelling approach.This also suggests that the use of hybrid numerical and analytical modelling techniques forms a viable method for efficient evaluation of coupled process problems, as has been shown by several other authors, e.g.(Moës et al., 1999;Pruess 2005, Qin 2005;McDermott et al. 2007, 2011, Mohammadi, 2008) The four processes included which led to an aperture change include 1. Elastic compression and expansion of the contacting fracture asperities 2. Elastic compression and expansion of the fracture wall between contacting asperities 3. Non-linear flow 4. Plastic deformation during closure The effect of including these processes is illustrated as paired points marked in Figure 16 and discussed below.The difference between these two points is that there has been some loading and subsequent unloading of the fracture surface which has led to a permanent change, or plastic deformation, in the hydraulic aperture.This is currently represented by a stress corrosion term.
The stress corrosion parameter represents the plastic closure of the fracture with time as a consequence of mechanical damage of the fracture surface and asperities.After the experiments had been undertaken, the fracture surfaces were rescanned.It is of note that mechanical damage and a general "flattening" of the small scale roughness of the surfaces could be observed.
All the parameters in A key consequence of the simulation was that it was impossible to fit the variation of the hydraulic aperture with the changes in the system without taking into account the elastic opening and closing of the fracture as a consequence of the channels available to fluid flow elastically expanding due to the high fluid pressure compressing the matrix between contact points.This is a feature which we are not aware of having been observed experimentally before.It was observed here due to the high fluid pressures which were required to maintain the supercritical CO2 conditions.
Using standard approximations of the relationship between confining pressure, effective stress and hydraulic aperture utilises the effective stress law such that Where '  is known as the effective stress across the fracture, here  represents the Biot's coefficient.This can be understood as the efficiency of the transfer of the fluid pressure to work against the normal confining pressure.The effective stress is then used via various analytical or empirical laws to predict the change in fracture permeability as a function of changing fluid pressure.Rutqvist (2015) and Jiang et al. (2010) discuss several depth normal stress permeability models of many of the effective stress permeability relationships used.However the impact of the fluid pressure in compressing the fracture wall between the contact points is not explicitly taken into account, although the form of the empirical fittings used will accommodate the increase in the channel space.
For the experimental results in this paper, if the elastic expansion of the flow channels was not taken into account it was impossible to model the experimental results.Trying to match the experimental results using (24) would require the higher fluid pressures to be associated with a Biot's coefficient of higher than 1 to maintain the required fracture aperture corresponding to the confining pressure.In reality what is occurring is there is still compression at the contacts of the two fracture surfaces, but there is also an increase in the fracture aperture due to the fluid pressure compressing the matrix between two contact points and creating extra space for fluid flow.The implication of a Biot's coefficient greater than 1 is that the fluid is exerting a higher pressure on the rock matrix than the fluid is at.
That is for instance a fluid at 20 MPa exerts 22 MPa on the rock, this is obviously physically incorrect.
The consequence of this observation is that at higher fluid pressures the compressive effect that the fluid pressure has on the matrix needs to be taken into account when simulating changes in fluid pressure.This type of inference has been made before on porous media, see Zoback andByerlee, (1975), andBerryman (1992) examined the possibility of using multiple continua to understand such behaviour.Biot's coefficients greater than one indicates that important processes controlling the flow rate are not being simulated in the model.To simulate the results it was necessary to further develop a hybrid analytical numerical finite element model so that it included elastic deformation, plastic deformation and nonlinear flow.Additionally fracture profiling data both before and after the experimental work was used to develop statistical models of the fracture surface, fit aperture distributions and test damage hypothesis.

Conclusions
The hybrid analytical numerical model was able to match the experimental results with typical parameters extremely well despite widely varying conditions of flow, fluid pressure and stress.The use of simplified effective stress models for the prediction of high pressure fluid flow through rocks can be shown to be inadequate for the in situ fluid and rock stresses under investigation, leading to implicit assumptions that the Biot's coefficient is greater than 1.
The experimental work and the numerical investigation has led to a significant increase in process understanding with respect to the flow of high pressure supercritical fluids through fractures in low permeability media at depth.The model will now be used as a benchmark to investigate the remaining data set presented in this work and several other comparable data sets recently derived from similar experimental investigation on a wide variety of caprocks.

Figure 1
Figure 1 Image of the naturally fractured sample recovered from the Wissey field.
was further developed to enable multiple high pressure, high temperature supercritical CO2 flow experiments to be undertaken through fractured, 38 mm diameter core samples (up to 75 mm in length) under various confining pressures.To assess the thermal and mechanical controls on the hydraulic behaviour of the fractured core sample to scCO2 flow, the rig facilitates control of: radial confining pressure acting on the sample (up to 69MPa), the downstream fluid (pore) pressure (up to 95MPa), steady state CO2 fluid flow rate through the sample (up to 20ml/m) and the sample and fluid temperature (up to 80 o C).Figures2 and 3.The hydraulic, thermal and mechanical behaviour is analysed through the measurement of the radial confining pressure, upstream and downstream fluid pressure, the differential fluid pressure across the sample, flow rate across the sample, and fluid and sample temperatures.

Figure 2 Figure 3
Figure 2 Experimental rig -CO2 flow-through analysis CO2 was tested under downstream fluid pressures of 10, 20 and 30 MPa; and confining pressures of 35, 45, 55 MPa.The sequence of experimental results presented and modelled in this paper are highlighted in purple in Figure 4.It can be seen that these are part of a larger sequence of experimental results which are under ongoing simulation.The experimental step noted on the abscissa represents changes in the fluid flow rate, the measurements being taken once steady state conditions were observed.

Figure 4 Figure 5
Figure 4 Pressure scenarios tested and associated sample (and fluid) temperatures during flow experiments (purple band = results modelled)

Figure 6
Figure 6 Effective stress during flow experiments (purple band = results modelled here)

Figure 7
Figure 7 Fracture scan data of one surface of the Wissey sample

Figure 8
Figure 8 Asperity reference surface represented as a grid

Figure 9
Figure 9 High resolution asperity surface distribution, and two statistically identical generated profiles fitting the data.

.e
The normal distributions are then fitted together as if each represented one surface of the fracture.The aperture is then evaluated as can be considered to be the average fracture aperture, and m is the mismatch parameter between the two surfaces.The effect of the parameters e0 and m on the simulated aperture distribution are illustrated in Figure10

Figure 10 Figure 11 Figure 11
Figure 10 Combination of two statistical distributions of the asperities on the fracture surface to create an aperture distribution accounting for mismatch of the surface profiles.
further develop a model described inMcDermott et al. 2015 for simulating pressure and chemical solution dependent closure of a discrete fracture whose aperture is given in x,y coordinates.Here we develop that approach to include the elastic response of the matrix and fracture to both the changes in confining pressure and the changes in the fluid pressure within the sample, we also include nonlinear fluid flow and include stress corrosion as a plastic response to confining pressure.To match the experimental results it was necessary to include 1. Elastic compression and expansion of the contacting fracture asperities 2. Elastic compression and expansion of the fracture wall between contacting asperities 3. Non-linear flow4.Plastic deformation during closureThe above processes were included in as simplistic fashion as possible and in response to our understanding of the fracture behaviour.The parameterisation of these processes remained static for the simulation.The implementation of these processes is described in the results section the impact of including the different processes on the simulation results is demonstrated.

Figure 12
Figure 12 Aperture distribution is mapped to a mesh representing the fracture plane, from (McDermott et al 2015).

Figure 13
Figure 13 Conceptual numerical model with discrete apertures P is the fluid pressure in pa,  is the viscosity [Pa s], v is the fluid velocity [m/s]  is the fluid density and  is the turbulence factor.This equation can be seen to show that up to a certain fluid velocity there is a quasi linear relationship between the flow and rate and the pressure gradient.As the flow velocity increases beyond a certain point the kinetic energy term becomes more important and the linear relationship breaks down.
regime divided by a Reynolds number such that flow above the Reynolds number is considered nonlinear.The individual element Reynolds number Re e is calculated from (11) compared with a user specified input parameter Reynolds number Re u above which turbulent flow is considered by the model.The equivalent value of * k given in (10) used in (8) for the simulation with a linear expression,

Figure 15
Figure 15 Analytical solutions used for evaluating the impact of rock stress and fluid pressure changes'Given the amount of contact area in the fracture, the characteristic length L of this is evaluated for the contact area (change of fracture aperture due to contact stress variation) and for the channel area (change of fracture aperture due to fluid pressure variation).
change in aperture due to stress corrosion and represents the aperture change due to mechanical damage and t  the time interval over which the change is being calculated for.In the model we use the value s  as a factor indicative of plastic deformation, from equations (20), and (21) this can be shown to be equivalent to hydraulic aperture and the numerical simulation of the results.The experimental hydraulic aperture derived from the cubic law is an average approximation of the measured fluid flow properties of the fracture.Oren et al. (1998)  show that local roughness features within the fracture surfaces can impact the validity of the cubic law to determine the hydraulic aperture, particularly where fracture surfaces can not be represented as parallel planes within certain geometrical criterial and where the areas of surface contact significantly influence the availability of flow paths in the fracture surface.In this work the simulated contact areas remain below 15%, according to Oren et al. a difference of a few percent due between the Navier Stokes solutions and the cubic law approximations can be expected.Brush and    Thompson (2003)  simulated several fractures surfaces of varying roughness's using average values of asperities on a 2D mesh similar to the numerical methodology presented in this paper and compared complete Navier Stokes simulations against local cubic law approximations.For surfaces with similar characteristics to those of the sample investigated here, they showed that there was a minimal difference between the Navier Stokes approximation and the cubic law approximation.Oren et al. (1998) and Bush et al. (2003)     both indicated that the cubic law approximation could not account properly for higher velocity flow rates where the inertia term of the Navier Stokes equations becomes significant, i.e. the onset of turbulent flow.Zou et al. (2015) examine fracture surfaces in detail and solve the Navier Stokes equations to demonstrate the impact of local scale roughness on the flow through the fracture, particularly the development of eddies due to the detailed structure of the fracture surface.In this work the fracture surfaces were approximated from a statistical representation of the fracture surfaces and not a discrete measurement of the aperture profile.Solving the Navier Stokes Equations in detail would exceed the accuracy of the geometrical data.The impact of the inertia term is accounted for by including a non-linear flow parameter discussed above in (10).

Point pair 1 :
Elastic compression and expansion of the contacting fracture asperities; Moving from left to right the first point 1 has a larger hydraulic aperture than the second point 1.Referring to the experimental conditions, only the confining stress has increased by 10MPa from 35 MPa to 45 MPa.This change is represented in the model as elastic closure or expansion of the fracture asperities.Point pair 2: Elastic compression and expansion of the fracture wall between contacting asperities; Moving from left to right the first point 2 has a smaller hydraulic aperture than the second point 2. Comparing these points the only difference in experimental conditions is that the downstream fluid pressure has increased from 10 MPa to 20 MPa, under a confining stress of 45 MPa.The response at Points 1 previously illustrates the size of the expected elastic rebound due to asperity compression.The response at Points 2 is clearly larger, even although the effective stress change is the same.The majority of the increase in fracture aperture is due to the compression of the fracture wall between the contacting asperities caused by the increased fluid pressure and leading to a larger hydraulic aperture.Point pair 3: Non-linear turbulent flow approximation; Moving from left to right, the only change experimentally between these points, and measurements in between them is that the flow rate in the fracture increases.This has a clearly non-linear impact on the derivation of the hydraulic aperture.The larger the flow rate, the smaller the hydraulic aperture appears to be.This is accounted for in the model by the effect of a turbulent flow approximation.Point pair 4: A plastic response to closure; Moving from left to right, there is no change in terms of experimental conditions of flow rate, downstream pressure and confining pressure.

Figure 16
Figure 16 Simulation of multiple loading cycles of confining pressure, fluid pressure and fluid flow rates.Marked points discussed in text.
Unique experimental results are presented for the flow of supercritical CO2 through a naturally fractured caprock sample from a southern North Sea reservoir.A series of experiments representing in situ conditions of fluid pressure and rock stress were carried out with cyclical loading of the sample, cyclical fluid flow rates and changes to the in situ fluid pressures.The fluid flow rates, confining pressures, and differential pressures were recorded providing an insight into the coupled hydraulic and mechanical behaviour of a fracture through which supercritical CO2 is flowing.The results are presented as hydraulic aperture as a function of time, fluid pressure, confining stress and flow rate.

Table 1
below contains a summary of the sample details.

Table 1
Wissey sample details

Table 2 Measurement instruments and accuracies Parameter Measuring instrument Instrument accuracy
table 3 are within the standard literature range they require little further comment, apart from the Reynolds number.Within the literature there is quite a disparity shown by authors simulating both fractured and porous media exactly where turbulent flow should be taken into account.Zeng and Grigg (2006) and various references therein for water and gas flow in porous media range give Reynolds numbers from 0.11 to ~1000, with an upper value being 2300 from De Marsily, 1986.Zhang and Nemcik 2013 present results for the onset of turbulent water flow in fractured sandstone which suggest values of 5 to 20.
The value of ~0.5 given in the table seems fully reasonable for supercritical CO2 flow in this experimental work.