Influence of Rock Properties on Structural Failure Probability—Caprock Shale Examples from the Horda Platform, Offshore Norway

Abstract

In any geological subsurface fluid injection, a viable top seal is required to contain the vertical movement of the injected fluid plume. However, seal integrity assessment is challenging because of the uncertainties possessed by various parameters. A probabilistic solution might be more appropriate when significant uncertainties are present. In this study, we evaluate Drake shale caprock structural reliability using a stochastic method. Drake shale is the primary top seal in the Aurora CO₂ storage site, located in the Horda Platform area in the northern North Sea. Based on the lithological variations, Drake caprock shale is divided into two parts designated by upper and lower units. Six model scenarios from the upper and lower Drake units have been tested. The probabilistic structural failures of varying model scenarios are estimated using the First-Order Reliability Method (FORM). Drake Formation shale shows a considerably low probability of failure (~0) with a high reliability index in the initial stress-state condition and after-injection scenarios. Moreover, the parameter sensitivity study indicates that horizontal stress and cohesion are the most influential input parameters during reliability estimation. Comparative analysis between the caprock properties and failure probability reveals that rock strength properties such as cohesion and friction angle strongly dictate the probability of failure estimation. Moreover, comparing two caprock shale formations indicates that the structural failure values are not correlatable; hence, a formation-specific failure assessment is recommended.

1. Introduction

The planned Aurora CO₂ storage site is located in the Horda Platform (south of Troll Field), in the northern North Sea, and was selected as a possible geological CO₂ storage site under the Longship full-scale CCS (carbon capture and storage) project. Longship is the first project of its kind, where a gigaton-level of CO₂ will be injected into the subsurface saline aquifers. Although the subsurface geological CO₂ sequestration has already been demonstrated as a safe and reliable solution by several pilot projects worldwide (i.e., Snøhvit, Norway, [1]; In Salah, Algeria, [2]; Sleipner, Norway, [3]; Ketzin, Germany, [4]; and Otway, Australia, [5]), the Longship project needs to evaluate potential geomechanical failure risks for safe CO₂ storage. Caprock integrity is a crucial parameter to assess during subsurface CO₂ injection because the caprock acts as a top seal and prevents the upward movement of the injected fluids. Mechanical fracturing by either shearing or tensile stresses can occur within caprock shale when reservoir pore pressure exceeds the formation rock fracture strength stress [6]. The risk can be even higher in any CO₂ injection project (i.e., the Longship project), because injecting CO₂ into a saline aquifer will produce an increase in the induced reservoir pressure and influences the mechanical rock deformation and failure. Therefore, caprock’s reliability assessment is crucial to prevent unexpected leakage risks of injected CO₂ and in assisting project decisions.

Caprock shales or mudstones mainly consist of clay and silt-sized particles and significantly differ from other clastic rocks in composition, fabric, heterogeneity, etc. [7,8,9,10]. Caprock shales can deposit in a broad range of environments (i.e., floodplains, lakes, shorefaces, prodeltas, abyssal plains, etc.), and based on the provenance, the type of clay minerals vary significantly [11]. These variations in paleodepositional conditions influence the diagenesis processes [12] and inherent uncertainties in subsurface rock properties. Moreover, laboratory-based caprock strength properties data are limited due to the lack of coring from the caprock section during drilling. Therefore, the number of input parameters is limited, adding additional uncertainties on top of the inherent one. The 3D seismic inverted properties might reduce this data uncertainty range compared to limited wells, and laboratory test data provide better control on lateral rock property variability.

It is highly unlikely that the Drake caprock shales will fail due to capillary breakthroughs because they have high percentages of clay minerals with considerable thickness [6,13]. Instead, injection-related top seal shear failure and tensile fracture are the main likely caprock failure risks in the Aurora storage site. The poroelastic effect of the injected fluid will also be influenced by the mechanical behavior of the rock and the stress path of any injection site [14,15,16,17]. The elastic behavior of the pore fluid under the drainage condition (i.e., fluid flow behavior) is also unknown in a saline aquifer due to no production history, which adds more uncertainty to the system. Generally, the local stress-state condition (i.e., the principal stress values and directions) is also very complex and challenging to assess.

The deterministic assessment of caprock integrity is somewhat questionable, considering the number of uncertainties that are present [18,19,20]. The deterministic approach might oversimplify the complex behavior of caprock as a top seal by assuming the high certainty and constant one-scenario-based input parameters. Generally, all the variables in subsurface structures are stochastic to a certain degree. A probabilistic approach can be more suitable, where the input uncertainty range can be used as an input with a probabilistic failure output [21,22]. The failure probability of an event is a numerical measure of the chance of its occurrence. This outcome will help in choosing the design alternative to improve structural reliability and minimize the failure risk [19]. Structural reliability (SR) can be illustrated by the ability of a structure that comprises given requirements under specific conditions during its intended design life [19,20,23]. In the 19th century, reliability theory was originally developed by maritime and life insurance companies to predict the probability of death for a given population. The intention was to estimate the profitable rates to charge customers. Biological organ failure is in many ways similar to structural system failure (i.e., aircrafts, cars, ships, bridges, etc.). Therefore, the same concept (i.e., reliability methods) were rapidly implemented in various applications in the multidisciplinary design environment, mainly in the engineering realm. This is because the engineering system requires stringent performance, a narrow safety margin, liability, and market competition [19,20].

In the engineering domain, structural reliability (SR) generally deals with load (S) and strength (R) in terms of forces, displacements, and stress acting on the structures. The SR theory is the estimation of failure probability considering the uncertainties originating from poor knowledge of design quantities, such as the likelihood of events, lack of knowledge, variability, degree of belief, inaccuracy, etc. [19]. These uncertainties are related to the loads and the strength of the studied structure. The failure mode can be estimated from the limit-state function, which is defined as the specific limit (borderline surface) between the safe and failure states. There are two types of general limit-state categories: (1) ultimate limit states indicate a structural collapse of part or all (i.e., corrosion, fatigue, deterioration, fire, plastic mechanism, progressive collapse, fracture, etc.); (2) serviceability limit states are related to disruption of the normal use of the structures (i.e., drainage, leakage, excessive deflection, excessive vibration, local damage, etc.) [19].

Based on the strength (S) and load (R) random variables, Cornell [24] proposed the relation for estimating the safety index or reliability index. This approach is a mathematical optimization problem to find the point on the structural response surface (i.e., limit-state approximation) that has the shortest distance from the origin to the surface in the standard normal space. Hasofer and Lind [25] improved the method and introduced a geographic interpretation of the reliability index, which improved the failure probability by transforming the limit-state function into the so-called standard space. Based on the response surface (g(u)) approximation at the most probable failure point, three (3) different methods were projected: (1) First-Order Reliability Method (FORM), (2) Second-Order Reliability Method (SORM), and (3) Higher-Order Reliability Method (HORM).

Structural reliability (SR) can be achieved by several methods (i.e., Monte Carlo, directional sampling method, FORM, SORM, numerical integration, machine learning approach, etc.) that have the solutions for different structural problems (i.e., slope stability, underground mining, fault reactivation, caprock integrity, etc.). For instance, the reliability of the slope of a dam/embankment has been analyzed by several authors using various methods [18,26,27,28,29,30]. At the same time, Nomikos and Sofiamos [31] introduced a new approach to the underground mining industry. Neural-network-based (machine learning) structural reliability estimation methods have been practiced for several decades [32,33,34,35]. In offshore marine structural design (i.e., offshore jacket structure, wind turbines, shipbuilding industry, etc.), reliability methods are also widely used [20,36,37]. The structural reliability assessment of containment for CO₂ sequestration (CCS) is relatively new. Pereira [38] introduced SR methods during the probabilistic estimation of fault reactivation in CCS injection-related pressure perturbation. Rahman et al. [39] estimated the fault system reliability using the event tree methods. Moreover, the same method was used for caprock integrity analysis [40]. The present study is a continuation of the earlier caprock characterization approach where the structural reliability of the Drake caprock shale was evaluated by using a Mohr–Coulomb-failure-criteria-based analytical model. The probability of shear failure and tensile fracture of Drake shale in different conditions (i.e., in situ stress state, injected, and theoretical failure) were estimated. The effect of caprock mineralogy and strength property (i.e., cohesion) were investigated. As this is a new approach for subsurface caprock characterization, a standard rock assessment chart needs to be established. In this study, an initial caprock failure assessment chart is proposed based on the available database such as Drake Formation (this study) and Draupne [40] Formation caprock shale reliability.

2. Study Area

The studied Aurora area is located south of the Troll Fields in the Horda Platform area, in the northern North Sea (Figure 1a). This area experienced two main rifting events during the Permo-Triassic and the Late Jurassic to Mid-Cretaceous times [41,42,43]. The first rifting event was centered in the Horda Platform area with deep-rooted faults and thick syn-depositional wedges, while the second event shifted westward into the deep basin with less influence on the study area [44]. However, the major faults (i.e., Tusse, Svartalv, Troll, etc.) were reactivated during the second rifting event and affected the paleodeposition [41,42,43,45]. N-S trending of those major faults is now demarcating the structural elements in the Horda Platform (Figure 1a). For instance, the studied Aurora area is bounded by Tusse and Svartalv Faults in the east and west, respectively (named the Tusse–Svartalv Fault block). Several NW-SE oriented minor faults are also present, where the number of faults increases towards the north (Figure 1c). Moreover, the overall structure is dipping in the southward direction.

The studied Drake caprock shale is part of the Dunlin Group, deposited in the Early Jurassic post-rift time (Figure 1b). The Drake Formation mainly consists of marine shales deposited in the distal part of the basin [45]. Coarser grain basin margin deposits are also found [46]. Based on the lithological variation, the Drake Formation is divided into two parts, where the upper unit consists of heterolithic deposits comprising sandstones, alternating with siltstone and claystone, while the lower unit mostly consists of claystone [47]. In this study, the silty upper Drake shale is named the upper unit, while the shaly lower Drake is called the lower unit. The main reservoir rocks in the Aurora injection site are Cook and Johansen Formation sandstones, which were also deposited in the Early Jurassic time. Although the upper Amundsen/Burton Formation is present between Cook and Johansen Formations, both reservoir sandstones are in communication due to erosion/non-deposition of the clay layers (Figure 1b). Acoustic properties (i.e., V_p and V_s) and density of upper and lower Drake shales show a considerable difference in well 31/5-7.

3. Materials and Methods

Top seal integrity depends on the geomechanical properties of shales and the stress-state condition of the study area. This study assesses Drake caprock’s structural integrity by using an analytical model defined by the Mohr–Coulomb failure criterion. We wanted to keep our model simple; hence, Mohr–Coulomb failure criterion has been implemented. However, we have a plan to test other criteria (i.e., Hoek–Brown, Barton, etc.) in the future. The workflow used in this research is illustrated in Figure 2. Input parameters (i.e., stress state, pore pressure, and rock strength) were scouted from the published database [48] and estimated using wireline logs from well 31/5-7 and seismically inverted property from the 3D seismic survey GN10M1. The Mohr–Coulomb-failure-criterion-based deterministic factor of safety (FoS) equation was used to define the limit-state function for the probabilistic model. The failure probability of the Drake caprock shale was estimated using the First-Order Reliability Method (FORM) proposed by Hasofer and Lind [25]. The relative importance of the input parameters was also assessed. A detailed description of each section is discussed below.

3.1. Model Parameters

3.1.1. In Situ Stress-State Condition

Leak off test data in the study area and Gassnova report [48] reveal that the horizontal stress gradient is lower than the vertical stress, which indicates a normal faulting regime (Figure 3). We assume isotropic horizontal stress conditions to make our model simple. The hydrostatic and overburden pressure gradients are calculated using the density logs from the well 31/5-7. The water depth at the well location is 307 m. The horizontal stress gradient is estimated using the leak off test data scouted from the Gassnova report [48]. From the report, the base case minimum horizontal stress profile is defined based on the lower bound of the majority of the LOTs, which is adapted as the minimum horizontal stress in this study.

3.1.2. Rock Strength Properties

Laboratory-based caprock strength properties such as cohesion and friction angle are limited due to the unavailability of core in the overburden section. However, the seismic-inversion-based 3D property cubes can be used to resolve the data limitation issue and have spatial control over it by estimating acoustic and elastic properties. The pre-stack simultaneous inversion method is used to estimate the P-wave velocity (${\mathrm{V}}_{\mathrm{p}}$), S-wave velocity (V_s), and density (ρ). Poisson’s ratio ($\mathsf{\nu }$) cube is calculated using the equation below.

$$\mathsf{\nu }=\frac{{\mathrm{V}}_{\mathrm{p}}^2-2{\mathrm{V}}_{\mathrm{s}}^2}{2\left({\mathrm{V}}_{\mathrm{p}}^2-{\mathrm{V}}_{\mathrm{s}}^2\right)}$$

(1)

The inversion methods are not described here, but readers are advised to read Fawad et al. [49] for further information. The referred literature described seismic inversion methodology using a database from the Smeaheia area. The same workflow was also implemented to estimate properties in the Aurora site, which was used in this study.

The ${\mathrm{V}}_{\mathrm{p}}$ cube is converted into unconfined compressive strength (UCS) using the Horsrud et al. [50] proposed empirical equation (Equation (2)).

$$\mathrm{UCS}=0.77{\left({\mathrm{V}}_{\mathrm{p}}\right)}^{2.93}$$

(2)

Cohesion (${\mathrm{S}}_0$) of rocks is estimated from UCS using the following equation:

$${\mathrm{S}}_0=\frac{\mathrm{UCS}}{2}$$

(3)

Top upper Drake, top lower Drake, and top Cook horizons were interpreted using the post-stacked 3D seismic survey GN10M1. The horizons were then converted into time surfaces and used to estimate the average cohesion strength of the upper Drake (from the top upper Drake to top lower Drake surfaces) and lower Drake (from the top lower Drake to top Cook surfaces) units. The distribution of cohesion property of upper and lower Drake shales is illustrated in Figure 4. Spreading of the frequency indicated normal to log-normal distributions. In contrast, lower Drake shale is slightly less cohesive than upper Drake shale. The average and standard deviation values illustrated in the figures are used as the model input parameter.

3.1.3. Horizontal Stress Path

Stress approximation in any sedimentary basin is essential to define stress-state conditions and estimate the effect of poroelasticity. Compaction-induced changes within total vertical stress are minimal for laterally extensive reservoirs, assuming equal elastic properties for the reservoir and overburden. Hettema et al. [51] show in a numerical solution that the stiffness contrast has a negligible effect on the vertical total stress change in a laterally extensive reservoir. However, the horizontal stress path is very sensitive to pore pressure changes. An explanation is that the earth’s surface is a free surface where strain is allowed to absorb any change in pore pressure, while there are constraints in lateral strains [52]. Therefore, the change in horizontal stress is proportional to the change in pore pressure, and the following is stated [51]:

$${\mathsf{\gamma }}_{\mathrm{h}}=\frac{\Delta {\mathsf{\sigma }}_{\mathrm{h}}}{\Delta \mathrm{P}}=\mathsf{\alpha }\left(\frac{1-2\mathsf{\nu }}{1-\mathsf{\nu }}\right)$$

(4)

where $\Delta {\mathsf{\sigma }}_{\mathrm{h}}$ is the change in horizontal stress, $\Delta \mathrm{P}$ is the pore pressure change, $\mathsf{\alpha }$ is Biot’s coefficient or effective stress parameter, and $\mathsf{\nu }$ is Poisson’s ratio.

Considering that Biot’s coefficient is equal to 1, Equation (3) can be rewritten:

$${\mathsf{\gamma }}_{\mathrm{h}}=\left(\frac{1-2\mathsf{\nu }}{1-\mathsf{\nu }}\right)$$

(5)

Using the seismic inverted Poisson’s ratio cube, the horizontal stress path volume was estimated. The frequency distribution of Poisson’s ratio (ν) of upper and lower Drake units also indicated normal to log-normal spreading (Figure 5). The average ν slightly varied between the upper and lower shale units, 0.3 and 0.32, respectively. The average and standard deviation of stress path values are estimated similar to the cohesion strength property that has been defined.

3.2. Model Definition

3.2.1. Reliability Analysis and Limit-State Function

The probability of failure does not characterize the deterministic structural safety index such as fail or non-fail; instead, it represents the failure likelihood. Hence, it indicates the reliability of any structure. It provides a rational framework for dealing with uncertainties of structural parameters and assists in making the project decision. The reliability of a structural component depends on active load (S) and resistance (R) by the structure itself [22]. Uncertainties of S and R are represented as random variables, assuming the load and the resistance variables’ areas are statistically independent. If both the resistance and load variables are normally distributed, the safety margin (random variable M) can be determined by the following:

$$\mathrm{M}=\mathrm{R}-\mathrm{S}$$

(6)

and the probability of failure (${\mathrm{P}}_{\mathrm{f}}$) is assessed through Equation (7), where ${\mathrm{P}}_{\mathrm{f}}$ is determined by the realization of the random variables R and S but not the overlapping of two curves.

$${\mathrm{P}}_{\mathrm{f}}=\mathrm{P}\left(\mathrm{R}-\mathrm{S}\le 0\right)=\mathrm{P}\left(\mathrm{M}\le 0\right)$$

(7)

where M is also normally distributed with the mean ${\mathsf{\mu }}_{\mathrm{M}}={\mathsf{\mu }}_{\mathrm{R}}-{\mathsf{\mu }}_{\mathrm{S}}$ and standard deviation ${\mathsf{\sigma }}_{\mathrm{M}}=\sqrt{{\mathsf{\sigma }}_{\mathrm{R}}^2+{\mathsf{\sigma }}_{\mathrm{S}}^2}$.

Failure probability is determined by the use of the standard normal distribution function shown in Figure 6a, which mathematically states the following:

$${\mathrm{P}}_{\mathrm{f}}=\mathsf{\Phi } \left(\frac{0-{\mathsf{\mu }}_{\mathrm{M}}}{{\mathsf{\sigma }}_{\mathrm{M}}}\right)= \mathsf{\Phi } \left(-\mathsf{\beta }\right)$$

(8)

where ${\mathsf{\mu }}_{\mathrm{M}}/{\mathsf{\sigma }}_{\mathrm{M}}=\mathsf{\beta }$ is called the safety/reliability index, which is the standard deviation by which the mean value of the safety margin M exceeds zero or most likely exceeds the failure point (Figure 6a).

Including uncertainties of the random variables, the probability of failure (${\mathrm{P}}_{\mathrm{f}}$) can be determined from the limit-state function (g(X)), which is a boundary between desired (g(X) > 0) and undesired (g(X) ≤ 0) performance of any structure and can be defined within a mathematical model for functionality and performance [53].

In this study, the Mohr–Coulomb-failure-criteria-based limit-state function is determined from the deterministic factor of safety (FoS) equation. Assuming an isotropic horizontal stress condition within a normal faulting regime, the factor of safety (FoS) for shear failure is defined as the following:

$$\mathrm{FoS}\_\mathrm{shear}=\frac{[\left(\frac{{\mathsf{\sigma }}_1^{\prime }+{\mathsf{\sigma }}_3^{\prime }}{2}\right)+\frac{{\mathrm{S}}_0}{\tan \mathsf{\varphi }}]\sin \mathsf{\varphi }}{\frac{{\mathsf{\sigma }}_1^{\prime }-{\mathsf{\sigma }}_3^{\prime }}{2}}$$

(9)

$${\mathsf{\sigma }}_1^{\prime }={\mathsf{\sigma }}_1-{\mathrm{p}}_{\mathrm{p}}$$

(10)

$${\mathsf{\sigma }}_3^{\prime }={\mathsf{\sigma }}_3-{\mathrm{p}}_{\mathrm{p}}$$

(11)

where ${\mathsf{\sigma }}_1^{\prime }$ is effective vertical stress, ${\mathsf{\sigma }}_1$ is vertical stress, ${\mathsf{\sigma }}_3^{\prime }$ is effective horizontal stress, ${\mathsf{\sigma }}_3$ is horizontal stress, ${\mathrm{p}}_{\mathrm{p}}$ is pore pressure, S₀ is cohesion, and ϕ is friction angle.

The state of the structure is safe when the factor of safety is greater than 1 and fails when it is less than 1. Therefore, the limit-state function is defined as the following:

$$\mathrm{g}\_\mathrm{shear}\left(\mathrm{X}\right)=\mathrm{FoS}\_\mathrm{shear}-1$$

(12)

where g(X) is the limit-state function, which is the boundary between the safe (g(X) > 0) and failure (g(X) ≤ 0) state. This limit-state function only evaluates the shear failure and fracture risks. There might be a tensile failure risk that must be considered when defining the probabilistic failure limit. Therefore, the tensile failure risk limit-state function is also defined as the following:

$$\mathrm{g}\_\mathrm{tensile}\left(\mathrm{X}\right)={\mathsf{\sigma }}_3-{\mathrm{p}}_{\mathrm{p}}$$

(13)

Equation (13) represents the tensile failure based on the Mohr–Coulomb failure envelope, which is a pessimistic estimation. However, considering the Mohr–Coulomb–Griffith failure envelope, we tested another two scenarios with a tensile strength cut-off. The tensile strength cut-off is considered half of the cohesion of that rock sample.

3.2.2. Reliability Method

The studied model is executed using the method called First-Order Reliability Method (FORM), which is proposed by Hasofer and Lind [25] and widely used in practical engineering problems [22,54]. Let us assume two random variables with the limit-state function g(X) are stated in physical space, where the grey zone represents the failure events and the white zone is the safe region (Figure 6b). FORM normalizes the random variables into normally distributed variables (Z₁ and Z₂) with zero mean and defines the limit-state function g(Z) = 0. After that, the g(Z) surface is linearized at point z*, which is the shortest distance from the zero mean and is called the design point (Figure 6c). The shortest distance is called the reliability index (β), which is an indicator of structural stability. The higher the distance, the safest the structure is. The normal vector direction to the failure surface is denoted as α and indicates the random variable’s relative importance. Mathematically the reliability index (β) can be expressed as the following:

$$\mathsf{\beta }= \overrightarrow{\mathsf{\alpha }}{\mathrm{z}}^{*}$$

(14)

where the normal vector to the failure surface $\overrightarrow{\mathsf{\alpha }}$ is denoted as

$$\overrightarrow{\mathsf{\alpha }}=\frac{\nabla \mathrm{g}\left({\mathrm{z}}^{*}\right)}{\left|\nabla \mathrm{g}\left({\mathrm{z}}^{*}\right)\right|}$$

(15)

where g(Z) is the gradient vector and assumed to exist:

$$\nabla \mathrm{g}\left(\mathrm{Z}\right)=\left(\frac{\partial \mathrm{g}}{\partial {\mathrm{z}}_1}\left(\mathrm{Z}\right),\dots \dots \dots \dots \dots .,\frac{\partial \mathrm{g}}{\partial {\mathrm{z}}_{\mathrm{n}}}\left(\mathrm{Z}\right)\right)$$

(16)

Therefore, an iterative method can be used to estimate the reliability index due to the non-linear optimization problem [55,56]. The Python-based First Order Reliability Method (FORM) can iterate and estimate the failure value with several iterative steps [39,40].

The method (i.e., FORM) is initiated and run by using a python-based open-source structural reliability analysis module called ‘PyRe’ [57], which was created using the core function of the Finite Element Reliability Using Matlab (FERUM) project started in 1999 at the University of California, Berkeley, for pedagogical purposes. In this study, only the core reliability functionality was used, which focuses on the analytical solution of the structural reliability assessment. The 2019 version of Petrel^TM was used for seismic interpretation and property extraction purposes.

3.3. Modeling Scenarios and Input Parameters

Based on various assumptions, six different scenario-based models were tested, where both upper and lower Drake units were analyzed (

Table 1. Different model scenarios were tested in this study to evaluate Drake caprock shale’s reliability.

	Unit	Assumptions
Case-1	Upper Drake	Initial stress-state condition
Case-2	Lower Drake
Case-3	Upper Drake	After the CO2 injection scenario
Case-4	Lower Drake
Case-5	Upper Drake	Theoretical shear failure due to decreasing σ3
Case-6	Lower Drake

). Case-1 and Case-2 are defined based on the initial stress-state condition, while Case-3 and Case-4 represent the pressure build-up scenario after CO₂ injection into the saline aquifer adapted from the Gassnova [48] reservoir simulation model. The Gassnova model used an injection rate of 3.2 MT per year over a period of 50 years, which accounted for a total of 160 MT of CO₂. This leads to a reservoir pressure increase of 3.2 MPa. The average permeability and porosity in the injection location are about 625 mD and 20%, respectively. The injection-induced pressure increase near the well location is estimated at 3.2 MPa.

Two theoretical failure scenarios were also analyzed, where the horizontal stress decreases up to the failure point while the other parameters (i.e., σ_v and P_p) remain in their initial condition. Please note that the shear failure horizontal stress is defined based on the Mohr circles and Coulomb–Griffith failure envelope. However, considering the tensile stress cut-off at zero (i.e., Coulomb failure envelope), we are within the tensile failure zone (i.e., horizontal stress is lower than the pore pressure). Therefore, different probabilistic tensile failure risk sensitivity was also evaluated using the limit-state function stated in Equation (13). Tensile failure was estimated only for Case-5 and Case-6, because considerable differences were observed between σ₃ and P_p in Case-1 to Case-4. Note that Case-5 and -6 are not real scenarios; instead, they are theoretical failure cases where σ₃ decreases hypothetically until achieving shear failure. When σ₃ reduced to a certain level closer to the P_p value, there was a possibility of tensile failure; hence, we tested tensile failure risk as well. Moreover, these theoretical failure models (i.e., Case-5 and -6) of failure probability are necessary to define the caprock failure assessment chart, which is not available for the subsurface characterization yet.

Stress parameters in the studied section were estimated based on well 31/5-7, considering this well as the potential injection point. The input database used in this study is shown in

Table 2. The average and standard deviation of the model input parameters are presented with the type of distribution and data sources. Please note that the superscript numbers in the property names represent the case numbers explained in Table 1.

Parameters	Average	Unit	Standard Deviation	Random Distribution	Sources
Initial vertical stress (σv)	50.50	MPa	3.21	Normal	[48]
Initial horizontal stress (σh1–4)	36.50	MPa	1.33	Normal	[48]
Initial horizontal stress (σh5)	20.53	MPa	0.00	Normal	[48]
Initial horizontal stress (σh6)	21.84	MPa	0.00	Normal	[48]
Initial pore pressure (Pp1,2,5,6)	25.70	MPa	0.28	Normal	[48]
Injected pore pressure (Pp3,4)	28.90	MPa	0.28	Normal	[48]
Horizontal stress path (γh)	0.60	-	0.11	Normal	Inversion cube
Cohesion, upper Drake (S01,3,5)	12.23	MPa	3.98	Normal	Inversion cube
Cohesion, lower Drake (S02,4,6)	11.26	MPa	3.06	Normal	Inversion cube
Friction angle (ϕ1–6)	21.63	Degree	2.33	Normal	[58,59]

. The properties also reflect spatial variability because of the use of seismic inverted properties while estimating the average and standard deviation of the horizontal stress path and cohesion strength. Stresses (i.e., vertical and horizontal) and pore pressure were calculated using the gradient illustrated in Figure 3. To our knowledge, there are no friction angle data available for Drake caprock shale; hence the Draupne caprock shale value [40] was used. The random variables such as vertical stress (σ_v), horizontal stress (σ_h), pore pressure (P_p), cohesion (S₀), and friction angle (ϕ) were used to run the stochastic model, where the arithmetic average with standard deviation was used to define the ranges. Three sigma rules ((maximum value-minimum value)/6) were implemented during estimating standard deviation, which covers 99.73% of all values of a normally distributed parameter [18]. Moreover, most geological processes follow a normal or log-normal distribution [21], which is also revealed in Figure 4 and Figure 5; hence, the normal random distribution of the property was used in the studied models. In addition, standard deviation serves as a measure of uncertainty, where a small value indicates clustered closely around the mean with more precision and vice versa.

3.4. Parameter Sensitivity

The OVAT (one variable at a time) technique [60,61] was used to estimate the deterministic sensitivity by alternatively assigning minimum and maximum values in each input parameter when other parameters remain in their mean values. Usually, the output is illustrated as a tornado diagram, where the relative influence of each parameter can be assessed. The sensitivity output from this study is shown in Figure 7, where the cohesion (S₀) has the most significant impact on the factor of safety in both upper and lower Drake units. In situ stresses (i.e., vertical and horizontal) also have considerable influence. Pore pressure itself has a negligible impact. A similar parameter impact trend was observed between upper (Case-1) and lower (Case-2) Drake units. Only the initial stress-state condition scenarios (Case-1 and -2) are presented here as an example.

4. Results

The Drake caprock shale structural reliability (based on shear failure) for different pressure conditions is summarized in

Table 3. The deterministic factor of safety (FoS) and the probability of failure (P_f) of the Drake caprock shale in different studied cases. Corresponding reliability index (β) values are also shown.

Model Scenarios	Factor of Safety	Probability of Failure	Reliability Index
Case-1 (UD)	2.57	0.00	10.80
Case-2 (LD)	2.44	0.00	11.08
Case-3 (UD)	2.40	0.00	9.72
Case-4 (LD)	2.27	0.00	10.00
Case-5 (UD)	1.00	1.10 × 10−2	2.29
Case-6 (LD)	1.00	2.52 × 10−4	3.48

. The corresponding reliability index (β) and deterministic factor of safety (FoS) for shear failure values are also presented. The estimated probability of failure (P_f) is close to zero in both in situ stress conditions (Case-1 and -2) and after-injection scenarios (Case-3 and -4). Based on the reliability index, the lower Drake unit represents higher safety than the upper Drake shale in both scenarios (i.e., before and after injection). However, the calculated deterministic factor of safety values demonstrates the opposite trend: the upper Drake shale is safer than the lower Drake unit. In addition, the reliability index in the after-injection scenario slightly decreases than the in situ stress-state condition.

Drake caprock shale’s theoretical shear failure scenarios (Case-5 and -6) represent significant variation between upper and lower shale units. The lower Drake shale has a substantially lower failure probability (2.52 × 10⁻⁴) than the upper unit (1.10 × 10⁻²), with both having the same deterministic factor of safety value of 1. Compared with the other cases, the trend between the upper and lower units also repeats (i.e., higher structural reliability in lower Drake than the upper unit) in theoretical failure scenarios. The corresponding reliability index in failure cases also shows considerable variation.

4.1. Tensile Failure

Tensile failure risk was also evaluated for theoretical failure scenarios (i.e., Case-5 and -6) using the same reliability method. Minimum (i.e., Coulomb failure envelope) and maximum (i.e., Griffith tensile extension) tensile stress cut-offs were assessed, where the minimum cut-off used is zero (0), and the maximum cut-off used is half of the cohesive strength of those specific intervals (Table 2). The initial stress state and injected scenarios were not analyzed, assuming there will be no tensile failure risk because the value of horizontal stress is significantly higher than the pore pressure even in the after-injection scenario.

The theoretical scenario models illustrated the high possibility of failure even in the maximum tensile stress cut-off cases (

Table 4. Tensile failure and reliability index of Case-5 and -6 in different tensile stress cut-off values.

Model Scenarios	Tensile Stress Cut-Off	Probability of Failure	Reliability Index
Case-5 (UD)	0	9.99 × 10−1	−3.89
	6.115	2.39 × 10−1	0.71
Case-6 (LD)	0	9.98 × 10−1	−2.9
	5.13	1.70 × 10−1	0.95

). Based on the Coulomb failure cut-off, the failure probability is ~100% with a negative reliability index value. However, using the Griffith tensile extension cut-off (half of the cohesion), the failure probability relatively increases with a positive reliability index value. Comparing the upper and lower Drake shale results, the lower Drake shale is slightly less failure-prone than the upper unit. Moreover, comparing the result between shear failure (Table 3) and tensile failure (Table 4) of Case-5 and -6 reveals that the tensile failure probability is significantly higher than the probabilistic shear failure in the case of caprock failure.

4.2. Probabilistic Sensitivity

The normal vector of the failure surface α (Figure 6 and Equation (15)), which is the relative design sensitivity factor or the relative parameters’ importance factor, is often referred to as probabilistic sensitivity and is very useful for the relative ranking of random variables. The vector can be obtained during performing the probabilistic analysis by treating every individual input parameter as a deterministic variable, where a positive value indicates a direct relationship between the variables and the response and a negative sensitivity suggests an inverse relationship [38,62].

The sensitivity factors of each case are illustrated in Figure 8. The effect of the vertical stress (σ_v) in all cases is close to zero. Pore pressure (P_p) and friction angle (ϕ) show a direct relation with the reliability value, while horizontal stress (σ_h), horizontal stress path (γ_h), and cohesion (S₀) are indicated as an inverse relation. The maximum influences during the in situ stress-state condition (Case-1 and -2) and after-injection scenarios (Case-3 and -4) in both upper and lower Drake units are represented by σ_h and S₀, while a significant negative relation is observed between S₀ and theoretical failure structural reliability (Case-5 and -6). In addition, overall, the variation between the upper and lower Drake units is negligible.

The input parameter contribution to the reliability function is represented by the square of the sensitivity factor (α²), where the sum of all the parameters is equal to 1 [63]. Because of the similar trend between the upper and lower Drake units, only the lower Drake unit (the primary caprock in the Aurora site) model contributions are illustrated in Figure 9. The probabilistic failure values of Case-2 and -4, which are estimated using the FORM technique, mainly depend on the horizontal stress, 76% and 70%, respectively. A substantial influence from cohesion is also observed. However, the cohesion strength property is controlled in the theoretical failure scenario (Case-6) by a 94% impact on overall structural failure.

5. Discussion

Failure probability of subsurface rock stability could be a possible approach to deal with highly uncertain input parameters [22], where the random variable method included the entire range within the final probabilistic output. The probability of failure of the Vette fault and Draupne caprock shale indicated this method’s practicality during quantifying the failure risks [39,40]. However, the difference in deterministic and probabilistic sensitivities within the previous analysis indicated the limitation of the number of input parameter samples and the method used to define the standard deviation. The probabilistic model is also susceptible to the input parameter ranges; hence, emphasis is needed to define the uncertain parameter range before use as an input parameter in the failure probability estimation. Additionally, statistical analysis with more data points is needed to define the input parameters’ range to obtain a field-scale level reliable assessment. In this study, using seismic inverted properties with the six-sigma rule to define standard deviation significantly improves the parameter’s uncertainty range. A consistency between deterministic and probabilistic input parameter sensitivities reflects the improvement of the studied method. According to the analysis, the most influential parameters are horizontal stress and cohesion strength of the rocks (Figure 8 and Figure 9). However, the most influential parameter changes dramatically in theoretical failure scenarios (Case-5 and -6), where the contribution of cohesion is 94% (Figure 9).

Cohesion of the upper and lower Drake shales is high; hence, the effective stress Mohr circle is away from the Coulomb failure envelope (Figure 10). As a result, the chances of shear failure are close to zero in the initial stress-state condition and injected scenarios (Case-1 to -4). FOAM estimation is also dictated by horizontal stress and cohesion (Figure 9). However, when we forced ${\mathsf{\sigma }}_3$ to have a theoretical shear failure (Case-5 and -6), the value declined below pore pressure. In this scenario, we not only have the increased shear failure but also the possibility of increased tensile failure risks. This explained the ~100% tensile failure probability when using the tensile cut-off of zero (0) (Table 4).

5.1. Effect of Caprock Properties

The effect of rock properties on probabilistic failure was assessed by comparing the studied upper and lower Drake shales with shallow organic-rich Draupne caprock shale. Cross plots of the reliability index and caprock shale properties (cohesion and bulk clay mineralogy) are illustrated in Figure 11. Draupne caprock shale results are adapted from Rahman et al. [40,64]. The clay mineral fraction within the Draupne and Drake caprock shales varies significantly, where the average Draupne clay mineral percentage is 65%, and 40% and 20% represent the lower Drake and upper Drake, respectively. Although the difference in cohesion strength of the upper and lower Drake shale is minimal, a considerable variation with Draupne shale was observed. Furthermore, the studied Draupne Formation is located within a mechanically compacted zone, while the Drake Formation is chemically compacted [64,65,66]. The reliability indexes of upper and lower Drake shales are significantly higher in initial stress-state conditions than Draupne shale, irrespective of the percentage of soft clay minerals (i.e., ductile minerals). Moreover, a minimal reliability index difference is also observed between upper and lower Drake units, though significant variation was identified in bulk clay mineralogy. This indicates that the probabilistic failure value does not directly correlate with ductile clay mineral percentages but depends on the caprock strength property such as cohesion.

On the contrary, the reliability index has a positive correlation with cohesion strength. However, compared with the parameter sensitivity analysis, the relative influence of the cohesion property varied significantly between Draupne and Drake caprock shales ([40]; Figure 9). Based on the probabilistic sensitivity, this reveals that caprock integrity in the mechanically compacted shale (i.e., Draupne shale) is more affected by principal stresses. In contrast, for the deeper caprock (Drake shale), the rock strength (i.e., cohesion) is a more influential parameter.

The factor of safety in both Draupne and Drake caprock shales has a similar range, while the reliability index dramatically increases in Drake shale (Figure 12). The differences in probabilistic failure numbers between Draupne and Drake caprock shales indicated the importance of the probabilistic approach over the deterministic while assessing the top sealing reliability. In addition, this specifies that the relationship between deterministic and probabilistic safety factors varied within different caprock shales. Different β values in failure scenarios indicate that the estimation of the probability of failure depends on parameter uncertainty, and the boundary between safe and failure conditions has a range instead of a constant number.

The probabilistic assessment of top seal integrity facilitates decision-making in uncertainty by explicitly showing the trade-off between the investment and the reduction in potential failure risk [67]. In connection with the CO₂ storage site, the probability of top seal leakage risk might simplify future project decisions.

5.2. Initial Safety Chart

The reliability index and probability of failure in any structure are relative measurements of the current condition and provide a qualitative estimation of the expected performance [29]. However, top seal integrity is complex and depends on many factors such as pressure, temperature, mineralogy, etc. There is no standard chart for subsurface rock structural reliability to compare, because this method is a novel subsurface structural reliability assessment approach. An initial safety boundary is proposed using the theoretical shear failure result from Draupne and Drake caprock shales (

Table 5. The performance chart with the unsatisfactory events and reliability index ranges are modified after [29]. The color represents the proposed subsurface structural safety boundary based on analysis.

Expected Performance Level	Probability of Unsatisfactory Performance	Reliability Index (β)	Proposed Safety Boundary
High	0.0000003	5.0
Good	0.00003	4.0	Safe
Above average	0.001	3.0	Boundary
Below average	0.006	2.5
Poor	0.023	2.0
Unsatisfactory	0.17	1.5	Unsatisfactory
Hazardous	0.16	1.0

). The boundary between safe and unsatisfactory performance is proposed based on reliability index values 2 to 3, based on the range obtained for theoretical failure cases. However, the proposed safety chart is not directly comparable with the published expected performance level for soil [29]. This chart is a proposal, but we need further work to use it for quantitative assessment. In addition, the proposed safe and unsafe zones are not valid for tensile failure probability. Nevertheless, a tensile failure estimation is highly recommended if the CO₂ injection reservoir pressure increases to a range close to the horizontal stress value.

This study reveals that the probabilistic approach used to assess structural reliability should be convenient and robust to capture the risks related to the uncertainty of the parameters. FORM seems to be one of the appropriate methods for static subsurface structural reliability assessment [55]. Therefore, the workflow used in this study can be effective and useful to evaluate the shared structural failure risks in any subsurface injection project, especially CO₂ injection site evaluation.

6. Conclusions

The reliability of the top seal is critical for any permanent storage of CO₂ into the subsurface. Therefore, Drake caprock shale’s low failure probability is crucial for safe CO₂ storage in the Aurora site. The main findings of this study are stated below:

▪

The Drake caprock shale failure chance is zero in both the in situ stress-state condition and after-injection scenarios. However, this method indicated a qualitative assessment and requires further investigation to quantify the results.

▪

Horizontal stress and cohesion are the most influential input parameters when estimating probabilistic failure.

▪

Shallow caprock shale units are more affected by the principal stresses, while the rock strength properties influence the deep caprock shale units while estimating the failure probability.

▪

A subsurface structural reliability safety chart is proposed to evaluate the safe and unsatisfactory caprock shales. However, as the failure probability values are formation-specific, a standard universal chart might be impractical and requires further investigation.

▪

Tensile failure estimation should be carried out if the injected pore pressure and horizontal stress have a similar range.

Although the Drake Formation shale indicates it is less likely to fail, the injection-related potential risks can be affected by other factors (e.g., stiffness contrast between reservoir and surroundings, geometrical effects, drainage condition, stress paths, etc.) and require further evaluation using a numerical simulation.