
==== Front
iScience
iScience
iScience
2589-0042
Elsevier

S2589-0042(24)01778-4
10.1016/j.isci.2024.110553
110553
Article
Effects of mechanical properties of the underburden on induced seismicity along a basement fault during hydrogen storage in a depleted reservoir
Burtonshaw James E.J. james.burtonshaw16@imperial.ac.uk
12∗
Paluszny Adriana apaluszn@imperial.ac.uk
1∗∗
Mohammadpour Aslan a.mohammadpourshoorbakhlou@imperial.ac.uk
1∗∗∗
Zimmerman Robert W. r.w.zimmerman@imperial.ac.uk
1∗∗∗∗
1 Department of Earth Science & Engineering, Imperial College London, South Kensington, London SW7 2AZ, UK
∗ Corresponding author james.burtonshaw16@imperial.ac.uk
∗∗ Corresponding author apaluszn@imperial.ac.uk
∗∗∗ Corresponding author a.mohammadpourshoorbakhlou@imperial.ac.uk
∗∗∗∗ Corresponding author r.w.zimmerman@imperial.ac.uk
2 Lead contact

24 7 2024
16 8 2024
24 7 2024
27 8 11055315 2 2024
21 5 2024
17 7 2024
© 2024 The Author(s)
2024
https://creativecommons.org/licenses/by-nc-nd/4.0/ This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).
Summary

This study models the geomechanical deformation of a depleted gas field, wherein gaseous hydrogen is stored in a North Sea reservoir, and is cyclically injected and withdrawn. A fault is modeled within the underburden, and its slip is investigated during a three year storage period. Parametric simulations are conducted to study the influence of the underburden mechanical properties, such as Young’s modulus, Poisson’s ratio, and permeability on induced seismicity. The fault is predominantly in stick during the bulk of the injection, storage, and withdrawal periods, but minor fault slip (<4 mm) occurs shortly after a change in operational regime. The Young’s modulus of the underburden unit has the strongest control on fault slip. To reduce the seismic hazard, an underburden with low Young’s modulus (<15 GPa), high Poisson’s ratio (>0.25), low Biot coefficient, and low permeability (<1×10−19 m2) is found to be most suitable for hydrogen storage.

Graphical abstract

Highlights

• H2 was injected, stored, and withdrawn from a faulted depleted gas field for 3 years

• Young’s mod (E) had key control on slip (significantly decreased as E increased)

• Max potential event had magnitude 3.09 when the Poisson ratio (ν) was smallest

• Low E (<15 GPa), high ν (>0.25), low Biot coeff./permeability minimized magnitude

Structural geology; Petrophysics; Energy engineering

Subject areas

Structural geology
Petrophysics
Energy engineering
Published: July 24, 2024
==== Body
pmcIntroduction

The long-term storage of gaseous hydrogen in depleted oil and gas fields could be an essential technology in meeting future energy needs in a renewable energy world. Hydrogen has emerged as a key player in the global energy transition due to its versatility,1 potential for decarbonization,2,3 storage potential,4,5 potential for long duration grid balancing,6,7 and its ability to address challenges in sectors that are difficult to electrify directly, such as heavy industry8,9,10 and long-haul transportation.11 In a low-carbon energy world, an annual energy dichotomy will be observed, in which for half of the year an excess of renewable energy is produced, and for the other half, a deficit of energy exists.12 In order to avoid the use of fossil fuels as seasonal fuels to satisfy this energy deficit, green hydrogen can be synthesized from excess renewable energy during the periods of energy abundance,4,12 and then injected into subsurface porous media such as depleted gas fields where it is stored. Subsequently, the hydrogen can be withdrawn during the periods of energy deficit, avoiding the need for fossil fuels. The hydrogen is then either combusted13 as a fuel or used as a reactant in a fuel cell14 where the hydrogen molecules are stripped of their electrons to generate an electric current that can be used to generate power.

Storage in depleted gas reservoirs is the only pragmatic storage technique capable of handling sufficient quantities of hydrogen to satisfy future demand. One typical depleted gas field with porosity, ϕ = 0.15, thickness, h = 150 m, area, A = 4,000,000 m2, hydrogen saturation, Sh = 0.75, hydrogen recovery factor, Rf = 0.8, and hydrogen density, ρh = 9 kg/m3, can store the same volume of hydrogen as nearly 120 typical salt caverns.4

Hydrogen can be produced via a number of physical processes and chemical reactions making it a highly versatile fuel, and is often given a color prefix to denote which chemical or physical mechanism has been implemented in its production (see Table 1). Hydrogen can be produced via clean pathways (green, blue, white, or pink) or via more traditional, greenhouse gas emissions-heavy routes (brown, black, or gray). Clean hydrogen is the main-focus of future geological hydrogen storage in the context of the energy transition. While green hydrogen is the ultimate goal for geological hydrogen storage, renewable energy needs to become more widespread to make it a more cost-effective method of hydrogen production. In the meantime, blue and turquoise hydrogen offer excellent storage alternatives. In particular, turquoise hydrogen can sustainably use the Earth’s substantial natural gas resources without releasing carbon dioxide, since the carbon atoms of the natural gas become carbon soot, offering an excellent route for early commercial-scale geological hydrogen storage.Table 1 Different types of hydrogen, the chemical or physical route by which they were produced, their cost, their percentage of the current global hydrogen mix and their cleanliness in terms of their CO2 emissions

Type	Reaction/Process	Reactants	Products	Cost ($/kg)	Current % of H2 supply	Cleanliness	
Green	Wind/Solar/Geothermal/Hydropower-derived Electrolysis	H2O	H2, O2	3.38–12.0015,16	0.517	Carbon-neutral	
Blue	SMR, WGS, CCS	CH4, H2O	H2, CO2	1.50-3.50)17,18,19	117	Carbon-neutral	
Brown	Lignite Coal Gasification	C,O2,H2O	H2, CO2	0.94–2.2519,20	117,21	Very heavy CO2 route	
Black	Bituminious Coal Gasification	C,O2,H2O	H2, CO2	0.94–2.2519,20	2117,21	Very Heavy CO2 route	
Gray	SMR, WGS, no CCS	CH4, H2O	H2, CO2	1.00–1.7522,23	7617,24	Heavy CO2 route	
Pink	Nuclear-derived Electrolysis	H2O	H2O2	2.24–5.9225	0	Carbon-neutral	
White	Serpentinization, Hydraulic Fracturing	Mg1.5 Fe0.5 SiO4, H2O	Fe3O4, Mg3Si2O5(OH)4	–	0	Carbon-neutral	
Purple	Nuclear-derived high-temperature electrolysis	H2O	H2, O2	2.24–3.7325	0	Carbon-neutral	
Red	Nuclear-derived High-temperature thermocatalytic/thermochemical H2O splitting	H2SO4, I2	O2, H2	2.18–5.6525	0	Carbon-neutral	
Yellow	Solar-derived Electrolysis	H2O	H2, O2	3.38–6.8426	unk.	Carbon-neutral	
Turquoise	Pyrolysis	CH4	C, H2	1.60–2.4727	< 0.5	Theoretically carbon-neutral	
Orange	Redox	CO2, MgO/CaO, FeO, H2	H2, Fe2O3, MgCO3/CaCO3	unk.	0	Carbon-neutral	

Green hydrogen is produced from the electrolysis of water, in which the source of electricity is derived from renewable energy such as solar, wind, or hydropower.22 In an electrolytic cell, water is split in the presence of electric current and an electrolyte such as sulfuric acid to produce hydrogen at the cathode and oxygen at the anode.28 Currently, three main types of electrolyzer are used: alkaline electrolyzers, polymer electrolyte membrane electrolyzers, and solid oxide electrolyzers. Blue hydrogen is produced from methane involving a two-stage scheme. The first stage is the steam methane reaction in which methane is converted to carbon monoxide and hydrogen in the presence of 700–1000°C steam and a nickel catalyst.29,30 The second stage is the water-gas shift reaction,30 in which the toxic carbon monoxide that poses a disposal problem, is reacted with further hot steam to yield CO2 and more hydrogen. Critically, the CO2 must then be captured and sequestered in depleted oil and gas fields, salt caverns, or lined rock caverns. Therefore, while blue hydrogen generates CO2 during its synthesis, by storing it, none is released to actively participate in the greenhouse effect. White hydrogen is naturally occurring hydrogen, usually produced from hydraulic fracturing (“fracking”) of subsurface geological deposits. It is unique within the hydrogen color spectrum because it is the only type of hydrogen which does not require production from a chemical or electrical process—it is natively ready for immediate use. The main mechanism for white hydrogen production is the serpentinization of ultramafic igneous rocks. Serpentinization refers to the redox reaction which occurs when Fe-rich minerals such as olivine are exposed to hot hydrothermal fluids, typically yielding serpentine and magnetite, and evolving hydrogen gas, which is then stored in natural fractures and any pore space which has been etched out from reactive subsurface fluids. However, there are currently limited plans to exploit white hydrogen deposits, and no substantial formal assessment of its resource potential has taken place.31 Interestingly, these white hydrogen deposits appear to be geographically diverse, increasing their future potential as a wide-scale source of hydrogen. Truche and Bazarkina31 provide an initial scoping of known white hydrogen deposits, including intra-cratonic seepages, peralkaline granites, ophiolite deposits, mid-ocean ridge deposits, and clay-rich adsorption deposits. Deposits appear to be found across the globe, although appear to be less concentrated in Asia and Eurasia (although this could merely be a result of less publications in these geographies and/or less interest in these resources). Artificial serpentinization can even be initiated in the laboratory to generate hydrogen using diamond-anvil cells and adding a source of aluminum to the serpentinization reaction, either in the form of ruby microspheres or aluminum hydroxide.32

Pink hydrogen is hydrogen produced through electrolysis in which the source of electricity is provided by nuclear power. In the case of SOE electrolyzers, hydrogen may also be termed pink hydrogen where the source of steam is produced by heating water with thermal energy derived from nuclear fission.33 Pink hydrogen is yet to see implementation at a commercial scale. Subsets of pink hydrogen include purple and red hydrogen. Purple hydrogen uses thermal energy from nuclear power to perform water splitting by high temperature electrolysis, while red hydrogen uses the nuclear-derived thermal energy to perform high-temperature thermochemical and thermocatalytic water splitting (i.e., via the sulphur-iodine cycle34). Turquoise hydrogen is a promising form, but its clean nature is still debatable. This hydrogen refers to the pyrolysis of natural gas,35 yielding hydrogen gas and solid carbon (“carbon black”) using a thermal cracking reaction. Natural gas is pumped into a quartz bubble column reactor, partially filled with a molten metal, typically gallium, magnesium, copper, or tin,36 through an inlet. The molten metal acts as a thermal transfer medium, heating the rising bubbles of natural gas at temperatures typically in excess of 1000°C. As the natural gas is cracked, the solid carbon produced is deposited on the bubble interface with the gaseous hydrogen encapsulated within. When the bubble reaches the molten metal-air interface at the top of the reactor, the three chemical species (the molten metal, the hydrogen, and the carbon black) are gravitationally segregated due to their density contrasts. The carbon is deposited as residue layers on the top surface of the molten metal and the hydrogen rises to the top of the reactor where it exits through an outlet.36 Since the carbon species involved in the production of turquoise hydrogen are never oxidized or combusted, no carbon dioxide is released despite the inherent use of fossil fuels.

Orange hydrogen37 is also a promising carbon-neutral fuel, although it is yet to be commercialized. Captured CO2 is dissolved into brine to form an enriched-CO2 brine injectant. This injectant is then pumped into specific reactive formations that contain iron (II) oxide minerals such as magnetite and wüstite. These minerals and the injectant react in a standard redox reaction to form hematite, which gives the rock a brick orange color. Hydrogen gas evolves as a by-product and in-so doing forms an artificial reservoir of hydrogen. Simultaneously, the CO2 in the brine reacts with magnesium and calcium oxide to precipitate new carbonate minerals such as calcite and magnetite that has two additional benefits for greenhouse gas emission reductions: (1) the CO2 has been converted from a gaseous form to a solid immobile form and (2) that solid precipitation clogs the pore space and reduces permeability which further traps the CO2. Other production routes that use the conditions within the Earth’s crust together with external engineering techniques such as microbial engineering of processes such as fermentation and nitrogen fixation have been proposed to essentially grow hydrogen reservoirs in the subsurface.38

The role of hydrogen in the energy transition and reaching global net-zero carbon emission ambitions is not only as a form of long-term mass energy storage. Hydrogen is also a crucial feedstock in various industrial processes, particularly in sectors such as refining,39 ammonia production,40 and steel manufacturing.41 By replacing fossil-derived hydrogen with green hydrogen produced from renewable sources, these industries can significantly reduce their carbon footprint. Furthermore, for transportation, hydrogen fuel cells can power vehicles,11 offering an alternative to traditional internal combustion engines and battery electric vehicles. Hydrogen fuel cell vehicles (FCVs) offer fast refueling and longer ranges compared to battery electric vehicles,42 making them suitable for long-haul transportation, heavy-duty vehicles, and applications where weight and space are constraints. Hydrogen can also be used directly in fuel cells to generate heat and electricity for residential, commercial, and industrial applications.43 This can help decarbonize heating systems in buildings, which rely heavily on natural gas, and provide backup power during grid outages. In order for hydrogen to accomplish its role as a major player in the energy transition, key challenges must still be overcome and a hydrogen economy must be established.44 The development of infrastructure for the production, distribution, and storage of hydrogen must be scaled up. This includes building electrolyzers for hydrogen production,45 expanding refueling stations for hydrogen vehicles,46 retrofitting existing pipelines for hydrogen transport,47 and repurposing a number of depleted hydrocarbon fields for geological hydrogen storage.4 Countries and organizations are collaborating on research, development, and deployment initiatives to accelerate the adoption of hydrogen as a clean energy solution.

Although hydrogen is attractive due to its large-scale mass storage potential, no experience exists of the injection and storage of a fluid comprised of pure hydrogen into a porous medium such as an aquifer, depleted gas field or depleted oil field.4,12 Therefore, numerical simulation of the injection, storage, and withdrawal of gaseous hydrogen is imperative in order to assess potential geomechanical consequences. These phenomena include induced seismicity,4,48,49,50 caprock integrity,51,52,53,54 and particularly for onshore storage, surface subsidence.12,51,55 Given the very-low viscosity of hydrogen (9–12 μ Pa s) and its ultra-low density (10–30 kg/m3),4,12 the problem is fundamentally different from the injection and storage of other geo-energy fluids such as supercritical carbon dioxide (190–850 kg/m3, 20–75 μ Pa s),4,12,56 wastewater (1070–1160 kg/m3, 400–1200 μ Pa s),4,57 hydraulic fracturing fluid (2000–1000000 μ Pa s),58 or natural gas (50–210 kg/m3, 15–27 μ Pa s).4,12 Furthermore, since the hydrogen must also be withdrawn, unlike in carbon capture and sequestration, the reservoir, caprock and faults experience cyclical stresses and fatigue,59 that also make hydrogen storage in porous media a unique geo-energy problem.

In order to ensure the commercial viability of future geological hydrogen storage, two criteria must be met: (1) the geomechanical integrity of the storage unit must be maintained, and (2) a public “license to operate” must exist. If injecting hydrogen leads to induced seismicity, pre-existing faults could be reactivated which could offset the caprock or well casing and cement could be damaged, allowing hydrogen to leak. Furthermore, injecting hydrogen could lead to the growth of fractures in the caprock and permeability alteration that enhances flow through the caprock and reduces caprock breakthrough pressure, also enabling hydrogen to leak. Onshore, hydrogen injection may lead to changes in surface elevation, which could potentially lead to damage to commercial and residential property. Thus, geomechanical analyses of these problems are essential in addressing the physical and commercial viability of the technology. In the United Kingdom, injection-induced earthquakes of local magnitude (ML) as low as 1.6 and 2.960 in 2018 and 2019 respectively, occurred at the Preston New Road shale gas fields in Lancashire, and resulted in a nationwide hydraulic fracturing ban from 2019 to 2022. Therefore, investigations into the magnitude of potential induced seismicity during hydrogen storage are critical to prevent H2 storage technology following the same fate as hydraulic fracturing.

A linkage between subsurface fluid injection and induced seismicity has been known since the early 1960’s. In 1961, a 12,000-foot well was drilled at the Rocky Mountain Arsenal Chemical Weapons facility in Denver, United States for the injection of hazardous wastewater into the subsurface.61 Within a few months of the onset of injection, tremors began and over 700 small to modest magnitude earthquakes occurred between 1962 and 1966.61 The U.S. Army later shut in the wellbore in 1966,62 but earthquakes continued to occur, including one of magnitude 4.8 in 1967,63 as the pressure front propagated toward other faults. Earthquakes continued until 1981 according to reports from local residents. In 1969, the United States Geological Survey (USGS) performed injection experiments in Chevron’s Rangely oil reservoir in Colorado.63 They found that earthquakes ceased when the injection rate was reduced below a critical level and were initiated when the injection rate exceeded this critical level. In the early 2000’s, a boom in shale gas exploration occurred.64 With this came the need to dispose of large quantities of waste brine. In 2008, in the proximity of the Barnett Shale formation in Texas, United States, a cluster of earthquakes in the nearby town of Cleburne was linked to the wastewater disposal. Seismometer readings showed that 180 small earthquakes occurred between late-October 2008 and late-May 2009.65 Earthquakes also resulted from operations south of Fort-Worth in June 2009,66 Guy and Greenbrier in Arkansas67 and Youngstown, Ohio in March 2011.68 It was clear that earthquakes were occurring in regions that are not prone to seismicity, and rates of seismicity were drastically enhanced in the years following injection procedures relative to the years prior. The injection of supercritical carbon dioxide for carbon, capture and storage resulted in a sequence of microseismic events at the In Salah site in Algeria in 2004.69 On November 6th 2011, wastewater injection in Oklahoma, the United States lead to the largest ever induced seismic event with a magnitude of 5.6 and known as the “M5.6 Prague earthquake”.70 Significant damage resulted with the collapse of building walls and chimneys and people were injured,71 showing that it is pivotal to understand the physical causes of induced seismicity and how it can be reduced or avoided.

Burtonshaw et al.4 conducted a thorough 3-D, fully coupled hydromechanical investigation that studied the effects of reservoir mechanical properties on fault slip during the multi-year injection, storage, and withdrawal of gaseous hydrogen in a model of a North-sea depleted gas field. In their model, the slip along a 2000 m × 198.5 m fault located within the injection (reservoir) unit was studied for an injection rate of 20 m3/s which was representative of a worst-case scenario. They found that the maximum potential event magnitude was 3.56 and that the reservoir mechanical properties had a strong effect on the induced fault slip (i.e., when the reservoir Young’s modulus was 15 GPa, a magnitude 3.56 event was induced, but when the reservoir Young’s modulus was 40 GPa, the fault was completely stable with no slip). Therefore, moderate induced seismicity was found to be possible in reservoir rocks with certain mechanical properties. Furthermore, it was concluded that reservoir rocks of high Young’s modulus (>40 GPa), high Poisson’s ratio (>0.30) and high Biot coefficient (>0.65) would be ideal for hydrogen storage. Meng et al.72 and Hui et al.73 also conducted smaller analyses of the influence of reservoir mechanical properties (Poisson’s ratio, Young’s modulus, and Biot coefficient) on reservoir fault slip. However, Hui et al.73 used indirect indicators of induced seismicity such as Coulomb failure stress (CFS) rather than direct fault slip. Fan et al.74 performed the only work found in the literature that directly varied the reservoir mechanical properties and assessed the impact on induced seismicity along a fault located in the underburden. They only performed one parametric study of the reservoir Biot coefficient and that was only inclusive of two cases: α = 0.44 and α = 0.79. However, no previous study has examined induced seismicity along underburden faults during hydrogen storage and no previous study has attempted to quantify the underburden mechanical properties which are optimal for reducing fault slip and seismic risk. This is critical to study because some of the largest induced seismic events have occurred along basement faults (as have the majority of induced events75), with small pressure changes of less than 0.2 MPa, such as near magnitude 5 and 6 seismicity in the basement in Milan, Kansas and below the Arbuckle formation in North-Oklahoma respectively, both where wastewater was being injected.76,77

Induced seismicity on basement faults due to the injection of a fluid in an overlying reservoir unit is a well-known phenomenon. Examples of induced seismicity along basement faults include those from wastewater injection (Hearn et al.77; Horton78; Zhang et al., 201379; Verdon,80 Hornbach et al.81; McGarr et al.82), CO2 injection (Cao et al.,69 Chang et al.83), and geothermal (Buijze et al.84; Megies and Wassermann,85 Küperkoch et al.86). When injecting a fluid, the pore pressure increases, which reduces the effective stress acting on the fault, enabling the fault to slip at stress conditions lower than it otherwise would.87,88,89,90,91,92,93 However, for faults in basement units, the low permeability of basement/underburden rocks means that the elevated pore pressure front often remains restricted to the overlying reservoir unit. Therefore, the most common mechanism for fault slip along basement faults during injection is the formation of localized aureoles of elevated shear stress, which form due to poroelastic stress and strain transfer beyond the portions of rock where the injectant pressure front has diffused, and which lead to failure.74,94,95,96,97 Sometimes, low permeability basement faults can slip due to the elevated pore pressure front, if the basement fault also propagates into the reservoir unit, thus allowing direct fluid communication between the reservoir and underburden.74,79,81,98,99

In the present paper, the Imperial College Geomechanics Toolkit (ICGT)100,101 is used to study fault slip and induced seismicity along a single fault located in the underburden unit (below the reservoir unit in which injection, storage, and withdrawal takes place) of a four-layer 4000×4000× 3300 m anticlinal depleted gas field, during hydrogen storage. Simulations cover a period of three years, in which each annual cycle consist of a 5-month injection period, 2-month storage period, and 5-month withdrawal period. A parametric analysis is undertaken for four underburden mechanical properties: Young’s modulus (5–40 GPa), Poisson’s ratio (0.10–0.35), Biot coefficient (0.60–1.00), and matrix permeability (7.0×10−24 m2–1.0×10−16 m2). These parameter ranges are designed to cover all realistic values for the underburden of various depleted North Sea candidate fields found in the literature.4 For each mechanical property, a total of four cases are run, in which only that property is varied. In each case, average fault slip versus time is shown for a time period of three years. This analysis particularly helps to understand two fundamental questions relating to the geological storage of hydrogen in depleted fields: (1) to what extent are underburden faults prone to slip and induce seismicity? and (2) what underburden mechanical properties would be most advantageous and disadvantageous in our pursuit to minimize seismic events during injection, storage, and production?

Model

Experimental setup

Description

A numerical model with dimensions 4000 m × 4000 m × 3300 m (Figure 1) was constructed to study the injection, storage, and withdrawal of gaseous hydrogen. The model is similar to that used in Burtonshaw et al.,4 except that there is no fault in the reservoir unit, and instead, a fault is emplaced in the underburden unit. The model represents a generalized North Sea reservoir, with typical lithologies (chalk overburden, shale caprock, sandstone reservoir, and shale underburden) and geometry (anticlinal structure, with the reservoir bound from above and below by a low permeability seal). The mechanical and fluid properties (Table 2) are prescribed based on a literature search of depleted North Sea oil and gas fields4—the majority of the parameters come from the Leman, Britannia, and Goldeneye fields for which data were available. The fault properties are generalized.Figure 1 Depleted gas field model

(A) A solid and wireframe view of the model showing the four rock layers (overburden is cream, caprock is light gray, reservoir is orange, underburden is dark gray).

(B) A perspective front-right solid and wireframe view of the fault (blue) contained with the reservoir unit and the injection well on the left of the domain (orange) and the withdrawal well on the right of the domain (red). The reservoir and underburden units are removed in order to observe these.

(C) A top view of the rectangular fault surface.

Table 2 The geological, mechanical, geometrical, fluid dynamic, and fault properties and the in situ stresses of the four different rock layers within the models

Property	Overburden	Caprock	Reservoir	Underburden	
Mechanical	
	
Lithology	Chalk	Shale	Sandstone	Shale/Coal	
k (m2)	3.0×10−15	1.0×10−21	1.0×10−13	1.0×10−21	
ϕ	0.30	0.05	0.15	0.05	
h (m)	2300–2400	100	150	650–750	
E (Pa)	4.5×109	1.5×1010	1.5×1010	4.0×1010	
ν	0.35	0.35	0.15	0.35	
α	0.87	0.93	0.35	0.93	
ρb (kg/m3)	2080	2500	2500	2500	
σ∗ (Pa)	1.0×106	1.3×106	2.0×106	1.3×106	
	
Fluid	
	
Fluid	Brine	Brine	Hydrogen	Brine	
ρf (kg/m3)	1058	1042	9	1037	
μ (μ Pa s)	561.0	357.5	10.4	296.0	
c (Pa−1)	4.47×10−10	4.57×10−10	6.67×10−8	4.68×10−10	
	
Stresses	
	
σh (MPa)	–	–	58.5	–	
σH (MPa)	–	–	58.5	–	
	
Fault	
	
L (m)	–	–	–	2000	
w (m)	–	–	–	1000	
af,r (m)	–	–	–	1.0×10−6	
θ (°)	–	–	–	53.1	
ψ(°)	–	–	–	0	
x (m)	–	–	–	(0,0,2475)	
The symbology in the table is as follows: k = Permeability, ϕ = Porosity, h = Thickness, E = Young’s modulus, ν = Poisson’s ratio, α = Biot coefficient, ρb = Rock Density, K = Bulk Modulus, σ∗ - Uniaxial Tensile Strength, ρf = Fluid Density, μ = Fluid Viscosity, c = Fluid Compressibility, σh = minimum horizontal in situ stress at base, σH = maximum horizontal in situ stress at base, L = fault length, w = fault width, af,r = residual aperture, θ = fault dip, ψ = fault orientation to North and x = fault centroid location (x,y,z). The properties herein are derived from the literature review (Table 4) in Burtonshaw et al.4

Gaseous hydrogen is injected at 40.0 m3/s into the reservoir through a horizontal wellbore that penetrates the whole x-direction of the model and is produced through a second wellbore of the same dimensions at 34.0 m3/s at the opposite end (y-direction) of the model. Both wellbores are at a depth of 2,475 m. Simulations are run for three annual cycles with each annual cycle consisting of a 5-month injection period, a 2-month storage period, and a 5-month withdrawal period. The model contains a solitary fault of dimensions 2000 m × 1000 m, that is centered at a depth of 2,950 m, and is fully encapsulated within the low-permeability underburden unit. The fault strikes parallel to the injection and withdrawal wellbores.

The boundary conditions for the simulation are defined as follows: In situ stresses are imposed as Neumann boundary conditions. Initially, matrix pore pressure is in equilibrium with mechanical stresses, and a series of hydrostatic steps are conducted to establish equilibrium before the injection phase commences. Injection and withdrawal processes are represented as volumetric fluxes applied at the respective well elements. Modeling injection and withdrawal as volumetric fluxes at discrete well elements assumes a simplified representation of fluid flow within the reservoir because it presumes homogeneity in fluid distribution and pressure propagation within the reservoir because of the assumed homogeneity in rock properties such as reservoir permeability. Furthermore, in reality wellbore perforations are never perfectly spaced and of identical size. Therefore, fluid delivery to the rock is unlikely to be perfectly distributed along the well length. This may lead to slightly different pressure distributions in the reservoir which may slightly influence stress and strain transfer from the reservoir into the underburden, leading to small differences in the magnitude and timing of fault slip.

Discretization

The domain is spatially discretized with isoparametric linear tetrahedra internally and with isoparametric linear triangles on the model boundaries and fault (Figure 1). The model contains 15,428 triangular surface elements, of which 3898 are on each fault face, and 238,281 tetrahedral elements in the volume.

The fault exists as two rectangular surfaces that are connected at the fault tips, making a negative volume within the mesh. The slip is computed at each node on the fault surfaces and a weighted-average is used to compute the average fault slip used in the results section:(Equation 1) d=∑n=1Nxn·AnAT

where d is the average fault slip at a given time step, xn is the slip at the node of interest, n, An is the area surrounding the node of interest and AT is the total fault surface area.

Numerical experiments

Four sets of simulations have been conducted in which four simulations are run per set. Each set corresponds to a particular mechanical property (i.e., Young’s modulus, Poisson’s ratio, and matrix permeability) of the underburden unit. Within each set of simulations, the property is parametrically varied in four separate simulations, keeping all other parameters constant as in the default case, in order to assess the influence of that underburden property on induced fault slip. Table 3 summarizes the test cases performed.Table 3 Mechanical properties of the different simulation cases

Simulation	Default Case	Case 2	Case 3	Case 4	
Eres (GPa)	40	5	15	30	
νund	0.35	0.10	0.25	0.30	
αund	0.93	0.60	0.80	1.00	
kund (ND)	1	0.007	100	100,000	
E = Young’s modulus, kund = Permeability of the matrix - underburden unit, α = Biot coefficient, and ν = Poisson’s ratio. The subscript “und” refers to the underburden unit.

Mesh and time step effects

In order to ensure that the fault slip solution as a function of time is not dependent on the refinement of the meshes used or the time step used, an extensive study of the effects of the bulk rock domain and the fault mesh refinements, and the time steps, are performed in the following sections.

Domain mesh refinement

A total of four simulations were performed to assess the impact of the bulk rock domain mesh refinement on the average fault slip as a function of time. The four cases are summarized in Table 4, and vary from a coarse case with approximately just 17,000 nodes and 95,000 elements to a very refined case with approximately 125,000 nodes and 750,000 elements. In cases 1 and 2, the entire domain is refined evenly, in case 3, a box surrounding the fault of dimensions 2500×1050 × 750 m is refined more than the remaining domain and in case 4, the same box is refined to a greater extent than in case 3. The domain meshes are shown in Figure 2A. The results are shown in Figure 3B. It is found that provided the bulk rock mass is discretized with approximately at least 43,000 nodes and 253,000 elements then the fault slip solution is typically the same to within 0.1 mm (100 microns). Case 1 was found to be too coarse to accurately yield the solution with very large numerical errors generated. Thus, for the remaining simulations in this work the bulk rock mass will be discretized according to case 2, which suffers from less numerical errors than cases 3 and 4, and is computationally less expensive.Table 4 Four simulation cases run to examine the effect of domain refinement on the slip solution

Simulation Case	Nodes	Elements	Fault Elements	
Case 1	16,532	94,943	3826	
Case 2	42,856	253,709	3898	
Case 3	51,664	307,095	3827	
Case 4	125,339	748,302	3865	

Figure 2 Comparison of different domain and fault meshes

(A) Wireframe views of the four versions of the model for each refinement case in Table 4. The dashed boxes denote the zones of increased refinement.

(B) Zoomed views of the mesh at the center of the underburden unit.

(C) Mesh of the fault surface for the five versions of the fault in the model for each refinement case in Table 5, for varying amounts of fracture elements Efrac.

Figure 3 Plots of average slip as a function of operation time for different fault refinements, domain refinements and time steps

(A) five cases of different fault refinements represented by different numbers of fault elements, (B) four cases of different bulk rock/domain refinements represented by different numbers of nodes and (C) four cases of different time step. The full discretization data for each case in (A) and (B) is reported in Tables 4 and 5 respectively.

Fault mesh refinement

A total of five simulations were performed to assess the impact of the fault mesh refinement on the average fault slip as a function of time. The five cases are summarized in Table 5, and vary from a coarse case with approximately 700 elements to a very refined case with approximately 6,250 elements. The fault meshes are shown in Figure 2C. The results are shown in Figure 3A. It is found that provided the fault is discretized by at least approximately 1500 elements then the fault slip solution is the same to within a maximum of 0.3 mm (300 microns). Only in the coarsest case with 706 fault elements was the fault slip solution significantly different with up to approximately 1.5 mm difference. Thus, for the remaining simulations in this work the fault will be discretized according to case 4, although cases 2 and 3 would have also been acceptable.Table 5 Five simulation cases run to examine the effect of fault refinement on the slip solution

Simulation Case	Nodes	Elements	Fault Elements	
Case 1	17,603	105,927	706	
Case 2	61,508	369,113	1526	
Case 3	23,363	138,125	2594	
Case 4	42,855	253,697	4028	
Case 5	61,445	363,785	6252	

Time step

A total of four simulations were performed to assess the impact of the time step on the average fault slip as a function of time. The time steps tested included 0.2628×107 s, 0.1314×107 s, 0.657×106 s, and 0.3285×107 s. Time steps smaller than this were not tested due to the computational constraint of realistic simulation run time. The results are shown in Figure 3C. It is found that provided a time step of at least 0.657×106 s is used, the fault slip solution converges to the same result. Time steps slightly larger than 0.657×106 s replicate the same trend in slip but there are large inaccuracies such as in the 0.1314×107 s case, whereas time steps much larger than 0.657×106 failed to replicate the same trend in slip and were extremely inaccurate. Given time steps of 0.657×106 s or lower yield the same result, the 0.657×106 s case is chosen for the simulations performed in this work in order to reduce computational time.

Results

This section outlines the results from the numerical experiments outlined in the STAR Methods section. It is found in every case that the fault in the underburden unit does not substantially slip during the bulk of the injection, storage, or withdrawal periods. However, the fault slips significantly at the cross-over points between operation regimes (i.e., injection, storage, and withdrawal) and at the very onset of injection at time t = 0, reflecting sudden changes in the stress field at these points. Only at these points, where slip is induced suddenly may the slip potentially be seismic. These findings are unsurprising given that abrupt increases in injection rates tend to shortly precede the occurrence of earthquakes.68,102,103,104 Alghannam and Juanes102 found that the likelihood of triggering earthquakes depends critically on the rate of increase in pore pressure and that fluid injection at constant rate acts to trigger seismic rupture at early times before then observing aseismic creep at late times. This supports the present findings, wherein the fault slips to a relatively large extent immediately following the onset of injection, due to the large injection volume of 40 m3/s, and then slips very minimally for the remainder of the injection periods. For their case of injection at constant rate, it was observed that the critical stiffness increased at early times, potentially triggering earthquakes, and decreased at late times, potentially resulting in the cessation of earthquakes. Additionally, a dramatic drop in critical stiffness was observed upon stopping injection, showing that the propensity for seismicity rapidly follows changes in injection rate.

Young’s modulus

The Young’s modulus of the shale underburden was varied between 5 GPa and 40 GPa, and the influence on average slip as a function of time is shown in Figure 4A. It is generally found that as Young’s modulus increases, the fault slip decreases. The initial slip of the fault at t = 0, at the onset of injection, for the 5, 15, 30, and 40 GPa cases is approximately 3.83 mm, 1.31 mm, 1.48 mm, and 1.20 mm, respectively. For the slip period at approximately 1.58 years at the transition between storage and withdrawal regimes, the fault slip for the 5, 15, 30, and 40 GPa cases is approximately 2.61 mm, 1.61 mm, 1.05 mm, and 0.88 mm, respectively.Figure 4 Average slip as a function of operation time for four cases of varying underburden properties

(A) Young’s modulus, (B) Poisson’s ratio, (C) Biot coefficient, and (D) permeability.

The fault aperture and the fault slip at each individual node on the fault surface at a time of 1.58 years is shown in Figures 5 and 6. No systematic change in aperture was found with increasing Young’s modulus, although generally it is seen that larger apertures are observed for smaller Young’s modulii. This broad change can be explained by the fact that a stiffer rock is less deformable and more resistant to changes in shape or size when subjected to stress. As a result, the rock adjacent to the fault experiences less deformation under the same stress conditions compared to rock with lower Young’s modulus, and this reduced displacement and contraction of the neighboring rock restricts the ability of the fault to widen, leading to smaller fault apertures. The trend in aperture is not perfectly systematic because a combination of other factors also influences the aperture. For example, the interaction between flow and mechanical deformation plays a crucial role in fault behavior. Changes in Young’s modulus can influence fluid flow patterns and pore pressure distribution, which in turn affect the mechanical response of the rock mass and ultimately the fault aperture.Figure 5 The apertures of the fault at a time of 1.58 years during the end of the second storage and beginning of the second withdrawal periods for cases of varying underburden properties

(A–D) Young’s modulus, (E–H) Poisson’s ratio, (I–L) Biot coefficient, and (M–P) permeability.

Figure 6 The individual slip at each node of the fault at a time of 1.58 years during the end of the second storage and beginning of the second withdrawal periods for cases of varying underburden properties

(A–D) Young’s modulus, (E–H) Poisson’s ratio, (I–L) Biot coefficient and (M–P) permeability. Note that the scales are not the same in each subfigure.

Effective stress and strain are shown in Figures 7 and 8, respectively. It is found that as Young’s modulus decreases, the strain in the underburden unit increases significantly. A rock with a lower Young’s modulus is more compliant, meaning it is easier to deform and has less resistance to applied forces. As a result, when subjected to stress, such a rock exhibits greater strain because it is more willing to undergo deformation. Furthermore, the vertical stress profile through the model shows that at the depths of the upper underburden (2650–2950 m), the stress is the lowest in the 5 GPa case, where the most fault slip was observed and the stress is the highest in the 40 GPa case, where the least slip was observed, due to the fact that the high strain in the 5 GPa case results in larger stress relative to the very low strain in the 40 GPa case. The lower the effective normal stresses on the fault and the higher the strains in the vicinity of the fault, the greater the propensity the fault has to slip, which is observed here.Figure 7 Plots of effective stress vertically through the center of the model in the z-direction for four cases of varying underburden properties

(A) Young’s modulus, (B) Poisson’s ratio, (C) Biot coefficient, and (D) permeability at a time of 1.58 years.

Figure 8 The strain over a vertical interval containing the anticlinal caprock, reservoir and underburden units at a time of 1.58 years during the end of the second storage and beginning of the second withdrawal periods for cases of varying reservoir properties

(A–D) Young’s modulus, (E–H) Poisson’s ratio, (I–L) Biot coefficient and (M–P) permeability.

Taking into account Equations 12 and 13, and making the assumption that all the slip in one time step is released instantaneously, the upper bound for the moment magnitude of the potential induced seismicity at the onset of the initial injection is 2.73, 2.74, 2.98, and 3.00 for the 5, 15, 30, and 40 GPa cases, respectively. Therefore, despite fault slip increasing as Young’s modulus decreases, the magnitude of the maximum potential induced seismicity decreases as Young’s modulus decreases. This is because as Young’s modulus decreases, the contribution of the smaller shear modulus, G, given by E/2(1+ ν) in Equation 12 outweighs the increasing fault slip, d, and thus the seismic moment release is smaller as Young’s modulus decreases. It is found that the Young’s modulus of the underburden can have a strong influence on the magnitude of the induced seismicity.

Poisson’s ratio

The Poisson’s ratio of the shale underburden was varied between 0.10 and 0.35, and the influence on average slip as a function of time is shown in Figure 4B. It is found that as underburden Poisson’s ratio increases, the fault slip generally decreases. The initial slip of the fault at t = 0, at the onset of injection, for the 0.10, 0.25, 0.30, and 0.35 cases is approximately 1.33 mm, 1.16 mm, 1.15 mm, and 1.20 mm, respectively. For the slip period at approximately 1.58 years at the transition between storage and withdrawal regimes, the fault slip for the 0.10, 0.25, 0.30, and 0.35 cases is approximately 0.92 mm, 0.91 mm, 0.90 mm, and 0.88 mm, respectively. Thus, the Poisson’s ratio of the underburden unit does not have a large influence on basement fault slip with the slip difference being a maximum of 125 microns between the cases.

The fault aperture and the fault slip at each individual node on the fault surface are shown in Figures 5 and 6. No systematic change in aperture is found with increasing Poisson’s ratio. Poisson’s ratio characterizes the ratio of lateral contraction to longitudinal extension in a rock when it is subjected to stress. Changes in Poisson’s ratio affect the distribution of stress within the rock, but it does not necessarily impact fault behavior directly because the stress distribution around the fault remains relatively symmetrical with changes in Poisson’s ratio. This results in relatively consistent fault apertures irrespective of Poisson’s ratio.

Effective stress and strain are showed in Figures 7 and 8, respectively. The difference in underburden strain between the cases is minor, which is reflected in the near-identical fault slip observed in the four cases. The vertical stress profile through the model shows that at the upper underburden (2650–3300 m), the stress is the highest in the highest Poisson’s ratio case of 0.35, followed by the 0.30 case, with lower stresses seen in the 0.10 and 0.25 cases, implying a greater propensity for the fault to slip in the lower Poisson’s ratio cases, which is what is observed. The broadly larger stress with increased Poisson’s ratio occurs because a higher Poisson’s ratio leads to greater lateral expansion of the underburden rock in response to the applied stress. This lateral expansion results in increased stress concentrations in the upper underburden because the underburden rock in these regions is accommodating the lateral expansion while still being subjected to the compressive stress in the vertical direction and the opposing horizontal in situ stresses.

Taking into account Equations 12 and 13, and making the assumption that all the slip in one time step is released instantaneously, the upper bound for the moment magnitude of the potential induced seismicity at the onset of the initial injection is 3.09, 3.01, 3.00, and 3.00 for the 0.10, 0.25, 0.30, and 0.35 cases, respectively. Therefore, just as with the trend in fault slip, the magnitude of the maximum potential seismicity decreases as Poisson’s ratio increases. The magnitude of the maximum potential seismicity is largely insensitive to changes in underburden Poisson’s ratio, but is slightly larger when the Poisson’s ratio was low (0.10).

Biot coefficient

The Biot coefficient of the shale underburden is varied between 0.60 and 1.00, and the influence on average slip as a function of time is shown in Figure 4C. It is found that Biot coefficient is not a strong control on fault slip, but that a small trend is observed in which fault slip increases as Biot coefficient increases. The initial slip of the fault at t = 0, at the onset of injection, for the 0.60, 0.80, 0.93, and 1.00 cases is approximately 1.15 mm, 1.18 mm, 1.20 mm, and 1.22 mm, respectively. For the slip period at approximately 1.58 years at the transition between storage and withdrawal regimes, the fault slip for the 0.60, 0.80, 0.93, and 1.00 cases is approximately 0.84 mm, 0.86 mm, 0.88 mm, and 0.89 mm, respectively. Thus, Biot coefficient does not have a strong influence on fault slip along basement faults for underburden rocks in the typical range of 0.6–1.0.

The fault aperture and the fault slip at each individual node on the fault surface at a time of 1.58 years is shown in Figures 5 and 6. No systematic change in aperture is found with increasing Biot coefficient. The Biot coefficient governs the extent to which changes in pore fluid pressure influence solid deformation and vice versa. Since the underburden has a very low permeability, fluid flow is somewhat decoupled from solid deformation in this unit at the timescale of the presented simulations, because pore fluids are unable to move freely through the rock mass and exert significant pressure changes. In the absence of significant fluid flow and pressure changes, mechanical deformation processes become the dominant mechanism controlling fault slip behavior. Variations in Biot coefficient therefore have minimal impact on fault behavior in terms of slip and aperture because pressure change within and near the fault is inhibited.

Effective stress and strain are showed in Figures 7 and 8, respectively. The vertical stress profile through the model shows that throughout the upper underburden (2650–2950 m), the stress decreases with increasing Biot coefficient. This is because higher Biot coefficients lead to greater porosity changes and volume dilation, which can reduce stress within the rock mass. The lower the effective normal stresses on the fault, the greater propensity the fault has to slip which is observed here. The strains in the underburden unit are similar in each case, although in the region directly around the fault, the strain is slightly larger in the 0.60 and 0.93 cases, relative to the 0.80 and 1.00 cases. It is also seen that the aperture is larger in the 0.60 and 0.93 cases, relative to the 0.80 and 1.00 cases, showing that the aperture is a function of the strain.

Taking into account Equations 12 and 13, and making the assumption that all the slip in one time step is released instantaneously, the upper bound for the moment magnitude of the potential induced seismicity at the onset of the initial injection is 2.99, 3.00, 3.00, and 3.01 for the 0.60, 0.80, 0.93, and 1.00 cases, respectively. Therefore, the magnitude of the maximum potential seismicity is insensitive to changes in underburden Biot coefficient.

Matrix permeability

The matrix permeability of the shale underburden is varied between 1×10−16 m2 and 7×10−24 m2, and the influence on average slip as a function of time is shown in Figure 4D. It is found that for permeabilities in the ranges 7×10−24 m2 to 1×10−19 m2, the fault slip is insensitive to changes in underburden permeability. However, when the permeability of the underburden unit is high (i.e., 1×10−16 m2), the initial slip of the fault at t = 0 is considerably larger with 1.66 mm slip, compared to 1.25 mm, 1.20 mm, and 1.06 mm in the 7×10−24 m2, 1×10−21 m2, and 1×10−19 m2 cases. Thus, in particular early seismic hazard increases along underburden faults when the underburden permeability is large. This is because with higher underburden permeability (1×10−16 m2), the rock mass allows for more efficient fluid flow from the reservoir to the underburden. As a result, pore fluids can migrate more rapidly into and within the underburden, leading to a buildup of pore pressure within the fault zone. For the slip period at approximately 1.58 years at the transition between injection and storage regimes, the fault slip for the 7×10−24 m2, 1×10−21 m2, 1×10−19 m−21, and 1×10−16 m2 cases is approximately 0.88 mm, 0.88 mm, 0.90 mm, and 0.82 mm.

Effective stress and strain are showed in Figures 7 and 8, respectively. The average strain in the underburden unit increases as permeability decreases. Furthermore, the strain at the top of the underburden unit increases significantly as permeability decreases. Lower permeabilities are associated with higher pore pressures within the underburden unit that can exert additional support against the applied stresses, reducing the effective stress within the rock mass. This reduction in effective stress enhances the potential for mechanical deformation, resulting in increased strain within the underburden unit. However, in the region of the fault (center of the underburden unit), the strain is very similar in all cases, supporting the fact that the permeability has no large effect on the induced fault slip after the initial slip at t = 0.

Taking into account Equations 12 and 13, and making the assumption that all the slip in one time step is released instantaneously, the upper bound for the moment magnitude of the potential induced seismicity at the onset of the initial injection is 3.01, 3.00, 2.96 and 3.09 for the 7×10−24 m2, 1×10−21 m2, 1×10−19 m2, and 1×10−16 m2 cases, respectively. Therefore, the magnitude of the maximum potential seismicity is largely insensitive to changes in underburden permeability when the permeability is between 10−24 and 10−19 m2, but can increase for larger permeabilities.

Discussion

By investigating the influence of the mechanical properties of the various lithological units within storage fields (during hydrogen injection, storage and withdrawal) on fault slip and potential induced seismicity, it is possible to understand ahead of the onset of hydrogen storage operations, which storage fields are most and least prone to triggering induced seismicity. Therefore, this work performed a thorough analysis of four underburden mechanical properties on fault slip and induced seismicity along a large-scale underburden fault that is hydraulically isolated from the overlying reservoir.

As underburden Young’s modulus increases, the average fault slip generally decreases. For Young’s modulii in the range 15 GPa–40 GPa, the fault slips were similar, but the fault slip increased very significantly when the underburden Young’s modulus was reduced to 5 GPa, with the 5 GPa case yielding over three times the fault slip of the 40 GPa case. However, the magnitude of the maximum potential induced seismicity does not correlate with the trend in fault slip. As Young’s modulus decreased, the magnitude of the maximum potential induced seismicity decreased due to decreasing shear modulus with decreasing Young’s modulus. The reduction in maximum magnitude was significant with an Mw 3.00 in the 40 GPa case decreasing to an Mw 2.73 in the 5 GPa case. The strain induced in the reservoir unit decreased with increasing Young’s modulus.

No published results exist for comparison that studies the influence of underburden Young’s modulus on induced underburden fault slip. For a case of hydrogen storage, Burtonshaw et al.4 varied the reservoir Young’s modulus between 10 GPa and 60 GPa and found that as Young’s modulus increased, the average reservoir fault slip decreased very significantly, with over 10 cm slip in the 10 GPa case and no slip when the reservoir Young’s modulus was 40 GPa or more. However, they found that as Young’s modulus increased, the magnitude of the maximum potential induced seismicity decreased because the increasing shear modulus as Young’s modulus increased was not sufficient to outweigh the lack of any slip. In the hydraulic fracturing work of Meng et al.,72 it was also found that reservoir fault slip decreased very significantly from 12 cm to 1–2 cm as Young’s modulus was increased from 5 to 50 GPa. Therefore, as Young’s modulus increases fault slip decreases for both reservoir and underburden faults. However, the trend with the maximum potential induced magnitude was different depending on whether the fault was in the underburden or reservoir unit, owing to the fact that the slip along the underburden fault is less sensitive to underburden Young’s modulus than the slip along the reservoir fault is to reservoir Young’s modulus, where absolutely no slip was observed for high Young’s modulii (>40 GPa).

As underburden Poisson’s ratio increases, the average fault slip decreases. However, the dependency of the underburden fault slip on Poisson’s ratio is weak, with a maximum slip discrepancy between the Poisson’s ratio cases of less than 0.2 mm (16% of the slip in the 0.10 case). The fault slip was very similar in the 0.10, 0.25, and 0.30 cases but decreased slightly more significantly when the Poisson’s ratio was increased to 0.35, suggesting that the average fault slip might decrease more significantly for Poisson’s ratio in the range 0.35–0.50. The magnitude of the maximum potential induced seismicity followed the same trend as the fault slip: as the Poisson’s ratio increases, the magnitude of the potential seismicity decreases. The magnitude decreased from a maximum of 3.09 in the 0.10 case to 3.00 in the 0.35 case.

No prior results were found that examined the effect of underburden Poisson’s ratio on induced underburden fault slip. Burtonshaw et al.4 varied the reservoir Poisson’s ratio between 0.00 and 0.30 and found that as Poisson’s ratio increases, the average reservoir fault slip increases, with 4.53 cm of slip in the 0.00 case increasing to 10.43 cm in the 0.30 case. However, they found that for the highest Poisson’s ratio case of 0.30, the magnitude of the maximum potential induced seismicity was the smallest, because a large proportion of the fault surface was now in stick. Therefore, as Poisson’s ratio increases, the trend in fault slip is different for both underburden and reservoir faults. Furthermore, the reservoir Young’s modulus exerts a much greater control on the average fault slip than the underburden Young’s modulus. The trend in the maximum magnitude of the potential induced seismicity is the same for both underburden and reservoir faults: as Poisson’s ratio increases, the magnitude decreases.

As underburden Biot coefficient increases, the average fault slip increases. However, the strength of the dependence is weak, with a maximum slip difference between the largest (1.00) and smallest (0.60) Biot coefficient cases of less than 0.1 mm. The moment magnitude of the maximum potential induced seismicity also increased as Biot coefficient increased, although it was largely insensitive to changes in Biot coefficient, with a magnitude of 3.01 in the largest Biot coefficient case (1.00) and 2.99 in the smallest Biot coefficient case (0.60). The effective stress around the fault decreased as Biot coefficient increased, promoting fault slip.

Again no published results were found in the literature that analyzed the effect of underburden Biot coefficient on underburden fault slip. However, for a case of wastewater injection, Fan et al.74 varied reservoir Biot coefficient in two cases (0.44 and 0.79) and computed the change in CFS. They found that as reservoir Biot coefficient decreased, the change in CFS became more positive and thus the likeliness of underburden fault slip increased. Although, just two test cases are likely insufficient to draw solid conclusions. Furthermore, the difference in changed CFS was between the two cases was less than 0.2 MPa, suggesting that the expected slip difference with differing Biot coefficient would be small. Burtonshaw et al.4 varied the reservoir Biot coefficient between 0.25 and 0.90 and found that the relationship between reservoir fault slip and reservoir Biot coefficient was non-systematic, with the magnitude of the fault slip similar for values between 0.25 and 0.65 but no slip occurred when the Biot coefficient was as high as 0.90. Therefore, the fault slip along underburden and reservoir faults responds differently to Biot coefficient depending on whether it is the reservoir or underburden Biot coefficient. The reservoir Biot coefficient has a much stronger influence on reservoir fault slip and event magnitude (varies between 3.56 and no event) than underburden Biot coefficient has on underburden fault slip and event magnitude (varies between 2.99 and 3.01).

In the range 0.007 nD to 100 nD, the underburden fault slip is insensitive to changes in underburden permeability. Even when the underburden permeability increases to as much as 100,000 nD, the underburden fault slip remained insensitive to underburden permeability throughout the bulk of the three-year injection-storage-withdrawal period, only differing in places (i.e., at the onset of the initial injection where the fault slip was 33% higher in the 100,000 nD case relative to the 100 nD case). Therefore, it can be concluded that although underburden fault slip is similar for most typical underburden permeabilities, the lower the underburden permeability, the lower the risk of underburden seismicity.

As before, no previous studies have examined the influence on underburden permeability on underburden fault slip. Burtonshaw et al.4 varied the intrinsic permeability of the reservoir between 5 mD and 800 mD and found broadly that as reservoir permeability decreased, the magnitude of the average fault slip increased. However, the strength of the dependence of the magnitude of the fault slip on the reservoir permeability was weak with only a 6% difference between the whole range of permeabilities. Additionally, they found that the reservoir permeability has an effect on the time in which the slip occurs. Slip was induced earlier for higher reservoir permeabilities due to the fact that the fault pressure at a given time varies depending on the permeability. For smaller permeabilities, the fluid front from the injection well takes longer to reach the fault, and subsequently leak-off into the fault or pressurize the surrounding area. Additionally, even though the magnitude of the average fault slip increased with decreasing permeability, the proportion of the fault area in stick increased as permeability decreased. As a result, it was found that the magnitude of the induced seismicity broadly decreased as reservoir permeability decreased. Therefore, for both underburden and reservoir faults, the magnitude of induced seismicity decreases as the underburden/reservoir permeability decreases. In this work, no temporal delay in the underburden fault slip was seen with changing permeability, like with the reservoir fault cases, because the underburden fault is hydraulically isolated from the injection unit where pore pressure is increasing. Thus, the seismicity is not a direct product of pore pressure change, since the underburden fault never comes into contact with the fluid pressure front and as such should not be temporally affected. In the case where the underburden permeability was extremely high (100,000 nD), some additional slip occurred at the onset because some of the injectant fluid front rapidly leaked into the upper underburden, allowing some minor pressurization of the underburden fault.

In order to assess the effect of the chosen fault orientation on the results, additional simulations were run varying the fault orientation (a 45° rotation and a 240° rotation from the direction of the x axis). It was found that the slip was maximal in the fault orientation case used in this paper (0° from the x axis). When the fault was rotated by 240° from the direction of the x axis, the fault slip was reduced by approximately 17% on average relative to the cases presented herein. Furthermore, when the fault was rotated by 45° from the direction of the x axis, the fault slip was reduced by approximately 46% on average relative to the cases presented herein. This was because strain in the center of the underburden was larger when the fault was oriented in the direction of the x axis (in the plane of the maximum horizontal in situ stress). The strain at the center of the underburden in the 45° orientation case was the lowest, reflecting the fact that the fault orientation is the furthest away from the maximum horizontal in situ principle stress in this case. Therefore, this further reinforces that our computed magnitudes in this work are upper bounds on the maximum potential magnitudes that may occur along basement faults of this size during geological hydrogen storage.

Figure 9 displays a plot of the maximum potential magnitude of the induced seismicity versus fault slip where each point represents the points where rapid and significant slip occurred within one time step from Figure 4. These are maximum potential magnitudes because it has been assumed that all the fault slip within one time step was released instantaneously as in a seismic event. Sensitivity lengths, SMw and Sd, can be defined as Mw,max – Mw,min and dmax – dmin, respectively. The value of the sensitivity lengths will then display which parameters have the greatest influence of fault slip and maximum potential magnitude. SMw for Young’s modulus, Poisson’s ratio, Biot coefficient, and matrix permeability are 0.53, 0.63, 0.50, and 0.61, respectively. Sd for Young’s modulus, Poisson’s ratio, Biot coefficient, and matrix permeability are 3.61, 1.16, 1.03, and 1.46, respectively. Thus, in descending order, underburden fault slip is most sensitive to: (1) Young’s modulus, (2) matrix permeability, (3) Poisson’s ratio, and (4) Biot coefficient. Also in descending order, maximum potential magnitude is most sensitive to: (1) Poisson’s ratio, (2) intrinsic permeability, (3) Biot coefficient, and (4) Young’s modulus.Figure 9 Maximum potential induced magnitude versus average underburden fault slip

Includes all potential events from the Young’s modulus, Poisson’s ratio, Biot coefficient, and permeability cases. Zoomed views of the two main event clusters from the left, are shown on the right.

To minimize seismic risk during underground geological hydrogen storage, the depleted gas fields shortlisted as candidates fields that have low underburden Young’s modulus (<15 GPa), high underburden Poisson’s ratio (>0.25), low underburden Biot coefficient and low underburden intrinsic permeability (<1×10−19 m2) should be selected for storage. Equally, those with high reservoir Young’s modulus (>40 GPa), high reservoir Poisson’s ratio (>0.30), and high Biot coefficient (>0.65).

Hydrogen storage differs from other geo-energy applications in that the adopted injection rates will need to be very high. The ultra-low density of hydrogen of less than 10 kg/m3 (at typical reservoir temperatures and pressures) implies that very large volumes of fluid must be injected into the reservoir rock in order to store the masses necessary to successfully provide a reasonable percentage of the future energy mix. The injection rate used herein of 40 m3/s corresponds to an injection volume of 525,600,000 m3 over the 5-month injection period, which itself corresponds to the period of time in which an energy excess from renewable energy will exist. Given the density of the hydrogen is 9 kg/m3, this corresponds to an injected mass of 4,730,400,000 kg. One kg of hydrogen can produce 33.3 kWh of energy,4,105 which means that in the present model 157,522,320,000 kWh of energy is stored (approximately 157 TWh). This corresponds to just less than 10% of the United Kingdom’s national annual energy demand today.4 Therefore, this injection rate was chosen to be a worst-case scenario, at the upper limit of the hydrogen volumes that may be stored in future geological hydrogen storage projects. This has enabled the maximum potential seismic magnitudes to be quantified.

Limitations of the study

One of the benefits of the ICGT is that it has the capabilities of full hydromechanical coupling and friction, and can fully represent faults and fractures as negative volumes within a three-dimensional rock continuum. Due to this, fault slip can be fully resolved across the whole 2-D fault surface by directly computing node-to-node displacement comparisons of matching pairs of nodes on opposing fault surfaces. Other codes and commercial software used to investigate induced seismicity during a fluid injection, storage, or withdrawal process typically lack this capability, and thus, the literature is dominated by indirect evaluations of seismicity, such as via CFS,73,74 as opposed to direct mechanical fault slips; or direct fault slip is indeed computed but along a 1-D fault within a 2-D model that results in the fault slip being computed at just a single point along the entire depth of the fault.72,106,107,108,109,110

It is assumed in the work herein that each rock layer is homogeneous, but in reality, rock heterogeneities would likely have some influence on slip. Fault geometries and slip distributions across the fault surface are likely to be influenced by heterogeneities due to more complex fault zone architectures and fault thickness variations, as well as the presence of fault breccias and shear bands which impact localized permeability and thus the fault pressurization and subsequently the induced slip. Furthermore, rock stiffness heterogeneities within the fault zone may result in non-uniform stress distributions that alter the magnitude and distribution of stress accumulations during, and stress release after, slip events. It is also assumed that each strata are isotropic but non-isotropic rocks can have different horizontal and vertical permeabilities due to various rock fabrics that may exist within the rock that may enhance or limit fluid flow and pressure propagation in the vicinity of the fault relative to a homogeneous rock. These pore fluid pressure changes surrounding the fault may affect slip behavior by influencing frictional resistance and altering effective stresses. Linear-elastic behavior of the rock matrix is assumed but in reality some minor plastic deformation and damage accumulation may occur along the fault. The model does not include dynamic effects such as transient slip behavior and seismic wave propagation which may lead to potential inaccuracies in predicting the magnitude and timing of fault slip events.

Given the complexity of fully coupling thermal effects with hydromechanical processes in numerical simulations, excluding thermal effects can serve as an initial simplification to isolate the mechanical aspects of fault slip behavior. This approach allows for a clearer and more focused analysis of the mechanical controls on fault behavior, providing valuable insights that can later be integrated into more complex models. It is not expected that the temperature of the reservoir should change significantly during the injection and withdrawal of the hydrogen. However, there may be some minor thermal effects that influence fault slip. For example, hydrogen can exhibit a negative Joule-Thomson effect,111 meaning that contrary to most gases, its temperature increases as it expands. Therefore, it would be expected that during periods of withdrawal the fault slip may be enhanced slightly relative to that predicted here. This is because as hydrogen is withdrawn, the reservoir pressure decreases, the hydrogen expands, and its temperature increases which can enhance pore fluid pressures due to increased fluid volumes within fault zones in a process known as thermal pressurization and therefore reduce effective stress acting on the fault. Contrastingly, during periods of injection the fault slip may be reduced slightly relative to that predicted here due to the fact the hydrogen contracts as the reservoir pressure increases and its temperature decreases, which minimizes the effect of thermal pressurization. Even though hydrogen has a negative Joule-Thomson effect at certain temperatures and pressures more akin to surface conditions, at the conditions of a typical reservoir such as that used in this work (15 MPa and 90 ° C), hydrogen has a Z-Factor of 1.04 showing that it essentially behaves as an ideal gas.4 Ideal gases have Joule-Thomson coefficients of zero under all conditions,112 meaning that ideal gases experience no change in temperature when undergoing an adiabatic expansion, regardless of the pressure or initial temperature. Therefore, the Joule-Thomson effect would not be expected to have any substantial influence on the fault slip behavior. Furthermore, this effect would be expected to be more substantial if the fault was hydraulically connected to the injection interval. Given the fault is in the basement unit, this thermally derived pore fluid pressure change should not disturb the fault. In some subsurface geo-energy activities, phase changes can result in significant increases in volume, which could lead to rapid rises in pore fluid pressure which would influence the fault slip behavior. However, this is not the case for hydrogen as hydrogen remains in the gaseous state at all subsurface temperatures and pressures. Injecting the fluid at a temperature significantly different from the reservoir temperature may induce temperature changes and thermal stresses depending on the thermal conductivity and heat capacity of the reservoir that can influence fault apertures and thus fault permeability and pressurization.113 However, these will occur distal from the basement fault. Ideally, the hydrogen should be injected at a temperature similar to the reservoir temperature as in this work to avoid potential additional thermal stresses.

The present study employs a static friction model based on the Coulomb criterion to characterize fault behavior. The model delineates two distinct states according to Equation 9: stick and slip. In the stick state, the frictional traction is constrained to be less than or equal to the combined effect of cohesion and the product of the friction coefficient and normal traction. Conversely, in the slip state, the frictional traction exceeds this sum. While some might perceive this as a modeling constraint, adopting the rate-and-state model114,115,116 introduces additional uncertainties. Notably, the parameters a and b required for the rate-and-state model exhibit significant variability across faults and geological formations, posing challenges for model calibration. Moreover, the availability of a and b values for the specific rock units studied here is limited. However, the present static friction model partly addresses the “evolution effect”117 inherent in the rate-and-state model, as it considers changes in fault aperture over time due to geomechanical-fluid interactions. As fault apertures evolve over time, both tangential and normal tractions undergo changes, which dictate the occurrence of frictional stick-and-slip states on the fault. Consequently, the frictional resistance is subject to fluctuations due to mechanical and hydraulic variations over time. Since the model does not account explicitly for phenomena such as fault gouge compaction118 and cementation of fault rocks119,120 and the sealing and closure of the fault through mineral precipitation,121 pressure solution,122,123 and stress-induced cementation,124 then full fault aging/healing is not accounted for which would be complex to model accurately numerically. Despite these limitations, the model still allows for speculation on post-slip fault behavior if a rate-and-state model was employed. Additional simulations were run that examined the influence of the static friction coefficient between values of 0.50 and 1.00 on the fault slip. It was found that as the static friction coefficient increased, the fault slip decreased. For values between 0.50 and 0.65, the fault slip was nearly identical, but drastically reduced as the friction coefficient increased from 0.65 to 0.99. If the fault had an a - b value greater than zero, then the fault would be fault-strengthening post-slip and if it had an a - b value less than zero, then the fault would be fault-weakening post-slip. If the fault was fault-strengthening, this would have a similar effect to a larger friction coefficient which from our results would suggest a reduction in the predicted fault slip. If the fault was fault-weakening, then this would have a similar effect to a smaller friction coefficient which from our results would suggest an increase in the predicted fault slip. It is worth noting that rate-and-state effects typically induce minimal changes in the effective friction coefficient (i.e., +/− 0.05). Therefore, the slip is not expected to differ significantly from our results herein. This assertion is supported by our observation from our additional simulations that slip profiles for different friction coefficients exhibit similarity within the range of 0.50–0.70 which most faults in nature lie within, suggesting that temporal changes induced by a rate-and-state model would likely be insignificant.

The fault is approximated as a smooth macrofracture, but in reality faults within the Mesozoic shales of the North Sea can be relatively mature and display complex internal and external structures. These complexities include fault gouge,125 roughness,126 damage zones,127 and induced localized fracture networks. These features can influence fluid flow and stress distribution near the fault plane, affecting fault slip behavior. Ignoring such internal heterogeneity in a smooth macrofracture representation may lead to small inaccuracies in estimating fault slip. Fluid flow along faults can induce the formation of secondary fractures and fracture networks. These features can influence fluid flow pathways, pressure distribution, and stress transfer near the fault. Smooth macrofracture representations neglect the presence of such induced fractures, leading to small inaccuracies in the timing and magnitude of fault slip.

Seismic slip occurs when movement along a fault line leads to an earthquake,128 resulting in the sudden release of accumulated stress. This release causes surrounding rocks to fracture, generating seismic waves that propagate through the Earth, causing the ground to shake. In contrast, aseismic fault slip involves movement along a fault line without significant seismic wave generation or earthquake occurrence. This movement can be slow and continuous or episodic, allowing for gradual stress release over time. As this study does not model wave propagation, distinguishing between seismic and aseismic slip with certainty is not possible. However, seismic slip rates can range from fractions of a millimeter to several meters per event, while aseismic creep rates typically range from fractions of a millimeter to several centimeters per year.128 This annual aseismic creep rate range corresponds to a maximum of approximately 0.58 mm per week, which is the time step used in this work. Observations indicate slip easily exceeding 3 mm within a single time step, which supports seismic activity. The moment magnitudes provided in this study represent the maximum potential moment magnitudes, assuming instantaneous occurrence of all slip within a given time step. They serve as constraints on the upper limits of potential magnitudes, rather than predictions of event magnitudes because the distribution of slip within a given time step is non-discernible.

A single-phase fluid model is used within the reservoir unit. The ICGT does not currently have capabilities for multiphase flow as it is primarily a mechanics framework. In order to assess the influence of the in situ reservoir fluid on fault slip, an additional simulation was performed in which the hydrogen was replaced by brine of density 1,040 kg/m3. In this sense, the two end-member situations have been examined: (1) where the reservoir fluid is entirely hydrogen and (2) where the reservoir fluid is entirely brine. The real situation of brine and hydrogen coexisiting within the reservoir system should therefore lie between these two end-member states. The case of brine in the reservoir showed that the fault slip was up to approximately 5% larger than in the case with hydrogen in the reservoir. This additional slip arises because brine has a higher density than gaseous hydrogen. As a result, when brine is injected into the reservoir, it exerts higher pressure on the surrounding rock matrix compared to hydrogen. This increased pressure gradient contributes to enhanced stress perturbations near the fault plane, promoting fault slip. The change in slip between the hydrogen and brine case is small because the underburden fault and the reservoir fluid are not in direct hydraulic communication. Therefore, changes in pressure distribution in the reservoir associated with the in situ fluid can only impact fault slip via stress and strain transfer through the stiff underburden unit (E = 40 GPa), which does not facilitate much poroelastic deformation anyway. If the fault was in direct hydraulic communication with the reservoir fluid, it would be expected that the type of in situ fluid has a greater influence on the slip dynamics.

Conclusions

In this study, a 3-D, linear-elastic fracture mechanics simulator, with monolithic hydromechanical coupling was implemented to perform a parametric study aimed at assessing the influence of underburden mechanical properties on fault slip and induced seismicity along a basement fault during hydrogen storage into a depleted gas reservoir. The mechanical properties investigated were the Young’s Modulus, the Poisson’s ratio, the Biot coefficient, and the intrinsic permeability. Gaseous hydrogen was injected for five months, stored for a further two months, and withdrawn for a further five months, with this three-stage cycle repeated three times. The reservoir is assumed to be initially saturated with hydrogen. It was found that the underburden Young’s modulus is the key parameter with the strongest control on basement fault slip: fault slip significantly decreased as Young’s modulus increased. However, the increasing shear modulus of the underburden as the Young’s modulus increases, outweighed the decreasing slip, such that higher Young’s moduli resulted in larger magnitude seismicity. The Biot coefficient and Poisson’s ratio of the underburden unit had no substantial effect on induced fault slip. The fault slip was insensitive to changes in underburden permeability in the range 7×10−24 m2 to 1×10−19 m2, but increased significantly at times, when the underburden permeability was as high as 1×10−16 m2. Although the magnitude of the fault slip was largely insensitive to changes in underburden Poisson’s ratio, the maximum potential magnitude of the induced seismicity was most sensitive to the Poisson’s ratio from the tested parameter space. The maximum potential magnitude of the induced seismicity was least sensitive to changes in Young’s modulus. The underburden rock with low Young’s modulus (<15 GPa), high Poisson’s ratio (>0.25), low Biot coefficient and low matrix permeability (<1×10−19 m2) is found to minimize the magnitude of induced seismicity during the simulation of hydrogen storage in a depleted gas reservoir.

STAR★Methods

Key resources table

REAGENT or RESOURCE	SOURCE	IDENTIFIER	
Software and algorithms	
	
Imperial College Geomechanics Toolkit (ICGT)	Dr Adriana Paluszny (Imperial College London)100,101	N/A	
Rhinoceros 3D	TLM, Inc. https://www.rhino3d.com	N/A	
MayaVi	MayaVi https://mayavi.sourceforge.net/download.html	N/A	
Paraview	Sandia National Laboratories, Kitware Inc, Los Alamos National Laboratory https://www.paraview.org/download/	N/A	

Resource availability

Lead contact

Further information and requests for resources should be directed to and will be fulfilled by the lead contact, James Burtonshaw (james.burtonshaw16@imperial.ac.uk).

Materials availability

This study did not generate any new materials.

Data and code availability

• All data reported in this paper will be shared by the lead contact upon request.

• This paper does not report original code.

• Any additional information required to reanalyze the data reported in this paper is available from the lead contact upon request.

Method details

Simulations have been conducted using the Imperial College Geomechanics Toolkit (ICGT),100,101 a 3-D, monolithically coupled, hydromechanical simulator that is based on linear elastic fracture mechanics (LEFM). It is a computational geomechanics framework, written in C++, and uses the Finite Element Method (FEM). The ICGT interacts with external libraries and software such as Rhino for geometry construction, ANSYS ICEM for octree-based remeshing, and the Intel Pardiso direct solver for the inversion of the accumulated FEM matrix.

Fault flow and mechanical deformation of the rock are coupled by hydraulically loading the fault faces and ensuring compatibility of fault volumetric strains. Matrix flow and mechanical deformation of the rock are coupled via effective stress. Finally, fault flow and matrix flow are coupled via a leak-off mass transfer term that accounts for the loss of fluid from the fault to the matrix.

Specifically with relation to this work, the code has been applied and validated in the context of mechanical fracture contact and friction,129 hydrodynamics of CO2 injection-induced changes in fault geometry,113 induced seismicity resultant from CO2 injection56 and H2 injection, storage and production,4 and poroelasticity.130

Burtonshaw et al.56 conducted a comprehensive validation study in which the fault slip solution as a function of time from the ICGT was compared against that from COMSOL Multiphysics for a published case of CO2 injection.106 The model consisted of a 5-unit layer-cake system consisting of an upper aquifer, an upper caprock, a reservoir, a lower caprock and a lower aquifer unit, with an 80 ° dipping fault that penetrates to some extent all five layers. The original model of Mortezaei and Vahedifard106 was 2-D and this was replicated and extended to 3-D in the ICGT in the work of Burtonshaw et al.56 Furthermore, the 2-D model in Mortezaei and Vahedifard106 was thermohydromechanical (THM) in nature, whereas the 3-D model in Burtonshaw et al.56 was hydromechanical in nature (HM) with the thermal aspect being accounted for by attributing different fluid properties in each unit according to the pressure and temperature at the mid-depth of each unit from Mortezaei and Vahedifard.106

In both studies, the permeability of the reservoir unit was varied in three separate simulations (10−12, 10−13 and 10−14 m2) and the porosity of the reservoir unit in three further simulations (0.10, 0.15, and 0.20). The effect on fault slip as a function of time was computed. It was found that both models yielded similar results.

These results further support the ICGT as a suitable computational tool for the computation of fault slip during the injection of a geo-energy fluid.

Stress and strain

The rock matrix is assumed to be isotropic, homogeneous, and linear-elastic. Stress and strain are computed according to:(Equation 2) σ=D(ε−ε0)+σ0

where σ is the Cauchy stress tensor, D is the linear elastic stiffness matrix, ε is the infinitesimal strain tensor, ε0 is the initial strain and σ0 is the initial stress.131

Cauchy’s first law of motion must also be obeyed:(Equation 3) ∇·σ+F=0

where F is the external body forces acting on the domain, per unit volume.

Matrix flow, fault flow and deformation

Fault flow assumes a laminar flow model based on Lubrication theory (Equation 4).132,133 Matrix flow is modeled by the combination of Darcy’s law and mass conservation of the fluid (Equation 5). Fault flow, matrix flow and poroelastic deformation of the rock are concomitantly accounted for by the ICGT. Equation 3 is transformed into Equation 6 which is the final governing equation, by combining effective stress and imposing a boundary traction of hydraulic loading on fault faces, and subsequently integrating over each element.130(Equation 4) ∇·(af312μf∇pf)=afcf∂pf∂t+∂af∂t−kmμf∂pm∂nc

(Equation 5) ∫Ω∇·(kmμf∇pm+ρg)dΩ=∫Γckmμf∂pf∂ncdΓ+∫Ω[α∂(∇·u)∂t+(ϕcf+α−ϕKs)∂pm∂t]dΩ

(Equation 6) ∫Ω[∇·(Dε−αpmI)+F]dΩ=∫ΓcpfncdΓ

The nomenclature for Equations 4 to 6 is as follows: af = fault aperture, μf = fluid viscosity, pf = fracture pressure, ρ = fluid density, cf = fluid compressibility, km = matrix permeability, pm = matrix fluid pressure, nc = outward unit normal from fault faces, α = Biot coefficient, I = second-order identity matrix, Ω = domain, Γ = domain boundary, Γc = fault boundary, ϕ = porosity, Ks = matrix bulk modulus and u = rock matrix displacement. Note that in Equation 6, it is assumed that initial stress and initial strain are already in equilibrium.

Friction, contact and fault slip

The fault is a discrete discontinuity in the rock mass, represented by two smooth surfaces that act as the fault faces. When the two surfaces contact one another, a frictional contact constitutive law governs the boundary conditions on the fault surface. The law is derived from the Amontons-Coulomb first law of friction.

A normal gap function, gN, detects whether contact between fault surfaces has occurred and enforces displacement constraints. It describes the perpendicular gap between nodes on the two opposing fault surfaces (the follower, f, and the leader, l, surfaces):(Equation 7) gN=(xf−xl)·n=(uf−ul)·n−gˆN

where x is a position vector, n is the unit normal vector, u is a displacement vector and gˆN is a scalar defining the normal gap prior to initial displacement.

Once first contact of the fault faces occurs, an initial tangential gap vector, gˆT is defined, which measures the tangential displacement between two matching nodes on the follower and leader surfaces before the point of contact, and is given by gˆT = (I−n⊗n)(uf−ul)c. The tangential gap function, gT, measuring the same phenomenon after the point of initial contact, and referred herein as fault slip, xn, is defined as:(Equation 8) gT=(I−n⊗n)(xf−xl)=(I−n⊗n)(uf−ul)−gˆT.

The frictional (tangential) friction, τ, is defined based on one of two potential friction states (stick or slip):(Equation 9) {gT=0and∂gT∂t=0ifτ≤μ|p|+τc(stick)τ=(μ|p|+τc)∂gT∂t|∂gT∂t|ifτ>μ|p|+τc(slip)

where p is the normal traction acting on the fault surfaces, τc is the fault cohesion and μ is the friction coefficient.

The normal and tangential Kuhn-Tucker inequality contact constraints, which are the boundary conditions on the fault, are written as follows:(Equation 10) gN≥0,p≤0,pgN=0onΓc×[0,T],

(Equation 11) γ˙≥0,|τ|−(μ|p|+τc)≤0,γ˙[|τ|−(μ|p|+τc)]=0onΓc×[0,T],

where γ˙ is the slip rate.

The problem is then solved using a four-stage gap-based double-loop Augmented-Lagrangian algorithm to fully resolve the aperture and traction distribution of the fault where gaps are augmented instead of tractions. The algorithm is described in more detail in Nejati et al.129,134

Seismic moment release and moment magnitude

The scalar seismic moment, M0, is given by:(Equation 12) M0=E2(1+ν)×A×d,

where E is the Young’s modulus, ν is the Poisson’s ratio, A is the rupture area of the fault and d is the average slip across the fault surface. The rupture area is computed as the areal region of the fault surface that experiences significant fault slip above 0.5 mm. For slip below this, the fault is said to be in stick, and the area does not contribute to the A term.

The moment magnitude, Mw, of an event, with M0 specified in units of Nm, is then computed as135,136:(Equation 13) Mw=23log10M0−9.1.

Quantification and statistical analysis

Visualization of the aperture, stress and strain results used the MayaVi and Paraview softwares. VTK’s of these parameters were internally created from running the Imperial College Geomechanics Toolkit, and subsequently visualized in the aforementioned external softwares. The fault slip at each node of the fault at each time step was outputted into a.txt file. This.txt file was read into a python script which applied an area-based weighted average of the fault slip at each node to obtain an averaged fault slip value at each time step. The exact formulation of this area-based weighted average is described in the STAR Methods section of this paper.

Acknowledgments

The authors thank the UK Natural Environment Research Council (NERC) for funding SeisGreen Project (grant no. NE/W009293/1 ) which supported this work. The authors also thank the 10.13039/501100000288 Royal Society UK for supporting this research, through fellowship UF160443 , awarded to Adriana Paluszny.

Author contributions

Conceptualization, J.E.J.B. and A.P.; methodology, J.E.J.B. and A.P.; software, A.P., A.M., and J.E.J.B.; validation, J.E.J.B., A.P., and A.M.; formal analysis, J.E.J.B. and A.P.; investigation, J.E.J.B.; resources, A.P. and R.W.Z.; data curation, A.P.; writing – original draft, J.E.J.B. and A.P.; writing – review & editing, J.E.J.B., A.P., R.W.Z., A.M.; visualization, J.E.J.B. and A.P.; supervision, A.P. and R.W.Z.; project administration, A.P. and R.W.Z.; funding acquisition, A.P. and R.W.Z.

Declaration of interests

The authors declare no competing interests.
==== Refs
References

1 Qyyum M.A. Dickson R. Ali Shah S.F. Niaz H. Khan A. Liu J.J. Lee M. Availability, versatility, and viability of feedstocks for hydrogen production: Product space perspective Renew. Sustain. Energy Rev. 145 2021 110843 10.1016/j.rser.2021.110843
2 Seck G.S. Hache E. Sabathier J. Guedes F. Reigstad G.A. Straus J. Wolfgang O. Ouassou J.A. Askeland M. Hjorth I. Hydrogen and the decarbonization of the energy system in Europe in 2050: A detailed model-based analysis Renew. Sustain. Energy Rev. 167 2022 112779 10.1016/j.rser.2022.112779
3 Griffiths S. Sovacool B.K. Kim J. Bazilian M. Uratani J.M. Industrial decarbonization via hydrogen: A critical and systematic review of developments, socio-technical systems and policy options Energy Res. Social Sci. 80 2021 102208 10.1016/j.erss.2021.102208
4 Burtonshaw J.E.J. Paluszny A. Mohammadpour A. Zimmerman R.W. Effects of Reservoir Mechanical Properties on Induced Seismicity during Subsurface Hydrogen Storage Phil. Trans. 382 2024 20230187 10.1098/rsta.2023.0187
5 Hassanpouryouzband A. Joonaki E. Edlmann K. Haszeldine R.S. Offshore geological storage of hydrogen: is this our best option to achieve net-zero? ACS Energy Lett. 6 2021 2181 2186 10.1021/acsenergylett.1c00845
6 Le Duigou A. Bader A.-G. Lanoix J.-C. Nadau L. Relevance and costs of large scale underground hydrogen storage in France Int. J. Hydrogen Energy 42 2017 22987 23003 10.1016/j.ijhydene.2017.06.239
7 Mayyas A. Wei M. Levis G. Hydrogen as a long-term, large-scale energy storage solution when coupled with renewable energy sources or grids with dynamic electricity pricing schemes Int. J. Hydrogen Energy 45 2020 16311 16325 10.1016/j.ijhydene.2020.04.163
8 Benavides K. Exploring the Role of Hydrogen in Decarbonizing Heavy Industry PhD thesis 2023 Massachusetts Institute of Technology https://hdl.handle.net/1721.1/151820
9 Neuwirth M. Fleiter T. Manz P. Hofmann R. The future potential hydrogen demand in energy-intensive industries-a site-specific approach applied to Germany Energy Convers. Manag. 252 2022 115052 10.1016/j.enconman.2021.115052
10 Liu W. Zuo H. Wang J. Xue Q. Ren B. Yang F. The production and application of hydrogen in steel industry Int. J. Hydrogen Energy 46 2021 10548 10569 10.1016/j.ijhydene.2020.12.123
11 Aminudin M. Kamarudin S. Lim B. Majilan E. Masdar M. Shaari N. An overview: Current progress on hydrogen fuel cell vehicles Int. J. Hydrogen Energy 48 2023 4371 4388 10.1016/j.ijhydene.2022.10.156
12 Heinemann N. Alcalde J. Miocic J.M. Hangx S.J.T. Kallmeyer J. Ostertag-Henning C. Hassanpouryouzband A. Thaysen E.M. Strobel G.J. Schmidt-Hattenberger C. Enabling large-scale hydrogen storage in porous media–the scientific challenges Energy Environ. Sci. 14 2021 853 864 10.1039/D0EE03536J
13 Shadidi B. Najafi G. Yusaf T. A review of hydrogen as a fuel in internal combustion engines Energies 14 2021 6209 10.3390/en14196209
14 Singla M.K. Nijhawan P. Oberoi A.S. Hydrogen fuel and fuel cell technology for cleaner future: a review Environ. Sci. Pollut. Res. Int. 28 2021 15607 15626 10.1007/s11356-020-12231-8 33538968
15 Fan Z. Ochu E. Braverman S. Lou Y. Smith G. Bhardwaj A. Brouwer J. McCormick C. Friedmann J. Green Hydrogen in a Circular Carbon Economy: Opportunities and Limits https://www.energypolicy.columbia.edu/publications/green-hydrogen-circular-carbon-economy-opportunities-and-limits/ 2021
16 WoodMackenzie The Future for Green Hydrogen https://www.woodmac.com/news/editorial/the-future-for-green-hydrogen/.WoodMackenzie 2019
17 Birol F. The Future of Hydrogen (IEA) 2019 International Energy Agency https://www.iea.org/reports/the-future-of-hydrogen
18 Zapantis A. Blue Hydrogen 2021 Global CCS Institute https://www.ourenergypolicy.org/resources/blue-hydrogen/
19 Powell D. Focus on Blue Hydrogen 2020 Gaffney Cline https://www.gaffneycline.com/sites/g/files/cozyhq681/files/2021-08/Focus_on_Blue_Hydrogen_Aug2020.pdf
20 IEA Hydrogen Production Costs by Production Source, 2018 https://www.iea.org/data-and-statistics/charts/hydrogen-production-costs-by-production-source-2018/ 2018
21 Prime J. Martínez L.M. Statistics Report: Coal Information - Overview Tech. rep. 2020 International Energy Agency (IEA) https://iea.blob.core.windows.net/assets/a5f208e9-f66b-4d31-b5af-87d581b70c18/Coal_Information_Overview_2020_edition.pdf
22 van Renssen S. The hydrogen solution? Nat. Clim. Chang. 10 2020 799 801 10.1038/s41558-020-0891-0
23 Flowers S. The Energy Transition: Future Energy from The Edge 2020 Tech. Rep. https://www.woodmac.com/news/the-edge/future-energy-from-the-edge/.WoodMackenzie
24 Howarth R.W. Jacobson M.Z. How green is blue hydrogen? Energy Sci. Eng. 9 2021 1676 1687 10.1002/ese3.956
25 Pinsky R. Sabharwall P. Hartvigsen J. O’Brien J. Comparative review of hydrogen production technologies for nuclear hybrid energy systems Prog. Nucl. Energy 123 2020 103317 10.1016/j.pnucene.2020.103317
26 Gielen D. Taibi E. Miranda R. Hydrogen: A renewable energy perspective Tech. rep. 2019 International Renewable Energy Agency (IRENA) https://www.irena.org/publications/2019/Sep/Hydrogen-A-renewable-energy-perspective
27 Timmerberg S. Kaltschmitt M. Finkbeiner M. Hydrogen and hydrogen-derived fuels through methane decomposition of natural gas–GHG emissions and costs Energy Convers. Manag. X 7 2020 100043 10.1016/j.ecmx.2020.100043
28 Coutanceau C. Baranton S. Audichon T. Hydrogen Electrochemical Production 2017 Academic Press
29 Chen L. Qi Z. Zhang S. Su J. Somorjai G.A. Catalytic hydrogen production from methane: A review on recent progress and prospect Catalysts 10 2020 858 10.3390/catal10080858
30 Meloni E. Martino M. Palma V. A short review on Ni based catalysts and related engineering issues for methane steam reforming Catalysts 10 2020 352 10.3390/catal10030352
31 Truche L. Bazarkina E.F. Natural hydrogen the fuel of the 21st century E3S Web Conf. 98 2019 EDP Sciences 03006 10.1051/e3sconf/20199803006
32 Andreani M. Daniel I. Pollet-Villard M. Aluminum speeds up the hydrothermal alteration of olivine Am. Mineral. 98 2013 1738 1744 10.2138/am.2013.4469
33 Harman J. Hjalmarsson P. Mermelstein J. Ryley J. Sadler H. Selby M. 1MW-Class Solid Oxide Electrolyser System Prototype for Low-Cost Green Hydrogen ECS Trans. 103 2021 383 392 10.1149/10301.0383ecst
34 Xu D. Dong L. Ren J. Chapter 2 - Introduction of Hydrogen Routines Scipioni A. Manzardo A. Ren J. Hydrogen Economy 2017 Academic Press 35 54 10.1016/B978-0-12-811132-1.00002-X
35 Newborough M. Cooley G. Developments in the global hydrogen market: The spectrum of hydrogen colours Fuel Cell. Bull. 2020 2020 16 22 10.1016/S1464-2859(20)30546-0
36 Leal Pérez B.J. Medrano Jiménez J.A. Bhardwaj R. Goetheer E. van Sint Annaland M. Gallucci F. Methane pyrolysis in a molten gallium bubble column reactor for sustainable hydrogen production: Proof of concept & techno-economic assessment Int. J. Hydrogen Energy 46 2021 4917 4935 10.1016/j.ijhydene.2020.11.079
37 Osselin F. Soulaine C. Fauguerolles C. Gaucher E.C. Scaillet B. Pichavant M. Orange hydrogen is the new green Nat. Geosci. 15 2022 765 769 10.1038/s41561-022-01043-9
38 Hassanpouryouzband A. Wilkinson M. Haszeldine R.S. Hydrogen energy futures–foraging or farming? Chem. Soc. Rev. 53 2024 2258 2263 10.1039/D3CS00723E 38323342
39 Moradpoor I. Syri S. Santasalo-Aarnio A. Green hydrogen production for oil refining–Finnish case Renew. Sustain. Energy Rev. 175 2023 113159 10.1016/j.rser.2023.113159
40 Saygin D. Blanco H. Boshell F. Cordonnier J. Rouwenhorst K. Lathwal P. Gielen D. Ammonia production from clean hydrogen and the implications for global natural gas demand Sustainability 15 2023 1623 10.3390/su15021623
41 Ledari M.B. Khajehpour H. Akbarnavasi H. Edalati S. Greening steel industry by hydrogen: Lessons learned for the developing world Int. J. Hydrogen Energy 48 2023 36623 36649 10.1016/j.ijhydene.2023.06.058
42 Yang Z. Li H. Zhang H. Dynamic Collaborative Pricing for Managing Refueling Demand of Hydrogen Fuel Cell Vehicles IEEE Trans. Transp. Electrific. 1 2024 10.1109/TTE.2024.3381236
43 Thomas G. Pidgeon N. Henwood K. Hydrogen, a less disruptive pathway for domestic heat? Exploratory findings from public perceptions research Clean. Prod. Lett. 5 2023 100047 10.1016/j.clpl.2023.100047
44 Groll M. Can climate change be avoided? Vision of a hydrogen-electricity energy economy Energy 264 2023 126029 10.1016/j.energy.2022.126029
45 Jesse B.-J. Kramer G.J. Koning V. Vögele S. Kuckshinrichs W. Stakeholder perspectives on the scale-up of green hydrogen and electrolyzers Energy Rep. 11 2024 208 217 10.1016/j.egyr.2023.11.046
46 Genovese M. Fragiacomo P. Hydrogen refueling station: overview of the technological status and research enhancement J. Energy Storage 61 2023 106758 10.1016/j.est.2023.106758
47 Morales-España G. Hernández-Serna R. Tejada-Arango D.A. Weeda M. Impact of large-scale hydrogen electrification and retrofitting of natural gas infrastructure on the European power system Int. J. Electr. Power Energy Syst. 155 2024 109686 10.1016/j.ijepes.2023.109686
48 Butler M. Underhill J.R. A Critical Geological Evaluation of the Hydrogen Storage Potential in the Cousland Gas Field, Midland Valley of Scotland Earth Sci. Syst. Soc. 3 2024 10076 10.3389/esss.2023.10076
49 Corina A. Zikovic V. Soustelle V. Moghadam A. Groenenberg R. Review of current practices and existing experimental data for well and rock materials under cyclic hydrogen injection and withdrawal H2020 HyUSPRe project report D 5 https://www.hyuspre.eu/wp-content/uploads/2022/06/HyUSPRe_D5.1_Review-of-current-practices-and-existing-experimental-data-for-well-and-rock-materials-under-cyclic-hydrogen-injection-and-withdrawal_2022.05.31.pdf 2022
50 Krevor S. De Coninck H. Gasda S.E. Ghaleigh N.S. de Gooyert V. Hajibeygi H. Juanes R. Neufeld J. Roberts J.J. Swennenhuis F. Subsurface carbon dioxide and hydrogen storage for a sustainable energy future Nat. Rev. Earth Environ. 4 2023 102 118 10.1038/s43017-022-00376-8
51 Bai T. Tahmasebi P. Coupled hydro-mechanical analysis of seasonal underground hydrogen storage in a saline aquifer J. Energy Storage 50 2022 104308 10.1016/j.est.2022.104308
52 Zeng L. Sarmadivaleh M. Saeedi A. Chen Y. Zhong Z. Xie Q. Storage integrity during underground hydrogen storage in depleted gas reservoirs Earth Sci. Rev. 247 2023 104625 10.1016/j.earscirev.2023.104625
53 Perera M. A review of underground hydrogen storage in depleted gas reservoirs: Insights into various rock-fluid interaction mechanisms and their impact on the process integrity Fuel 334 2023 126677 10.1016/j.fuel.2022.126677
54 AlDhuhoori M. Belhaj H. AlHameli F. A Mathematical Model for Formation Caprock Integrity Incorporating Creep Deformation Mechanism: A Hydrogen Storage Seasonal Case Study Abu Dhabi International Petroleum Exhibition and Conference 2023 SPE 10.2118/216991-MS D041S141R006
55 Kumar K.R. Honorio H.T. Hajibeygi H. Simulation of the inelastic deformation of porous reservoirs under cyclic loading relevant for underground hydrogen storage Sci. Rep. 12 2022 21404 10.1038/s41598-022-25715-z
56 Burtonshaw J.E. Paluszny A. Zimmerman R.W. Numerical Modelling of Induced Seismicity Along a Fault During CO2 Injection into a Subsurface Reservoir Proceedings of the 15th ISRM Congress 2023 ISRM PaperISRM-2023–1688 https://onepetro.org/isrmcongress/proceedings-abstract/CONGRESS23/All-CONGRESS23/539787
57 Francke H. Thorade M. Density and viscosity of brine: An overview from a process engineers perspective Geochemistry 70 2010 23 32 10.1016/j.chemer.2010.05.015
58 Burtonshaw J. Thomas R. Paluszny A. Zimmerman R. Effect of Viscosity and Injection Rate on the Tensile and Shear Deformation Experienced by Interacting Hydraulic and Natural Fractures ARMA US Rock Mechanics/Geomechanics Symposium 2021 ARMA,ARMA–2021 https://onepetro.org/ARMAUSRMS/proceedings-abstract/ARMA21/All-ARMA21/468161
59 Ramesh Kumar K. Honorio H. Chandra D. Lesueur M. Hajibeygi H. Comprehensive review of geomechanics of underground hydrogen storage in depleted reservoirs and salt caverns J. Energy Storage 73 2023 108912 10.1016/j.est.2023.108912
60 Kettlety T. Verdon J.P. Fault triggering mechanisms for hydraulic fracturing-induced seismicity from the Preston New Road, UK Case Study Front. Earth Sci. 9 2021 670771 10.3389/feart.2021.670771
61 Evans D.M. The Denver area earthquakes and the Rocky Mountain Arsenal disposal well Mt. Geol. 1966 25 32 10.1130/Eng-Case-8.25
62 Whitney J. Resistance in the Midst of the Anthropocene: The Rise and Fall of Artificial Earthquakes at the Rocky Mountain Arsenal, 1962–1966 2019 Arcadia 10.5282/rcc/8880
63 Kuchment A. Drilling for earthquakes Sci. Am. 315 2016 46 53 10.1038/scientificamerican0716-46
64 Wang Z. Krupnick A. A retrospective review of shale gas development in the United States: What led to the boom? Econ. Energy Environ. Policy 4 2015 5 18 10.5547/2160-5890.4.1.zwan
65 Frohlich C. Two-year survey comparing earthquake activity and injection-well locations in the Barnett Shale, Texas Proc. Natl. Acad. Sci. USA 109 2012 13934 13938 10.1073/pnas.1207728109 22869701
66 Frohlich C. Hayward C. Stump B. Potter E. The Dallas–Fort Worth earthquake sequence: October 2008 through May 2009 Bull. Seismol. Soc. Am. 101 2011 327 340 10.1785/0120100131
67 Park Y. Mousavi S.M. Zhu W. Ellsworth W.L. Beroza G.C. Machine-learning-based analysis of the Guy-Greenbrier, Arkansas earthquakes: A tale of two sequences Geophys. Res. Lett. 47 2020 e2020GL087032 10.1029/2020GL087032
68 Kim W.-Y. Induced seismicity associated with fluid injection into a deep well in Youngstown, Ohio JGR. Solid Earth 118 2013 3506 3518 10.1002/jgrb.50247
69 Cao W. Shi J.-Q. Durucan S. Korre A. Evaluation of shear slip stress transfer mechanism for induced microseismicity at in Salah CO2 storage site Int. J. Greenh. Gas Control 107 2021 103302 10.1016/j.ijggc.2021.103302
70 Yeck W.L. Hayes G.P. McNamara D.E. Rubinstein J.L. Barnhart W.D. Earle P.S. Benz H.M. Oklahoma experiences largest earthquake during ongoing regional wastewater injection hazard mitigation efforts Geophys. Res. Lett. 44 2017 711 717 10.1002/2016GL071685
71 Boak J. Patterns of induced seismicity in central and northwest Oklahoma SEG Technical Program Expanded Abstracts 2016 2016 Society of Exploration Geophysicists 5039 5042 10.1190/segam2016-13960154.1
72 Meng H. Ge H. Fu D. Wang X. Shen Y. Jiang Z. Wang J. Numerical investigation of casing shear deformation due to fracture/fault slip during hydraulic fracturing Energy Sci. Eng. 8 2020 3588 3601 10.1002/ese3.766
73 Hui G. Chen S. Gu F. Pang Y. Yu X. Zhang L. Insights on controlling factors of hydraulically induced seismicity in the Duvernay East Shale Basin Geochem. Geophys. Geosyst. 22 2021 e2020GC009563 10.1029/2020GC009563
74 Fan Z. Eichhubl P. Newell P. Basement fault reactivation by fluid injection into sedimentary reservoirs: Poroelastic effects JGR. Solid Earth 124 2019 7354 7369 10.1029/2018JB017062
75 Hemami B. Feizi Masouleh S. Ghassemi A. Basement fault reactivation in response to injection into overlying layers ARMA US Rock Mechanics/Geomechanics Symposium 2021 ARMA https://onepetro.org/ARMAUSRMS/proceedings-abstract/ARMA21/All-ARMA21/468085.ARMA
76 Zoback M. Smit D. Meeting the challenges of large-scale carbon storage and hydrogen production Proc. Natl. Acad. Sci. USA 120 2023 e2202397120 10.1073/pnas.2202397120
77 Hearn E.H. Koltermann C. Rubinstein J.L. Numerical models of pore pressure and stress changes along basement faults due to wastewater injection: Applications to the 2014 Milan, Kansas earthquake Geochem. Geophys. Geosyst. 19 2018 1178 1198 10.1002/2017GC007194
78 Horton S. Disposal of hydrofracking waste fluid by injection into subsurface aquifers triggers earthquake swarm in central Arkansas with potential for damaging earthquake Seismol Res. Lett. 83 2012 250 260 10.1785/gssrl.83.2.250
79 Zhang Y. Person M. Rupp J. Ellett K. Celia M.A. Gable C.W. Bowen B. Evans J. Bandilla K. Mozley P. Hydrogeologic controls on induced seismicity in crystalline basement rocks due to fluid injection into basal reservoirs Groundwater 51 2013 525 538 10.1111/gwat.12071
80 Verdon J.P. Significance for secure CO2 storage of earthquakes induced by fluid injection Environ. Res. Lett. 9 2014 064022 10.1088/1748-9326/9/6/064022
81 Hornbach M.J. DeShon H.R. Ellsworth W.L. Stump B.W. Hayward C. Frohlich C. Oldham H.R. Olson J.E. Magnani M.B. Brokaw C. Luetgert J.H. Causal factors for seismicity near Azle, Texas Nat. Commun. 6 2015 6728 6811 10.1038/ncomms7728 25898170
82 McGarr A. Bekins B. Burkardt N. Dewey J. Earle P. Ellsworth W. Ge S. Hickman S. Holland A. Majer E. Coping with earthquakes induced by fluid injection Science 347 2015 830 831 10.1126/science.aaa0494 25700505
83 Chang K.W. Yoon H. Martinez M.J. Potential seismicity along basement faults induced by geological carbon sequestration Geophys. Res. Lett. 49 2022 e2022GL098721 10.1029/2022GL098721
84 Buijze L. van Bijsterveldt L. Cremer H. Jaarsma B. Paap B. Veldkamp H. Wassing B. van Wees J. van Yperen G. Ter Heege J. Induced seismicity in geothermal systems: Occurrences worldwide and implications for the Netherlands European Geothermal Congress 2019 11 14 https://www.researchgate.net/profile/Loes-Buijze/publication/334520139_Induced_seismicity_in_geothermal_systems_Occurrences_worldwide_and_implications_for_the_Netherlands/links/5d2f2902458515c11c37d304/Induced-seismicity-in-geothermal-systems-Occurrences-worldwide-and-implications-for-the-Netherlands.pdf
85 Megies T. Wassermann J. Microseismicity observed at a non-pressure-stimulated geothermal power plant Geothermics 52 2014 36 49 10.1016/j.geothermics.2014.01.002
86 Küperkoch L. Olbert K. Meier T. Long-term monitoring of induced seismicity at the Insheim geothermal site, Germany Bull. Seismol. Soc. Am. 108 2018 3668 3683 10.1785/0120170365
87 Ellsworth W.L. Injection-induced earthquakes Science 341 2013 1225942 10.1126/science.1225942
88 Wang L. Kwiatek G. Bohnhoff M. Rybacki E. Dresen G. Injection-induced fault slip and associated seismicity in the lab: Insights from source mechanisms, local stress states and fault geometry Earth Planet Sci. Lett. 626 2024 118515 10.1016/j.epsl.2023.118515
89 Scuderi M. Collettini C. Marone C. Frictional stability and earthquake triggering during fluid pressure stimulation of an experimental fault Earth Planet Sci. Lett. 477 2017 84 96 10.1016/j.epsl.2017.08.009
90 Bao X. Eaton D.W. Fault activation by hydraulic fracturing in western Canada Science 354 2016 1406 1409 10.1126/science.aag2583 27856850
91 King Hubbert M. Rubey W.W. Role of fluid pressure in mechanics of overthrust faulting: I. Mechanics of fluid-filled porous solids and its application to overthrust faulting Geol. Soc. Am. Bull. 70 1959 115 166 10.1130/0016-7606(1959)70[115:ROFPIM]2.0.CO;2
92 Shapiro S. Patzig R. Rothert E. Rindschwentner J. Triggering of seismicity by pore-pressure perturbations: Permeability-related signatures of the phenomenon Thermo-Hydro-Mechanical Coupling in Fractured Rock 2003 1051 1066
93 McGarr A. Maximum magnitude earthquakes induced by fluid injection J. Geophys. Res. Solid Earth 119 2014 1008 1019 10.1002/2013JB010597
94 Barbour A.J. Norbeck J.H. Rubinstein J.L. The effects of varying injection rates in Osage County, Oklahoma, on the 2016 M w 5.8 Pawnee earthquake Seismol Res. Lett. 88 2017 1040 1053 10.1785/0220170003
95 Segall P. Lu S. Injection-induced seismicity: Poroelastic and earthquake nucleation effects JGR. Solid Earth 120 2015 5082 5103 10.1002/2015JB012060
96 Chang K.W. Segall P. Injection-induced seismicity on basement faults including poroelastic stressing JGR. Solid Earth 121 2016 2708 2726 10.1002/2015JB012561
97 Chang K.W. Segall P. Seismicity on basement faults induced by simultaneous fluid injection–extraction Pure Appl. Geophys. 173 2016 2621 2636 10.1007/s00024-016-1319-7
98 Ogwari P.O. DeShon H.R. Hornbach M.J. The Dallas-Fort Worth airport earthquake sequence: Seismicity beyond injection period JGR. Solid Earth 123 2018 553 563 10.1002/2017JB015003
99 Langenbruch C. Weingarten M. Zoback M.D. Physics-based forecasting of man-made earthquake hazards in Oklahoma and Kansas Nat. Commun. 9 2018 3946 10.1038/s41467-018-06167-4 30258058
100 Paluszny A. Zimmerman R.W. Numerical fracture growth modeling using smooth surface geometric deformation Eng. Fract. Mech. 108 2013 19 36 10.1016/j.engfracmech.2013.04.012
101 Thomas R.N. Paluszny A. Zimmerman R.W. Growth of three-dimensional fractures, arrays, and networks in brittle rocks under tension and compression Comput. Geotech. 121 2020 103447 10.1016/j.compgeo.2020.103447
102 Alghannam M. Juanes R. Understanding rate effects in injection-induced earthquakes Nat. Commun. 11 2020 3053 10.1038/s41467-020-16860-y 32546793
103 Tang L. Lu Z. Zhang M. Sun L. Wen L. Seismicity induced by simultaneous abrupt changes of injection rate and well pressure in Hutubi gas field JGR. Solid Earth 123 2018 5929 5944 10.1029/2018JB015863
104 Cuenot N. Dorbath C. Dorbath L. Analysis of the microseismicity induced by fluid injections at the EGS site of Soultz-sous-Forêts (Alsace, France): implications for the characterization of the geothermal reservoir properties Pure Appl. Geophys. 165 2008 797 828 10.1007/s00024-008-0335-7
105 Ayodele T.R. Mosetlhe T.C. Yusuff A.A. Ogunjuyigbe A. Off-grid hybrid renewable energy system with hydrogen storage for South African rural community health clinic Int. J. Hydrogen Energy 46 2021 19871 19885 10.1016/j.ijhydene.2021.03.140
106 Mortezaei K. Vahedifard F. Numerical simulation of induced seismicity in carbon capture and storage projects Geotech. Geol. Eng. 33 2015 411 424 10.1007/s10706-015-9859-7
107 Deng Q. Blöcher G. Cacace M. Schmittbuhl J. Modeling of fluid-induced seismicity during injection and after shut-in Comput. Geotech. 140 2021 104489 10.1016/j.compgeo.2021.104489
108 Rutqvist J. Rinaldi A.P. Cappa F. Moridis G.J. Modeling of fault reactivation and induced seismicity during hydraulic fracturing of shale-gas reservoirs J. Petrol. Sci. Eng. 107 2013 31 44 10.1016/j.petrol.2013.04.023
109 Zbinden D. Rinaldi A.P. Urpi L. Wiemer S. On the physics-based processes behind production-induced seismicity in natural gas fields JGR. Solid Earth 122 2017 3792 3812 10.1002/2017JB014003
110 Mazzoldi A. Rinaldi A.P. Borgia A. Rutqvist J. Induced seismicity within geological carbon sequestration projects: Maximum earthquake magnitude and leakage potential from undetected faults Int. J. Greenh. Gas Control 10 2012 434 442 10.1016/j.ijggc.2012.07.012
111 Zhang Z. Gao M. Chen X. Wei X. Liang J. Wu C. Wang L. The Joule–Thomson effect of (CO2+ H2) binary system relevant to gas switching reforming with carbon capture and storage (CCS) Chin. J. Chem. Eng. 54 2023 215 231 10.1016/j.cjche.2022.03.017
112 Saygın H. Şişman A. Joule–Thomson coefficients of quantum ideal-gases Applied energy 70 2001 49 57 10.1016/S0306-2619(01)00018-6
113 Salimzadeh S. Paluszny A. Zimmerman R.W. Effect of cold CO2 injection on fracture apertures and growth Int. J. Greenh. Gas Control 74 2018 130 141 10.1016/j.ijggc.2018.04.013
114 Rice J.R. Lapusta N. Ranjith K. Rate and state dependent friction and the stability of sliding between elastically deformable solids J. Mech. Phys. Solid. 49 2001 1865 1898 10.1016/S0022-5096(01)00042-4
115 Hosseini N. Priest J.A. Eaton D.W. Extended-FEM analysis of injection-induced slip on a fault with rate-and-state friction: Insights into parameters that control induced seismicity Rock Mech. Rock Eng. 56 2023 4229 4250 10.1007/s00603-023-03283-6
116 Garagash D.I. Fracture mechanics of rate-and-state faults and fluid injection induced slip Philos. Trans. A Math. Phys. Eng. Sci. 379 2021 20200129 10.1098/rsta.2020.0129
117 Liang C. Ampuero J.-P. Pino Muñoz D. The paucity of supershear earthquakes on large faults governed by rate and state friction Geophys. Res. Lett. 49 2022 e2022GL099749 10.1029/2022GL099749
118 Proctor B. Lockner D.A. Kilgore B.D. Mitchell T.M. Beeler N.M. Direct evidence for fluid pressure, dilatancy, and compaction affecting slip in isolated faults Geophys. Res. Lett. 47 2020 e2019GL086767 10.1029/2019GL086767
119 Williams R.T. Fagereng Å. The role of quartz cementation in the seismic cycle: A critical review Rev. Geophys. 60 2022 e2021RG000768 10.1029/2021RG000768
120 Romano C. Williams R. Cheng F. Evolution of fault-zone hydromechanical properties in response to different cementation processes Lithosphere 2022 2022 1069843 10.2113/2022/1069843
121 Yehya A. Rice J.R. Influence of fluid-assisted healing on fault permeability structure JGR. Solid Earth 125 2020 e2020JB020553 10.1029/2020JB020553
122 Bedford J.D. Hirose T. Hamada Y. Rapid fault healing after seismic slip JGR. Solid Earth 128 2023 e2023JB026706 10.1029/2023JB026706
123 Hunfeld L.B. Chen J. Hol S. Niemeijer A.R. Spiers C.J. Healing behavior of simulated fault gouges from the Groningen gas field and implications for induced fault reactivation J. Geophys. Res. Solid Earth 125 2020 e2019JB018790 10.1029/2019JB018790
124 Ji Y. Hofmann H. Rutter E.H. Zang A. Transition From Slow to Fast Injection-Induced Slip of an Experimental Fault in Granite Promoted by Elevated Temperature Geophys. Res. Lett. 49 2022 e2022GL101212 10.1029/2022GL101212
125 Buijze L. Guo Y. Niemeijer A. Ma S. Spiers C. Effects of heterogeneous gouge segments on the slip behavior of experimental faults at dm scale Earth Planet Sci. Lett. 554 2021 116652 10.1016/j.epsl.2020.116652
126 Wang L. Kwiatek G. Renard F. Guérin-Marthe S. Rybacki E. Bohnhoff M. Naumann M. Dresen G. Fault roughness controls injection-induced seismicity Proc. Natl. Acad. Sci. USA 121 2024 e2310039121 10.1073/pnas.2310039121
127 Thakur P. Huang Y. Kaneko Y. Effects of low-velocity fault damage zones on long-term earthquake behaviors on mature strike-slip faults JGR. Solid Earth 125 2020 e2020JB019587 10.1029/2020JB019587
128 Wallace L.M. Slow slip events in New Zealand Annu. Rev. Earth Planet Sci. 48 2020 175 203 10.1146/annurev-earth-071719-055104
129 Nejati M. Paluszny A. Zimmerman R.W. A finite element framework for modeling internal frictional contact in three-dimensional fractured media using unstructured tetrahedral meshes Comput. Methods Appl. Mech. Eng. 306 2016 123 150 10.1016/j.cma.2016.03.028
130 Salimzadeh S. Paluszny A. Zimmerman R.W. Three-dimensional poroelastic effects during hydraulic fracturing in permeable rocks Int. J. Solid Struct. 108 2017 153 163 10.1016/j.ijsolstr.2016.12.008
131 Cook R.D. Concepts and Applications of Finite Element Analysis 2007 John wiley & sons
132 Lugt P.M. Morales-Espejel G.E. A review of elasto-hydrodynamic lubrication theory Tribology Trans. 54 2011 470 496 10.1080/10402004.2010.551804
133 Zimmerman R. Kumar S. Bodvarsson G. Lubrication theory analysis of the permeability of rough-walled fractures Int. J. Rock Mech. Min. Sci. Geomech. Abstracts 28 1991 325 331 10.1016/0148-9062(91)90597-F
134 Nejati M. Finite Element Modeling of Frictional Contact and Stress Intensity Factors in Three-Dimensional Fractured Media Using Unstructured Tetrahedral Meshes PhD thesis 2015 Imperial College London
135 Bormann P. Wendt S. DiGiacomo D. Seismic sources and source parameters New manual of seismological observatory practice 2 (NMSOP2) 2013 Deutsches GeoForschungsZentrum GFZ 1 259 10.2312/GFZ.NMSOP-2_ch3
136 Bormann P. Di Giacomo D. The moment magnitude M w and the energy magnitude M e: common roots and differences J. Seismol. 15 2011 411 427 10.1007/s10950-010-9219-2
