E ﬀ ect of cold CO 2 injection on fracture apertures and growth

The injection of cold CO 2 is modelled in three dimensions using a two-stage coupled thermoporoelastic model, with the aim of evaluating changes in apertures and potential growth of fractures. Non-isothermal ﬂ ow is considered within the fractures and the rock matrix, and the two ﬂ ow domains are coupled through a mass transfer term. The numerical model has been developed using standard ﬁ nite elements, with spatial discretisation achieved using the Galerkin method, and temporal discretisation using ﬁ nite di ﬀ erences. A full-scale ﬁ eld case geometric model, based on the Goldeneye depleted hydrocarbon reservoir in the North Sea, is developed and used for simulations. The in situ faults are modelled discretely as discontinuous surfaces in a three- dimensional matrix, including basement, reservoir, caprock and overburden layers. The faults are assumed initially to be low-permeable faults, with the same permeability as the caprock. However, the simulations show that their apertures (and as a result, their permeabilities) vary due to the thermoporoelastic e ﬀ ects caused by the injection of the relatively cold CO 2 . The change in the fracture apertures is mainly due to thermal e ﬀ ects; the reservoir layer undergoes contractions due to the cooling, signi ﬁ cantly increasing fault aperture in the region of the fault within the reservoir, whereas the fault aperture is reduced in regions within the caprock. Propagation of fractures under thermoporoelastic loading is investigated. Results show that the distance to the injection well, as well as spatial orientation of fractures with respect to the injection well, a ﬀ ect aperture evolution and potential growth of fractures. A sensitivity analysis is performed on the parameters a ﬀ ecting the fracture growth: minimum normal stress acting on the fracture plane, dip angle of the fracture, and the contact friction coe ﬃ - cient. It is found that low friction, low normal contact stress, or high in situ shear stress on the fracture surfaces may trigger propagation under combined mode II and III within the reservoir layer, or at the interface of the reservoir and caprock.


Introduction
Carbon Capture and Storage (CCS) is a key technology for addressing climate change and maintaining the security of energy supplies, while potentially offering important economic benefits (IPCC, 2005). In particular, depleted hydrocarbon reservoirs have the potential capacity to store significant quantities of carbon dioxide, making use of natural subsurface storage space which had previously successfully stored hydrocarbons for millions of years.
CO 2 is effectively sequestered by a combination of four principal mechanisms: dissolution, mineral trapping, capillary trapping, and structural trapping. CO 2 is stored at a supercritical state, having liquidlike density and gas-like viscosity, and can be found in three forms within the reservoir: CO 2 -rich phase, dissolved in pore water, and immobile by geochemically interacting with in situ minerals (Hangx et al., 2013). Thus, numerical models for geological storage sites not only require a detailed reservoir-caprock geological representation, but also require accurate modelling of the multiphase nature of CO 2 -brine displacement and trapping processes in order to determine the long-term fate of the injected CO 2 (Meer Van der, 1995). Injected CO 2 is modelled at the micro-scale incorporating relative permeability curves to describe the displacement mechanism and account for pore capturing (Yang et al., 2013). In contrast, when examining the caprock trapping mechanism, macroscopic, thermoporoelastic effects are often the focus, allowing the CO 2 modelling to be more simplistic (Jing, 2003). If simplified to a single-phase regime, and assuming equilibrium between rock and fluids, such an approach allows for the energy equation to be solved with reference only to temperature and pressure.
Geomechanics plays a crucial role in the performance assessment of the sequestration seal (Pan et al., 2013). During injection, CO 2 will reach the formation at lower temperatures than the reservoir, and so numerical models of injection should capture the effects of temperature contrast in addition to the fluid pressure increase (Vilarrasa et al., 2014). The importance of geomechanics during the injection of the cold CO 2 is two-fold: (i) During the injection time and around the injection point, the injection pressure may exceed the minimum in situ stress, which then may lead to the creation of hydraulic fractures. Nucleation of new fractures is a result of mechanical failure of the continuum reservoir/caprock (Rutqvist et al., 2008;Mathias et al., 2009;Vilarrasa et al., 2010;Pan et al., 2016). This failure can be avoided by controlling the injection pressure. For instance, U.S. regulatory standards restrict the maximum injection pressure to be less than the measured fractureclosure pressure of the caprock (USEPA, 1994). In depleted hydrocarbon reservoirs, such as in the Goldeneye field in the UK North Sea, the target pressure after injection is close to the initial reservoir pressure prior to hydrocarbon withdrawal, so it is conservative regarding the tensile failure of the rock. Therefore, a depleted reservoir is a safer option for CO 2 storage, as compared to an aquifer (UKCCS Geomechanics Summary Report, 2011). (ii) After injection ceases, the cold front continues to advance into the reservoir, inducing thermal contraction of the reservoir matrix, and unloading the normal compression on the existing fractures and faults. Thus, the aperture of the faults and fractures increases . This process can transform otherwise impermeable in situ faults and fractures into potential leakage pathways. Furthermore, such a reduction of the normal compression can reactivate the fault due to the shear, and may lead to induced seismicity (Rutqvist and Tsang, 2002;Lucier et al., 2006;Rutqvist, 2012).
The security of the CO 2 storage in the target site relies mainly on the seal capacity of the caprock, considering the buoyancy of CO 2 . The International Panel on Climate Change recommends that 99% or more of the injected CO 2 should remain in the targeted reservoir for at least a thousand years (IPCC, 2005). Caprock continuity depends on the depositional environment, and is a concern of geologic characterisation and site selection (Chapuis et al., 2011), whereas caprock integrity depends on the material properties and heterogeneities, as well as stress perturbations due to pressure and temperature changes. Caprocks are initially of very low permeability by definition, as they have sealed the hydrocarbons in place for many years (Bifani, 1986;Ketter, 1991;Ritchie and Pratsides, 1993). Thus, CO 2 seepage through the caprock can be considered negligible for depleted reservoirs, unless the caprock has become damaged during hydrocarbon withdrawal, or by leakage through existing wells (Nordbotten et al., 2008). The risk for permeability increase as a result of a volumetric strain during production is expected to be negligible (UKCCS Geomechanics Summary Report, 2011). Depleted reservoirs can store buoyant fluids, and also if they are un-faulted, the overlying strata are capable of trapping significant gas accumulations (Williams et al., 2013). However, caprocks are susceptible to fracturing or fracture/fault opening, which can act as an upward migration path for stored carbon dioxide. The success of CO 2 containment relies on maintaining caprock integrity by ensuring that fluctuations in pressure or temperature, or any other state-changing events, do not create new leakage pathways, or enhance pre-existing pathways. To analyse these processes, coupled thermoporoelasticity (THM) equations need to be solved. This can most conveniently be done using a robust numerical model such as the finite element method (FEM), which can easily account for material inhomogeneities, as opposed to alternative methods such as the boundary integral or displacement discontinuity methods.
CO 2 leakage through the caprock undermines the caprock seal, and may lead to the mobilisation of stored gas. The leakage flow can be a result of fracture network connectivity in the caprock prior to, or after, injection (Smith et al., 2011), or the result of hydraulic fracture propagation induced during CO 2 injection (Pan et al., 2013). Existing permeable fractures and faults in the caprock can greatly affect the reservoir CO 2 containment capacity (Elkhoury et al., 2015) and facilitate the build-up pressure transmission to upper layers (Mbia et al., 2014) which may cause surface uplift above and around the injection wells during CO 2 injection (Rutqvist et al., 2010). Nucleation of new fractures, propagation of existing fractures, and reactivation of existing sealed faults due to stress perturbation from pore pressure changes in the underlying reservoir can enhance the caprock permeability (Candela et al., 2011;Renard et al., 2012;Mason et al., 2013;Gheibi et al., 2017). Propagation of existing fractures and faults through the reservoir and caprock is more critical, as it may jeopardise the caprock integrity and seal. Therefore, monitoring seismicity in order to track the opening of fractures has been recommended during CO 2 injection (Vilarrasa et al., 2013). Seismic activities have been observed during fluid injection into underground reservoirs (Healy et al., 1968;Miller et al., 2004;Cappa and Rutqvist, 2011;Wang et al., 2016;Gheibi et al., 2018) and during hydrocarbon production (Grasso and Feignier, 1990). Therefore, it is also important to assess the propagation potential of existing faults and fractures during reservoir production and CO 2 injection using a three-dimensional model.
A three-dimensional coupled thermoporoelastic (THM) model is utilised in this research to study the reservoir compaction and fracture deformation during the injection of cold carbon dioxide. The fully coupled FEM model has been previously used for modelling hydraulic fracturing Salimzadeh et al., 2017a,b) and for heat extraction from fractured geothermal reservoirs (Salimzadeh et al., 2018a,b). In this model, fractures are represented discretely using discontinuous surfaces in a 3D mesh, and the contact stress between fracture surfaces is resolved using the gap-based Augmented Lagrangian (AL) method for accurately enforcing the contact constraint in fractured media (Nejati et al., 2016). The contact model has been extended to account for thermal and hydraulic loadings. In this research, the mechanical deformation-contact model is iteratively coupled with nonisothermal flow model (Salimzadeh et al., 2018b). Iterative coupling between flow and deformation is used here, rather than full coupling, due to its lower computational cost, and greater stability when simulating processes in low-permeability caprock. The large contrast between the permeability of the reservoir and caprock layers creates solver convergence problems in the fully coupled system, which are manageable for smaller systems, but become an impediment for larger and longer simulations. In this research, the onset and direction of the fracture propagation under THM loading are evaluated using the three stress intensity factors (mode I, II, III), computed using the displacement correlation (DC) method. A full-scale model of the Goldeneye reservoir, located offshore in the UK North Sea, is developed as a study case for CO 2 injection simulations. Injection of CO 2 is performed over ten years, and the simulation of storage is continued for 160 years to study the thermal effects that usually arise at longer times, compared to the poroelastic effects.

Model setup -Goldeneye reservoir
The Goldeneye depleted gas condensate field is located offshore in the UK North Sea. The reservoir rock consists of the weakly consolidated Captain Sandstone, overlain by the low-permeability Rødby unit as caprock. In this study, a full-scale model of the Goldeneye field is built using the data available in the literature (UKCCS Geomechanics Summary Report, 2011;Hangx et al., 2013), as shown in Fig. 1. Several formations are grouped into larger layers to reduce the simulation difficulties and cost. Four layers are included in the model: 'overburden' or top layer (Nordland group; Coals; Dornoch; Ekofisk), caprock (Rødby), reservoir (Captain E, D, C, and A; Valhall and Scapa), and 'underburden' or basement layer (Humber; Heron). The total size of the model is 10 × 6 × 4 km. The properties of the layers are listed in Table 1. Vertical stress is calculated from weight of the overburden, whereas the minimum principal stress is estimated from the leak-off tests (LOT), which show a normal stress regime. Vertical stress is the maximum principal stress, and both horizontal stresses are assumed to be equal to the minimum principal stress. The Captain sandstone is located at a depth of 2600 m; hence, the vertical stress is set to 55 MPa, and the horizontal stresses are set to 45 MPa.
Production from the Goldeneye field started in 2004 and ceased in 2010, during which time the reservoir pressure decreased from an initial value of 25 MPa to around 15 MPa. Injection of CO 2 is planned for a period of ten years, with a maximum of injection pressure of 28 MPa. The reservoir and its surrounding medium underwent cycles of loading (during hydrocarbon production) and unloading (during CO 2 injection). Therefore, the properties of the rock are expected to vary during hydrocarbon production and CO 2 injection phases (UKCCS Geomechanics Summary Report, 2011). To account for this "hysteresis" effect, the Young's modulus of the reservoir and surrounding layers is increased to 20 GPa as given in Table 1. The in situ temperature at the reservoir depth is expected to be 83°C, and the CO 2 is planned to be injected at 20°C (UKCCS Geomechanics Summary Report, 2011). The thermal conductivity, heat capacity and the rock density are set to 1.5 W/mK, 900 J/kg°C, and 2100 kg/m 3 , respectively (see Table 1). A total of seventeen major faults have been mapped in the field, penetrating the caprock and the reservoir. The faults are near vertical, with a dominant NW-SE trend. Four of the faults, named Fault04, Fault05, Fault06 and Fault11, are included in the full-scale model to investigate their deformation, aperture variation, and propagation during CO 2 injection. These faults are positioned at a variable distance from the injection well, and have different orientations with respect to the injection well, which is the source of the stress perturbation. The case assumes a worst-case scenario, in which a set of faults has already formed at the interface connecting reservoir and caprock. Sets of discontinuities tend to enhance stresses, and tend to be more prone to growth, as compared to stand-alone fractures and faults (Thomas et al., 2017). The size and vertical positioning of the faults are uncertain, so in the present simulations two configurations are considered: (i) the faults extend from the top layer to the basement layer crossing through caprock and reservoir layers; (ii) the fractures extend half-way through the caprock and halfway through the reservoir. The former configuration is used to study the aperture variation and the potential of faults to act as leakage pathways, whereas the latter configuration is used to study the chance for slip and further propagation. It is assumed that the thickness of the faults and fractures are negligible compared to their dimensions, so they can be represented as surfaces in the three-dimensional model. The orientation and the position of the faults and fractures in the horizontal plane for both scenarios are similar, as given in Fig. 1.

Computational model
In the present approach, fractures are modelled as surface discontinuities in the three-dimensional matrix. Fluid flow through fractures is treated separately from flow through the matrix, thus, five interacting conceptual models are considered: (i) mechanical deformation, (ii) flow through the matrix, (iii) flow through the fracture, (iv) heat transfer through the matrix, and (v) heat transfer through the fracture. The deformation model is expressed satisfying the condition of equilibrium on a representative elementary volume (REV) of the porous medium. Fracture surfaces are not traction-free in the present model, and hydraulic loading, as well as the tractions due to the contact between fracture surfaces, are applied on the fracture walls. Assuming negligible shear tractions exerted from the fluid on the fracture surfaces, the fluid pressure is applied only in the normal direction of the fracture wall. The tractions on the fracture surfaces are where σ c is the contact traction on the fracture surfaces, p f is the fluid pressure in the fracture, and n c is the outward unit normal to the fracture surface (on both sides of the fracture). When two surfaces of a fracture are in partial contact at the micro-scale, the mean aperture of the fracture is a function of the normal contact stress. In this study, the classic Barton-Bandis model (Bandis et al., 1983;Barton et al., 1986) is used to calculate the fracture aperture under contact stress as where a f is the fracture aperture, σ n is the normal component of contact stress, a 0 is the fracture aperture at zero contact stress, and, a and b are model parameters. The normal contact stress is directly given by the contact tractions in the contact mechanical model. Fractures are modelled as surface discontinuities within the matrix; therefore, the contact problem arises and the contact stresses (normal and shear) are required to be computed in order to avoid the inter-penetration of fracture surfaces under compressive loading. The Augmented Lagrangian (AL) method has been utilised for enforcing the contact constraint, by combining the Lagrange multiplier and penalty methods to exploit the merits of both approaches (Wriggers and Zavarise, 1993;Puso and Laursen, 2004). A sophisticated algorithm is used for the treatment of frictional contact between the fracture surfaces, based on isoparametric integration-point-to-integration-point discretisation of the contact contribution. Contact constraints are enforced by using a gap-based AL method (Nejati et al., 2016). In this model, penalties are defined at each timestep as a function of local aperture, so that they are larger away from the fracture tips, and reduce to zero at the tips. For the frictional contact, the classical Coulomb friction condition is considered, in which the maximum shear traction equals the normal traction multiplied by the friction coefficient. Thus, the contact mode can evolve from the "stick" to "slip" mode if the shear component of the contact traction exceeds the maximum shear traction (shear resistance). In this work, the contact model in Nejati et al. (2016) is further extended to account for hydraulic and thermal loadings.
Flow through the fracture is defined using the cubic law while flow through the matrix is based on Darcy's law (Zimmerman and Bodvarsson, 1996). The mass transfer between the fracture and the rock matrix is defined as a leakoff term Khalili, 2015, 2016 where k n is the intrinsic permeability of the rock matrix in the direction normal to the fracture (in the direction of n c ). The governing equation for heat transfer through the rock matrix can be obtained by combining Fourier's law with an energy balance for saturated rock. It is assumed that the fluid velocity in the rock matrix is low enough such that the solid grains and the fluid in the rock matrix are always in local thermal equilibrium. Heat is also exchanged between matrix and fracture fluid by conduction through the fracture surfaces, and by advection through the leakoff mass exchange term, as (Salimzadeh et al., 2018a) where λ n is the average thermal conductivity of the rock matrix along the direction normal to the fracture (in the direction of n c ). The final flow and heat transfer models are derived combining the conservation equation with linear momentum balance equation (Salimzadeh et al., 2018b). More details on the fully coupled and iteratively-coupled thermoporoelastic (THM) models for deformable matrix with discrete fractures can be found in Salimzadeh et al. (2018a,b), respectively. In the iteratively-coupled model, used in the present study, TH and M are coupled sequentially to improve the performance and the solver convergence for low-permeability rocks.
Concentrations of stress around the fracture tips due to pressure and temperature perturbations induce fracture propagation. Stress concentrations can be quantified locally at each tip using stress intensity factors (SIFs). The SIFs are key parameters in evaluating and predicting the onset of fracture growth, and growth direction. In this study, three stress intensity factors (SIFs) for three modes of fracture opening are computed using the displacement correlation (DC) method. The DC method is computationally cheap and is able to yield very good approximations of the SIFs (Kuna, 2013). The three SIFs are mode I (K I ) for opening due to tensile loading, mode II (K II ) for in-plane shearing due to sliding, and mode III (K III ) for out-of-plane shearing due to tearing. The crack grows if and when the equivalent stress intensity factor K eq overcomes the material toughness (K ic ). The equivalent SIF in the direction of propagation (θ p ) is calculated as (Schöllmann et al., 2002) ⎜ , and θ p is the propagation angle.
The in-plane propagation angle (θ p ) and out-of-plane deflection angle (ψ p ) are determined using a modified maximum circumferential stress method that takes into account modal stress intensity factors under mixed loading (Schöllmann et al., 2002). The in-plane propagation angle (θ p ) is assumed to be perpendicular to the maximum circumferential stress σ 1 , thus, θ p can be calculated by The out-of-plane deflection angle ψ p is defined by the orientation of σ 1 as At every timestep, the SIFs and growth computations are performed at fifty locations along the fracture tip.
The governing equations are solved numerically using the finite element method. Spatial and temporal discretisations are realised using the Galerkin method and finite difference techniques, respectively. In the mechanical deformation-contact model, the displacement vector u, and in the TH model the fluid pressures p m and p f , and matrix and fracture fluid temperatures T m and T f , are taken as primary variables. The discretised equations are implemented in the Complex Systems Modelling Platform (CSMP++, also known as CSP), an object-oriented application programme interface (API), for the simulation of complex geological processes and their interactions (formerly CSP, cf. Matthäi et al., 2001;Paluszny and Zimmerman, 2011). Quadratic unstructured elements are used for spatial discretisation of surfaces (quadratic triangles) and volumes (quadratic tetrahedra). The triangles on two opposite surfaces of a fracture are matched with each other, but do not share nodes, and duplicate nodes are defined for two sides of a fracture. The triangles are matched with faces of the tetrahedra connected to the fractures, and they share the same nodes. Fracture flow and heat equations are solved only on one side of the fracture, whereas, matrix deformation, fluid flow and heat transfer equations are accumulated over the volume elements. The ensuing set of linear algebraic equations is solved at each timestep using the algebraic multigrid method, SAMG (Stüben, 2001). The proposed THM model has been validated extensively against several analytical, experimental and numerical results

Simulation results and discussions
The full-scale model described in Section 2 has been used to investigate the effect of cold CO 2 injection on fracture aperture variation and fracture propagation. Constant horizontal (45 MPa) and vertical (55 MPa) far field stresses are applied, the reservoir pressure prior to the CO 2 injection is set to 15 MPa, and the initial temperature is set to 83°C. The injection is simulated through prescribed injection pressure that linearly increases to 28 MPa over ten years, and is then kept constant and equal to 28 MPa. The injection is performed through a single well that is located in the centre of the model, whereas in the real situation, injection may be performed through several wells. After ten years, the pressure at the well is kept at 28 MPa to avoid reservoir pressure drop, and to increase pressure in those parts of the reservoir far from the well. Throughout the simulation, the injection well is kept at a constant temperature that differs (20°C, 40°C, 83°C) for each simulated case. The top of the model is open for flow, while the vertical displacements are constrained at the bottom of the model.

Fault deformation during cold CO 2 injection
In this scenario, faults extend from the basement to the top layer of the geometric model, crossing through both reservoir and caprock layers. The faults are considered to be in "stick" mode by assigning a large friction coefficient. Sliding of the faults is not analysed in this scenario. Faults are initially assigned an equivalent permeability similar to the caprock so that they are not initial leakage pathways. The contact stress-fracture aperture relationship is defined by the Barton-Bandis model with model parameters of a = 1.6 × 10 −10 m/Pa, b = 1.333 × 10 −7 1/Pa, and a 0 = 9.65 × 10 −4 m. The fracture contact stress prior to the injection is 30 MPa, thus the initial aperture is a f = 4.81 × 10 −6 m which gives an equivalent permeability-thickness of = × − a /12 9.26 10 f 3 18 m 3 . Four cases are simulated: (i) injection temperature T inj = 20°C and Young's modulus E = 20 GPa, (ii) injection temperature T inj = 40°C and Young's modulus E = 20 GPa, (iii) injection temperature T inj = 20°C and Young's modulus E = 10 GPa, and (iv) isothermal case with Young's modulus E = 20 GPa. The Young's modulus and Poisson's ratio for all layers are assumed to be the same to avoid stress concentration in a particular layer. The model is discretised using 76,679 quadratic tetrahedral and triangular elements and 100,229 nodes. The simulations are performed for 160 years. The timestep starts initially from ten days and increases by a factor of 1.1 till it reaches the prescribed maximum timestep of ten years. A total of seventy-nine steps are run to reach the 160 years of simulation after the injection starts. The results for the fluid pressure in the faults, the fluid temperature in the faults, the faults' aperture, and the temperature plume in the reservoir at the end of the simulation are shown in Fig. 2-5 for different cases.
In the first case, while the pressure distribution in the faults is relatively similar to each other, the temperature distributions are quite different. Fault11 and Fault05 that face the well and are closer to the well, experience higher temperature change compared to Fault04 which is located further away from the well, and Fault06 that has one side towards the well. The top layer acts as drainage pathway for the CO 2 flow, thus the pressure does not build-up in this layer, instead it remains close to the prescribed pressure at the top boundary, 15 MPa. The aperture distribution in the faults closely follows the temperature distribution pattern, indicating that fault deformation is mainly induced by thermoelastic stresses. Fault aperture increases locally up to 0.5 mm from the initial value of less than 0.005 mm, by almost 100 times. However, it is important to notice, that this increase in the aperture is restricted to the regions located within the reservoir layer, and the temperature plume does not propagate significantly into the caprock layer during the 160 years of simulation. The maximum aperture is observed in the Fault05 and Fault11 that are submerged into the cold plume as shown in Fig. 2. The profile of the fluid pressure along a vertical line passing through each fault and through the intact matrix show that the faults are facilitating pressure diffusion through the caprock, and the CO 2 pressures are higher in the faults within the caprock. This is reflected by the aperture of the faults in the caprock, which has been increased to about 0.04-0.05 mm, roughly ten times the initial value. According to the Barton-Bandis fracture stiffness model used in the simulations, an increase of 5 MPa in the fluid pressure inside the faults will result in an aperture increase to 0.05 mm. Therefore, the elevated fluid pressure can be considered responsible for the aperture increase in the faults. Interestingly, Fault04, less affected by the thermoelastic deformations within the reservoir, shows higher aperture and higher fluid pressure in the caprock. This is due to the additional compression of the caprock layer from the cooling of the reservoir layer. Cooling down the reservoir layer results in the contraction of the matrix and reduction of the matrix volume. The reduction of the matrix volume reduces the layer's share of the applied far-field stress and thus the share of other layers such as the caprock layer increases. The compaction of the caprock due to the injection of cold CO 2 into the reservoir layer is reported in other studies (e.g. Vilarrasa and Laloui, 2015). In this regard, the pressurising of the reservoir matrix and the cooling of the reservoir matrix have opposite effects. Therefore, the orientation and the distance of the faults with regard to the temperature plume control the extent of faults deformation and aperture variation.
When the temperature of the injected CO 2 increases, in the second case, the thermoelastic loading reduces. As a result, the magnitude of the aperture increase in the faults within the reservoir layer reduces, as shown in Fig. 3. The highest aperture reduces to less than 0.3 mm, again focusing in Fault05 and Fault11, both faults with less favourable orientation and closest to the injection well. The aperture field still closely follows the temperature distribution, reflecting the fact that the thermoelastic deformation continues to play the main role in increasing fault aperture. Unlike the regions of the faults located within the reservoir, the regions within the caprock show little sensitivity to the temperature change in the injected fluid. This is mainly because the deformation in the faults within the caprock is due to the increased fluid pressure. Fault04 is an exception, as its aperture in the caprock layer reduces further due to lower thermoelastic deformation in this case.
In the third case, all layers are assumed to have a lower elasticity coefficient of E = 10 GPa. As compared to the previous cases, the deformable matrix reduces its thermoelastic deformation contribution. The maximum aperture reduces to 0.2 mm (again) on Fault05 and Fault11. However, the aperture field in Fault04 is less affected by the lower Young's modulus, as shown in Fig. 4.
In the fourth and last case, CO 2 is injected under isothermal conditions. Thus, the response of the reservoir and caprock is purely due to poroelastic stresses. In this case, the aperture field on all faults reduces to less than 0.1 mm. In this case, all faults approximately exhibit similar responses to the pressure perturbation (Fig. 5). This case clearly demonstrates that the HM model underestimates the deformations induced during CO 2 injection, considering that the temperature of the injected fluid is not the same as the reservoir temperature. It also shows that the isothermal HM model consistently overestimates fault aperture in the caprock and overburden layers, as compared to non-isothermal models. The profile of the aperture distribution along a vertical line passing through the faults is compared for all cases in Fig. 6. Fault05 and Fault11 are most affected by the thermoelastic deformation in the reservoir, while Fault04 is least affected. The thermoelastic deformation increases the faults' aperture within the reservoir layer, but it reduces the faults' aperture within the caprock. Thus, the HM model shows higher fault aperture within the caprock. This is due to the additional compression that the caprock layer receives from the thermal contraction in the reservoir layer. Therefore, the higher the aperture evolution in the reservoir layer, the lower the aperture evolution in the caprock. Growth is not observed in any of the faults, at any stage, primarily due to the high friction coefficient assigned to the faults.

Fracture propagation during cold CO 2 injection
In this section, smaller fractures are considered in the location as the faults shown in Fig. 1. The elliptic fractures are named: Fracture04, Fracture05, Fracture06, and Fracture11 which correspond to Fault04, Fault05, Fault06, and Fault11, respectively. The fractures, each of radius 250 m, are considered vertical (similar to the faults) and extend halfway into the reservoir layer and halfway into the caprock. The vertical direction of the fractures reduces the in situ shear stresses on the fractures. As compared to the case in which faults are modelled, for fractures, the contact model is run to evaluate "slip" mode and the friction coefficient (μ s ) is varied between 0.1 and 0.8 for different cases (as opposed to the high friction coefficient assumed for the faults case). The fluid-filled fractures and faults can have very small friction coefficients. The contact tractions are resolved using two Uzawa iterations, and three to five Newton-Raphson iterations per Uzawa iteration. The simulations are run for fifty years, and the stress intensity factors are computed at fifty locations on the fracture tip, and used to evaluate the onset and the direction of propagation. The fracture toughness is set to K ic = 1 MPa m 0.5 . This value is reasonable for sandstones (Senseny and Pfeifle, 1984). A sensitivity analysis is performed on the minimum in situ stress, dip angle of fractures and the friction coefficient. The minimum in situ stress (along the y-axis) in the Goldeneye reservoir is expected to be 45 MPa; however, it is reduced to 40 and 35 MPa in these analyses, whereas the vertical and maximum horizontal in situ stresses are kept constant at 55 and 45 MPa, respectively. Lowering the minimum in situ stress increases the in situ shear stress acting on the fracture planes. The fractures' dip angle is also varied between vertical and θ = 80°. Tilting the fractures from the vertical plane induces initial shear stress on the fracture planes. The contact friction coefficient is varied from 0.1 to 0.8 to investigate whether fractures are stable during simulation time of fifty years, or propagation occurs. Different cases are summarised in Table 2, varying between a best case scenario (σ y = 45 MPa, θ = 90°) and a worst case scenario (σ y = 35 MPa, θ = 80°). In the first case, in order to allow the fracture propagation, the friction coefficient is reduced to μ s = 0.1. As expected, the propagation occurs in Fracture11. Fracture05 and Fracture11, similar to Fault05 and Fault11, undergo the most thermoelastic deformation; however, the affected area in Fracture11 is closer to the tip of the fracture, whereas in Fracture05, it is close to the centre.
As shown in Fig. 7a, the location of the fracture propagation is at the interface of the reservoir and caprock layers, and the propagation occurs mainly under mode III i.e. out-of-plane shear (tearing). Propagation is observed in a very restricted region in Fracture11, in the region closest to Fracture06. The equivalent stress intensity factor (K eq ) exceeds the fracture toughness only at three locations (out of fifty) as given in Table 2. The small values of K I reflect the numerical accuracy of the contact model in resolving the inter-penetration. In the second case, the fractures are tilted 10°from the vertical plane (θ = 80°), and the simulations are repeated for μ s = 0.2. Tilting the fractures from the vertical plane increases the in situ shear stresses on the fracture surfaces; as a result, the propagation occurs within fifty years of simulation. The propagation is under combined mode II and III, and the contribution of mode II has increased as shown in Fig. 7b. The location of propagating tips has moved downward into the reservoir layer. The SIFs at the propagating tips are given in Table 2. The fractures in both cases (case I: θ = 90, μ s = 0.1, case II: θ = 80, μ s = 0.2) are stable without the thermal contribution. As shown in Fig. 7c, in case I, the in situ shear on the fracture is almost zero, and in case II, the in situ shear (1.71 MPa) is less than the μ s σ = 4.06 MPa. The contribution of the thermal stress is two-fold: (i) it reduces the normal traction on the fracture surfaces by cooling down the matrix, (ii) it increases the shear traction on the fracture by pulling the reservoir matrix towards the cold plume. The contribution of the thermal stresses causes the fracture to slide in both cases. The components of the contact traction (x, y, z components) along a vertical line passing through the centre of Fracture11 for two cases of the vertical fracture with μ s = 0.1, 0.2 at the exact timestep when the case with lower friction fails are plotted in Fig. 8. Results show that the shear component of the contact traction (z component) in case I: θ = 90, μ s = 0.1 has reached the maximum friction according to the classic Coulomb law. This occurs at the interface between reservoir and caprock layers and results in the stress redistribution over the fracture.
Four other cases are simulated in which the minimum horizontal stress has reduced to 40 and 35 MPa, while the friction coefficient is gradually increased to 0.8 for the worst case. The cases in which propagation of Fracture11 occurs within fifty years of simulation are shown in Table 2. As the in situ conditions (minimum horizontal stress and dip angle of the fracture) worsen, the in situ shear stress acting on the fracture plane increases, the contribution of mode II increases. The propagating tips move downward into the reservoir layer and propagation direction is similar to the second case (Fig. 7b). It is worth mentioning that no propagation in the above-mentioned cases occurs when the injection temperature increases to 40°C.

Conclusions
The injection of cold CO 2 into the Goldeneye field in the UK North Sea is simulated using a coupled thermoporoelastic model. The fullscale model contains four layers: top, caprock, reservoir, and basement layers, as well as four faults that have different orientation and distance from the injection well. Two sets of simulation are performed: in the first set, the large in situ faults crossing reservoir and caprock layers are modelled.
Simulations demonstrate that aperture distributions along faults that lie across rock layer boundaries are strongly local, with higher aperture change in parts within the reservoir layer. The magnitude of the fault aperture change during CO 2 sequestration depends on rock   S. Salimzadeh et al. International Journal of Greenhouse Gas Control 74 (2018) 130-141 properties, injection pressure and temperature, the location and orientation of the fault with respect to the injection well as the source of cold plume. Thus, only a three-dimensional model is capable of capturing the effect of geometry of fractures. The thermoelastic stresses dominate the reservoir deformation and strongly affect the distribution of fault apertures within the reservoir, while the aperture in the caprock is controlled primarily by poroelastic deformation. Matrix contraction due to the injection of cold fluid in the reservoir layer induces compression on the caprock layer that reduces the aperture of the faults within the caprock. Injecting warmer fluid, or having a more compliant  matrix (modelled by a lower Young's modulus) reduces thermoelastic effects.
In the second set, shearing of the existing fractures due to thermoporoelastic deformations is also investigated by modelling smaller fractures. A sensitivity analysis on the minimum in situ stress, fracture dip angle, and friction coefficient is performed. Results show that the propagation occurs in the most critically located and oriented fracture, and the mode of fracture propagation varies between mode II and III. Mode III occurs in most stable case with very low friction coefficient (μ s = 0.1), while mode II occurs in cases with higher in situ shear stress acting on the fracture prior to CO 2 injection. Higher in situ shear stress can be a result of lower minimum in situ stress and/or dip angle that deviates from vertical (θ = 80). For the worst case (σ y = 35 MPa, θ = 80), a friction coefficient as high as 0.8 cannot prevent fracture failure due to thermoporoelastic deformation of the reservoir matrix within fifty years of simulation.