Characterizing along- and across-fault fluid-flow properties for assessing flow rates and overburden fluid migration along faults: a case study from the North Sea

Implementation of large-scale CO₂ storage will require utilization of a wide range of storage reservoirs, including faulted reservoirs with structural traps. Structural traps have been proven to withhold sizeable hydrocarbon columns (Spencer and Larsen 1990; Knott 1993; Yielding et al. 1999, 2010; Jolley et al. 2007; Osmond et al. 2022); however, to ensure safe utilization of structural traps in saline aquifers for CO₂ storage, fault risk assessment needs to be extended from the reservoir across fault (Manzocchi et al. 1999; Yielding et al. 2010) to include the potential for fluid migration up along faults and into the overburden. The Smeaheia fault block in the Horda Platform, offshore Norway is bounded by the Vette Fault Zone (VFZ) to the west (Fig. 1). Along the VFZ, the Alpha closure in the Upper Jurassic succession (Fig. 1) has been identified for additional storage volumes as part of the upscaling of the Norwegian Longship carbon capture and storage (CCS) project (Wu et al. 2021). To develop and qualify the Alpha structure for CO₂ injection, the containment must be verified, and recent work has focused on describing the Alpha closure in detail and identifying potential risks in the Smeaheia fault block (Mulrooney et al. 2020; Wu et al. 2021). A comparison of the structural setting for the Smeaheia fault block and cross-fault sealing along the reservoir closures in the Troll oil and gas field indicates that a similar sealing potential may be expected for the VFZ (Osmond et al. 2020). Demonstrating that the overburden side seal of the Alpha closure along the VFZ exhibits sufficient sealing properties for storing CO₂ is an important step towards qualification of a storage site. The hypothesis for the current model is that combining existing knowledge of the relevant overburden formation permeability and geological variations in the VFZ demonstrates fluid-migration pathways with very limited possibilities for CO₂ migration during static fault conditions.

Faults in siliciclastic successions can provide conduits or barrier to flow, with hydraulic properties dependent on factors such as the fault juxtaposition, internal fault architecture, temperature and stress history, as reviewed by Pei et al. (2015). In the oil and gas industry, special attention has been given to faulted reservoirs and faults as barriers to flow in compartmentalized reservoirs. Several methods have been developed to quantify the potential for clay to reduce the permeability within fault zones compared to the reservoir permeability. Existing models include the shale smear factor model (SSF: Lindsay et al. 1993), the clay smear potential model (CSP: Bouvier et al. 1989; Lehner and Pilaar 1997), the shale gouge ratio model (SGR: Yielding et al. 1997) and extensions of these: simple shear zone (Childs et al. 2007), probabilistic SSF (Childs et al. 2007) and effective SGR (Knipe et al. 2004; Freeman et al. 2010), into fault-seal risk assessments (Færseth et al. 2007). These types of models are typically used to predict the sealing potential for across-fault hydrocarbon migration using correlations between SGR and fault-core permeability (e.g. Sperrevik et al. 2002; Childs et al. 2007), as well as SGR and maximum hydrocarbon column heights (Yielding et al. 2010). Despite this large number of fault models, the implementation of faults in full reservoir simulations is limited, and the most common method is to include the effects of faults on reservoir permeability as transmissibility multipliers (Manzocchi et al. 1999). Only a limited number of studies have tried to fully describe the fault-zone permeability and calculate fluid migration within the fault zone (Bense and Van Balen 2004; Bense and Person 2006; Fredman et al. 2007; Braathen et al. 2009; Fachri et al. 2013; Bjørnarå et al. 2021, 2022). Detailed fault permeability models suitable for direct flow simulations remain elusive. A good workflow for establishing such models needs to incorporate available site-specific data, as well as generic knowledge about faults. Conservative modelling approaches need to be replaced by uncertainty quantification and probabilistic assessments. One of the main obstacles, that of the lack of direct measurements of fault-zone permeability, sparse data from overburden logs and few direct tests, has the potential to be implemented into the existing juxtaposition and mixing models for sand–clay sequences to produce a likely fault permeability. Together with models for throw/fault width (e.g. Bense and Person 2006; Torabi and Berg 2011; Alaei and Torabi 2017), the sparse data can be combined into a realistic overburden fault model with the potential to directly address the flow field in the overburden above a CO₂ storage site and to quantify the fluid-migration potential using models that include uncertainty.

The current study has two objectives. The first objective of this paper is to present a reliable and robust fault characterization method addressing both along- and across-fault flow properties and the resulting leakage risk into the overlying seal formations of the VFZ. The characterization is based on geological data such as seismic imaging, gamma-ray well logs and available petrophysical data for the strata surrounding the fault. There are uncertainties in the data used in the characterization; thus, the second objective is to evaluate the sensitivity and variations in leakage potential due to these uncertainties in a parametric study. In parallel, the leakage potential is also evaluated in a model where the fault characteristics are ignored, and the fault outline represents a discontinuity, or juxtaposition, of the lithologies.

The flow model designed for describing the leakage risk along the VFZ in the overburden of the Alpha closure (Fig. 1) is based on the extensive geological characterization of the area, and the fault interpretation and modelling targeting this fault (Mulrooney et al. 2020; Osmond et al. 2020; Michie et al. 2021b; Rahman et al. 2021; Wu et al. 2021).

The Alpha closure is in the footwall of the north–south-trending VFZ within the Smeaheia fault block in the Horda Platform (Fig. 1a). The Smeaheia fault block is the easternmost rotated fault block in the Horda Platform bounded by the Øygarden Fault Complex and basement to the east (Fig. 1b); whereas to the west, the rotated fault blocks of the Tusse and Svartalv fault zones bound the Troll East and West hydrocarbon reservoirs, respectively. The main reservoir for the Alpha closure is the Sognefjord Formation, with the upper Sognefjord Formation displaying the best reservoir properties (Gassnova 2012; Statoil 2016; Ringrose et al. 2017; Mondol et al. 2018). However, in the flow model, the reservoir section consists of the entire Viking Group, including the Fensfjord Formation and the Krossfjord Formation, which are interlayered with the lower-permeability Heather Formation (Gassnova 2012; Sundal et al. 2014; Mondol et al. 2018). The reservoir is capped by the Draupne Formation, a well-documented sealing unit in the North Sea (e.g. Skurtveit et al. 2012, 2015; Soldal et al. 2021), whereas the juxtaposed Cromer Knoll Group provides the lateral seal for the closure within the VFZ (Fig. 1) and can also form a secondary top seal where thinning of the Draupne Formation occurs. The units below the Viking Group represent the underburden in the model. Since the focus for the flow model is to define the migration pathways for the side and overburden of the CO₂ storage reservoir in the Alpha closure, efforts to map the internal detail within the Viking Group (Sognefjord, Fensfjord and Krossfjord formations) and to distinguish formations in the underburden are not included in this study.

The geometry in the 2D flow-model profile (Fig. 1c) was digitized from the interpreted regional cross-section from 2D seismic profile NNST84-05 (Fig. 1b) (Mulrooney et al. 2020; Michie et al. 2021b) and shows key lithological units. The properties of the stratigraphic units within the model are reported in Table 1. The data are based on lithological descriptions obtained from well logs (wells 32/4-1 and 31/6-6; well locations are shown in Fig. 1a) and from published and unpublished sources. For overburden stratigraphy, the permeability represents low-permeability sections tested for seal capacity (e.g. Worden et al. 2020); hence, these values are used in the model as the lower-bound permeability. To assess the leakage risk related to the fault, the permeability of the units above the cap rock (Table 1) were varied between a lower-bound and an upper-bound permeability, together with parameters that correlate the fault width and the clay content (V_shale) (Sperrevik et al. 2002).

The VFZ throw and juxtaposition have been described from seismic interpretations (Mulrooney et al. 2020; Michie et al. 2021b). The clay content (V_shale) of the host-rock formations has been used to characterize the hydraulic properties and width of the VFZ. In this study, the clay content has been estimated from the gamma-ray (GR) log from well 31/6-6 (Fig. 2b), which is located relatively close to the seismic profile in the hanging wall of the fault (Fig. 1a). The GR log is related to the seismic horizon boundaries in the wells and has been projected onto the horizon boundaries along the VFZ within the model for the hanging wall (Fig. 2c) and footwall (Fig. 2d). The GR log was converted to V_shale using high GR values (223 gAPI; the maximum value in the dataset was 218 gAPI) in the Draupne Formation as the value for 100% clay, and the low GR value (33 gAPI; the minimum value in the dataset was 36 gAPI) in the reservoir section as the value for 0% clay.

The mathematical model for the VFZ is described by a zero-thickness element, or an edge, and the complex architecture is described mathematically. Here, we used the conceptual model, and methodology, proposed by Bense and Person (2006), which was originally designed with a hydrogeology application in mind in near-surface sediments, and is now considered to address fault flow in the overburden of a fault-bounded CO₂ storage reservoir. In their approach, Bense and Person (2006) considered that the fault width and permeability structure varied with the throw D (m) of the fault, and with the clay content V_clay (%) and permeability of the host rock, as described below.

The methodology of Bense and Person (2006) used in this study provides a static fault description and represents an initial characterization of the hydraulic behaviour of the fault. It excludes dynamic changes due to injection and post-injection processes (e.g. pore pressure and stress changes, and strain). An example of the fault flow algorithm is provided in Figure 4 to illustrate the conceptual idea. In the example, a high-permeability layer (sandy yellow, zero clay content) and a 20 m-thick low-permeability layer (shaly dark grey, 100% clay content) is sheared where the left is displaced downwards relative to the right. The equations in the methodology simplify as we have $dwd=dwu$⁠, $kd=ku$ and D_u = 0 m. The throw of the fault is a constant 15 m for illustration purposes only. The initial permeability profile is shown by the red line (relative to the right side of the fault in Fig. 4b). In this example, the sand is given a fault-width growth coefficient value of dw = 1 and the shale is given a value of dw = 0.5; hence, where shale is only juxtaposed against shale, the fault width is 7.5 m and where sand is only juxtaposed against sand, the fault width is 15 m, and where the layers mix the width will be between these two bounds (thick blue line in Fig. 4c). Figure 4d shows the smearing of the shale material, and the magenta line shows the location of the smearing profile in the triangle diagram for a constant throw of 15 m. Note that the initial permeability of the low-permeability layer is 10 times lower than the high-permeability layer.

Using equations (1), (2), (4) and (5) for this simplified example, the along-fault permeability (thick black solid line in Fig. 4b) is slightly higher compared to the clay content weighted average of the mix of the two materials (shown by thin solid back line, equivalent to dw = 1 for both materials), which is associated with the dw value for the low-permeability material being lower than the high-permeability material; thus, the high-permeability value is slightly favoured in the arithmetic mean. The across-fault permeability (thick black dotted line in Fig. 4b) is mainly dominated by the low-permeability values from the harmonic averaging technique. Again, the thick dotted line is slightly higher compared the thin dotted line due to different values of dw.

To quantify the leakage potential along faults into overlying seal formations, a simplified flow model was designed based on the regional cross-section in Figure 1c. To focus the work towards understanding the effects of the fault on fluid migration, the following assumptions and simplifications were made for the flow model.

The flow model was simplified to a 2D representation of the subsurface system under steady-state conditions. A 2D model implies no hydraulic variation in the direction perpendicular to the 2D model cross-section. Real systems may have geological heterogeneities in this direction, with pore pressure/fluid flow relief that will affect the flow field and is not captured in a 2D model. Because we assumed steady-state conditions, the timescale of the system was ignored and transient processes were not considered. For instance, due to the low permeability of most of the lithological groups and formations (see Table 1), it will take a long time for the pore pressure to dissipate and reach steady-state flow conditions in these areas, even millions of years in a system with the spatial scale considered here. The other transient processes that were ignored are the potential hydraulic activation of the fault and fracturing. Steady-state analysis is beneficial when identifying areas of pressure relief or recharge/discharge along a fault, and although the 2D model does not capture the full complexity of a real subsurface system, both simplifications are justified to provide insights into the effect on the flow field, including various internal hydraulic structures of the VFZ.

When injecting CO₂ into a storage formation the injection pressure is not constant in time, rather at a constant injection rate the injection pressure will steadily increase with time. However, under the conditions considered here, the injected CO₂ was in a super-critical phase that is much more compressible compared to water. Hence, the injection pressure will typically rise quickly, over the first few years of storage, and then level out and be relatively stable for the remaining time of the injection project. Therefore, to evaluate the leakage potential out of the reservoir, a steady-state scenario was assumed, with a constant injection pressure of a unit pressure in MPa (1 MPa) above the initial in situ pore pressure in the storage formation and the slow increase in pore pressure during early stages of injection was ignored. This injection pressure value is not based on any specific calculation or scenario, as it depends on several factors such as the permeability of the reservoir, the injection rate, the number of injection wells, the length of the injection interval (perforated section of the injection well) and even the stiffness of the reservoir formation. It is therefore beyond the scope of this study to estimate the injection pressure precisely. A fundamental principle in fluid mechanics (Darcy's law) states that the rate of fluid flow through a porous medium is directly proportional to the pressure gradient across the medium, which scales linearly with the injection pressure, thus the leakage rates calculated here will also scale linearly with the injection pressure, relative to the unit pressure applied here.

Although the scenario is related to CO₂ storage, the CO₂ phase is not considered here, only the flow of the brine. This is because the leakage potential under steady-state conditions is being considered, and to consider the CO₂ phase in addition would require transient two-phase fluid-flow simulations to capture the migration and buoyancy of CO₂ in the formations, and the need to consider the CO₂ entry pressure for the sealing units. For any potential leakage, from single- or two-phase fluid flow, a delay in leakage is expected in addition to dynamic leakage rates.

Since steady-state conditions are considered, there is no transient term in the mass conservation equation in equation (6) and it is not necessary to specify properties such as fluid and rock compressibility and porosity. Note that the analysis is also simplified by assuming that the fluid properties, such as density and viscosity, are constant.

Discrete features such as fractures and faults may be modelled using so-called zero-thickness joint or interface elements. Zero-thickness interface elements are introduced as geometrical elements with one less dimension compared to the adjacent continuum. In the 2D model presented here, the fault is thus described by an edge. Here a so-called triple-node interface element is used where two of the nodes describe the dependent variable (pore pressure, p) in the adjacent continuum on each side of the interface, and the third and middle node describes the average hydraulic potential in the fault (expressed by the pore pressure inside the fault, p_f). This type of model description is referred to as a discrete fracture model, by explicitly considering each individual fault and the average hydraulic potential in the faults and the fluid exchange between the faults and the surrounding continuum (e.g. Dietrich et al. 2005).

To evaluate the leakage potential associated with the VFZ, we further define a subsection of the seafloor where the calculated flow rate, Q_F (kg s⁻¹), can be associated with fluid flow up along the fault; this interval is indicated by the orange line in Figure 5.

To solve the boundary-value problem described by equations (6–9) and the boundary conditions shown in Figure 5, the finite-element method was used, discretized with triangular elements, in the commercially available software COMSOL Multiphysics.

To estimate the leakage potential associated with the VFZ, the uncertainty in the hydraulic properties of the geomodel were considered by solving many model realizations with different model parameter combinations. Three lithological groups in the overburden, the Rogaland Group, the Shetland Group and the Cromer Knoll Group, are considered to have higher permeability variations than the best-estimate values in Table 1, and as the hydraulic structure of the VFZ is directly affected by the permeability of the surrounding formations, the permeability values for these three lithological groups were varied in a parametric study (Table 2).

In the parametric study, models were solved where these three chosen lithological groups had five permeabilities ranging from the best-estimate value that represents the low-bound estimate of permeability (Table 1) to a value that is 100 times higher (high-bound estimate). This resulted in 5³ (125) possible model realizations. In addition, in the parametric study, the constants a and b in the correlation were also varied for the fault-width growth coefficients to estimate the fault width (equations 1 and 2). For both coefficients a and b three values were used, corresponding to best-estimate values from Sperrevik et al. (2002), and two values that were ±50% of the best-estimate value, resulting in 3² (nine) possible combinations. In total, 1125 model realizations were solved in the parametric study, the range of these values are summarized in Table 2.

In addition, all of the model realizations based on the permeability in the host-rock lithological groups (125) were solved without including the Vette Fault. In this no-fault model, the VFZ outline only represents a discontinuity, and the effect of the fault is reduced to how the various formations on the hanging wall and footwall are juxtaposed.

To describe the inherent physics in the leakage scenario and the effect of a highly variable fault architecture, both in width and anisotropic permeability values, four model realizations were highlighted out of the 1125 models solved for in the parametric study. These highlighted models were models #251, #375, #751 and #875. Model realizations with a parameter index of 251–375 represent the models with the lowest fault widths but varying permeability in the formations, where in model #251 the overburden formations in the parametric study have the lowest permeability in the study and in model #375 the formations have the highest permeability. Model realizations with a parameter index of 751–875 represent the models with the highest fault widths but varying permeability in the formations, where in model #751 the formations have the lowest permeability in the parametric study and in model #875 the formations have the highest permeability. This is summarized in Table 2.

First, the hydraulic structure of the fault is presented in terms of the fault width and permeability along and across the fault. The fault width varies with throw and the fault-growth coefficients according to equation (2), and the fault-growth coefficients vary with the clay content of the sheared host rock according to equation (1). All combinations of the parameters a and b in Table 2 produces nine different fault-width profiles, these are shown in Figure 6b (coloured lines, top red axes), together with the fault-throw profile (thick black line, bottom black axes). The calculated width of the VFZ was between 2.5 and 7.5% of the fault throw (Fig. 6b), which is in the typical range of reported values (Foxford et al. 1998; Faulkner et al. 2010; Torabi and Berg 2011; Schueller et al. 2013; Alaei and Torabi 2017).

The anisotropic permeability components vary with both the permeability of the sheared host rock and the fault throw; hence, all 1125 permeability profiles will be different, and these are shown in Figure 6c and d. However, the mixing effect on the permeability structure is limited and dominated by the host-rock permeability. When considering the red segment, there are five distinct permeability profiles (Fig. 6c, d) because the hydraulic structure is dominated by the five different permeability values for k in the Rogaland Group However, when considering the magenta segment (highlighted in Fig. 6d), the number of distinct permeability profiles is three because the hydraulic structure is now dominated by the three different fault-width growth coefficients dw (this is due to the large contrast in the permeability of the faulted host rock where the Cromer Knoll Group is mixed with the cap rock Draupne Formation and the reservoir section of the Viking Group).

The fault-width and permeability profile of the four highlighted cases are shown in Figure 7.

A close-up of the VFZ around the interface with the injection reservoir for two of the models in Table 2 is shown in Figure 8 to highlight the permeability structure. Figure 8a shows the low-permeability case, model #751, where the permeability in the three lithological groups have the lowest/best estimate, and Figure 8b shows the high-permeability case, model #875, where the permeability in the three lithological groups have the highest estimate. Note that both models #751 and #875 have the highest fault widths. The grey-scale colour of the background is proportional to the permeability of the various lithologies presented in Figure 3. The permeability structure of the VFZ is illustrated with two sides: the left side corresponds to the along-fault permeability and the right side corresponds to the across-fault permeability. The main difference between the along- and across-fault permeabilities is along the sheared parts from the bottom of the Draupne Formation in the footwall to the bottom of the Viking Group in the hanging wall, highlighting the part of the fault where reservoir sand is mixed into the fault zone, so increasing the permeability along the fault.

The effect of various fault-zone width and permeability profiles for the selected fault model realizations listed in Table 2 and a comparison with the no-fault flow model is presented as Darcy's flux magnitude (apparent velocity) in Figures 9 and 10, and then compared with flow rates out of the model boundaries in Figure 11. A detailed look at the apparent velocities presented in the full profile in Figure 9 shows a larger flow pattern for selected models and outlines the position of the selected section for the VFZ presented in Figure 10. Note that the velocities are highly affected by the adoption of steady-state conditions. In lithological groups and formations with very low permeability (see Table 1), it can take a long time for the flow field to reach steady-state conditions, even millions of years in the spatial scale considered here.

The low-permeability model presented in Figure 10a, c and e shows systematically lower apparent velocity compared to the high-permeability model in Figure 10b, d and f. The thin fault in Figure 10a shows a flow pattern draining the footwall reservoir as well as the overburden formations down the fault and towards the open boundary on the left side of the model (Fig. 11a). The wider fault with the same permeability in Figure 10c drains a similar area but stabilizes at a lower rate, as seen from Figure 11c. When the fault is not included (Fig. 10e), there is a change in the direction of the fluid flow for the footwall overburden draining upwards, although the rate remains low, as seen in Figure 11e. For the high-permeability model, a different pattern is observed for the thin fault, thick fault and no-fault models, showing upwards drainage in the footwall overburden (Fig. 10b, d, f). The corresponding flow rates are similarly shifted to be lower on the left side and higher on the top boundary above the reservoir and fault (Fig. 11b, d, f). An inflection point for the drainage inside the fault zone marks the change from upwards flow in the fault for most of the injection reservoir for the thin fault (Fig. 10b) and downwards flow from somewhere near the bottom of the reservoir section for the thick fault model (Fig. 10d). For the no-fault model, the dominating flow pattern is across the fault and into the Cromer Knoll Group

A summary of the aggregated flow rates from the scenarios in Figure 11 is given in Table 3, providing an overview of the percentage of flow for the various boundaries. In general, it should be noted that the right boundary has a limited flow rate due to the low permeability provided for the basement rock compared to Cromer Knoll Group and the Quaternary succession (Table 1). The left boundary is very distant from the VFZ and detailed geological parameters for this part of the geological section are not the focus here. The most relevant boundary for evaluating the leakage potential is the top surface above the fault. As can be seen from Figure 11, the maximum flow rates on the top surface are found where the Cromer Knoll Group occurs below the Quaternary succession for all the low-permeability cases (Fig. 11a, c, e). From Table 3, we see that the flow rates out of the boundary on top of the fault zone are around one-tenth of the rates for the total boundary (compare the two last columns) for low-permeability cases. As a comparison, the flow rate on top of the fault for high-permeability cases is relatively high, and is around half of the total top boundary. This relatively high flow along the fault for the high-permeability cases is also observed in Figure 11, where the maximum flow rate for the top surface is above the fault for the high-permeability models (Fig. 11b, d).

Flow rates for the models with a fault were compared to the flow rates calculated from the equivalent no-fault model (Fig. 12). The red line in each plot indicates the 1:1 relationship; hence, if a data point (grey dots) is above the red line then the flow rate in the no-fault model is higher than the flow rate for the model with a fault. In Figure 12 the grey-scale colour of the flow-rate points is proportional to the massiveness or area A (equation 3) of the fault, where the light grey-scale colour has large area A and the dark grey-scale colour has low area A. The highlighted cases are indicated with the corresponding colours indicated in the legend in Figure 12a. For the total flow rate (Fig. 12a), the models are spread on both sides of the 1:1 relationship, with more of the no-fault models providing the highest flow rate. Models that include faults dominate for the higher flow rates to the left boundary (Fig. 12b), whereas flow rates to the top surface (Fig. 12c) and the top surface above the fault (Fig. 12d) are systematically higher for the no-fault models. The comparison shows that all flow rates increase with the massiveness (area A) of the fault, and the models with higher permeability in the overburden provide the highest flow rates to the top surface above the fault and top surface in general.

The methodology for deriving the along- and across-fault permeability in this paper was based on algorithms from Bense and Person (2006). Their permeability calculation model incorporates a variety of mechanisms, creating a hydraulic anisotropy in the fault zone such as clay smear, drag of sand and grain reorientation (e.g. Faulkner and Rutter 1998; Faulkner 2004), and vertical segmentation of fault planes, all processes that are found to dominate faulting in the shallow subsurface (<500 m) and which produce an anisotropy in the fault permeability. Applying this similar model to the VFZ can be considered reasonable for the Viking Group and the overburden fault under consideration here as the burial depth of the Sognefjord Formation during faulting is estimated to have been of the order of 0.1–0.5 km (Wu et al. 2021). However, it might be a limitation that the reduction in permeability due to cataclasis and cementation (e.g. Fossen et al. 2007) is not included for the deeper parts of the fault. However, the modelling workflow demonstrates the capacity to utilize well-log data to derive clay content and to utilize the mixing of sand and clay in the fault zone to derive anisotropic permeabilities inside the fault as a function of throw, as shown in Figure 6. More focus on the logging and sampling of overburden formations could reduce the uncertainties in the permeability values for the overburden. The extrapolation of well data into the fault zone could be further supported by more detailed seismic mapping of the overburden.

One challenge with the VFZ permeability model is the calcareous muds and limestones described for the overburden sediments, especially the Cromer Knoll Group and Shetland Group (Table 1). There is limited research available on the key controls on hydraulic properties of carbonate faults in heterogeneous sequences containing carbonates and therefore limited understanding of fault permeability developments in these types of sediments (Michie et al. 2021a). In our model here, this uncertainty has been addressed by increasing the permeability in these formations and generating fault model scenarios with correspondingly higher permeability (Fig. 6c, d). Also, for all the fault permeability scenarios, the minimum and maximum fault permeabilities are limited to the minimum and maximum host-rock permeabilities. This means that increased fault permeability due to fractures that induce permeabilities higher than the porous mixing permeability is assumed an unlikely scenario for the hydrostatic pore pressure gradient and a limited pressure increase in the reservoir (1 MPa) was used in this study.

It should be noted that the permeability anisotropy depends on the contrast in the surrounding host-rock permeability, and we see that this can create a very high anisotropy ratio when very low-permeability material (e.g. from the Cromer Knoll Group and Draupne Formation) is mixed with high-permeability material (e.g. Viking Group) of up to seven orders of magnitude, see Figure 7b–e. In addition to anisotropy due to the mixing of different lithologies, there is also an anisotropy due to heterogeneity in the different lithologies. In the model definition here, this is not considered as the lithologies are defined as homogeneous with a constant effective material permeability representative of each of the lithological formations and groups, and, therefore, there is no anisotropy in the permeability profiles obtained when the same material is juxtaposed; see profiles in Figure 7 and an example in Figure 4. However, when the formations are heterogeneous, the methodology of Bense and Person (2006) will also provide an anisotropy ratio from the permeability in the fault zone; thus, the anisotropy ratio can be a good measure of the heterogeneity of the formations, although this is not illustrated here.

To discuss the implications of a permeable fault in the overburden, the fault permeability model (Fig. 6) is included in a simple 2D flow-simulation model (Fig. 1c). Comparing the apparent velocity, flow field and flow rate for the 1125 model realizations, it can be seen that including a fault with an anisotropic fault permeability significantly alters both the apparent velocity (Figs 9–11) and the flow rates (Figs 11 and 12). The main effect of including the fault is that an along-fault pathway is introduced into the flow model that allows for vertical communication, bypassing sealing units. Variable fault width has the greatest impact down the fault from the footwall reservoir and into the downfaulted reservoir unit in the hanging wall (Fig. 7). This can be attributed to the scaling relationship of the fault width with the throw and that the fault throw maximum is located below the injection reservoir.

In this paper we are comparing the leakage potential to the overburden in a model that includes the VFZ and a no-fault model where the fault outline only represents juxtaposition or a discontinuity in the material properties. The methodology presented here has also previously been compared to the geocellular method (Bjørnarå et al. 2021, 2022). In the geocellular model the geometrical and hydraulic structure of the fault is approximated by an arrangement of rectangular geocells with variable length and width, and the isotropic hydraulic property of each geocell is populated stochastically from a reference value that is based on the shale gouge ratio (SGR) of the fault. Compared to the fault model presented here, the geocellular model did not provide an effective seal and behaved more like the no-fault model. However, the VFZ is surrounded by low-permeability formations and, thus, the fault alone does not alone dictate the leakage potential to the surface. Further, the fault throw drops towards zero at the tips, meaning that there will be very limited mixing and smearing in these parts of the fault, resulting in a permeability structure of the fault that is very similar to the adjacent host rock (red section in Fig. 6c, d). Due to these two circumstances, despite the effective sealing property of the VFZ presented here, when compared to the geocellular model (Bjørnarå et al. 2022) the flow rate to the surface was reduced by less than one order of magnitude, which is similar to what is observed in this study. The main reason for this is the high along-fault permeability in the VFZ that stimulates along-fault flow in the down-fault direction, as shown in Figure 10a and c.

The total flow rate of brine out of the surface in the static calculations here range from 5.4 to 24 kg a⁻¹ m⁻¹ fault length, and for the no-fault model the range is from 5.6 to 93 kg a⁻¹ m⁻¹ fault length (see Fig. 12c). These flow rates are over an area of 49 km (total width of the model), noting that the main contribution to the elevated flow-rate values is from above the reservoir and VFZ (tubes in Fig. 11). Even if the reservoir was filled with CO₂, these rates would not be representative for CO₂ leakage for several reasons. First, the calculated flow rate across the top surface describes the fluid flow of brine due to the pressurization of the reservoir; it does not discriminate from where the fluid comes: for example, the Alpha closure in the Smeaheia fault block (Fig. 1a), where the CO₂ plume would be accumulating. Secondly, the source of brine can in practice be considered infinite, while the source of CO₂ would be limited to how much is stored and will therefore have a leakage profile that increases with time, reaches a steady-state level and then declines as: (1) the availability of the CO₂ is reduced; (2) the driving force, besides buoyancy, for leakage goes down (pore pressure increase due to injection dissipates with time, and eventually injection stops); and (3) the stored CO₂ becomes increasingly stabilized due to trapping mechanisms (residual, solubility and mineralization). Thirdly, the entry pressure for CO₂ in low-permeability and clayey material is high. For the pressurization considered here (1 MPa), data from Sperrevik et al. (2002) suggest that the fault permeability needs to be higher than 10⁻¹⁸ m² to have an entry pressure below a pressurization of 1 MPa. It can be seen from the permeability profiles in Figure 7 that this limits the migration of CO₂ to the fault section, which is mixed with host rock from the reservoir formations, and some high-permeability Cromer Knoll Group model realizations, thus strongly limiting CO₂ migration and the likelihood of CO₂ leakage to the surface. This low potential for leakage is in-line with the risk assessment presented by Wu et al. (2022).

Qualification of structural traps for the injection and storage of CO₂ is important to accelerate the deployment of CO₂ storage projects in saline aquifers. To improve our ability to assess and quantify the sealing properties of faults, we have presented a method to characterize subsurface faults in clastic sedimentary sequences in terms of fault width and fault permeability, based on data that are often available when characterizing CO₂ storage prospects.

The method has been applied to the Horda Platform in the Northern North Sea to describe the hydraulic structure of the Vette Fault Zone (VFZ) and to address the sealing and fluid migration potential near the Alpha structural closure. The example is implemented in a 2D cross-section and a simplified steady-state brine (single-phase) flow simulation with an applied unit (1 MPa) injection pressure in the reservoir section. To analyse the effects of including the fault zone as a discrete hydraulic structure, a total of 1125 fault model realizations were defined to capture variations in the fault width and permeability, and these were compared to no-fault models. Our main findings for the VFZ can be summarized as follows:

• Including the VFZ as a discrete structure in the flow simulation has most impact on the downward flow in the fault and across the fault due to the mixing of sand into the fault zone, which increases the along-fault permeability for this part of the fault. The drainage pathway for the reservoir is down the fault for most realizations. In addition, the overburden drains downwards for many model realizations, whereas for the no-fault models the drainage is upwards towards the top surface.

• In the overburden above the storage reservoir, the fault throw is small and the lithological units are low-permeability clay-rich units; hence, the modelled fault zone is narrow, and the fault-zone permeability is low. The parametric study on the effects of increasing the permeability of the overburden units and fault show increasing flow rates to the top surface and the zone directly above the faults, draining more of the reservoir upwards.

• The modelled brine flow rate towards the surface when the fault model is applied range between 5.4 and 24 kg a⁻¹ m⁻¹ fault length for the model realizations tested in the current study, and for the no-fault model the range is from 5.6 to 93 kg a⁻¹ m⁻¹ fault length. Although these results are similar in magnitude, including the fault reduced the flow rates for all tested model realizations. Note that the calculated rates can be highly biased due to the adoption of a 2D representation and steady-state conditions. A conversion of brine flow rates into CO₂ flow rates is beyond the scope of this work as further understanding of the CO₂ capillary sealing and entry pressure for the fault zone needs to be implemented in the model to extend it for two-phase flow consideration.

• The modelled fault flow rates are highly dependent on the available permeability data and the level of detail for the stratigraphy used for the fault modelling. We experienced best data availability from the reservoir and immediate cap-rock unit, whereas for the overburden and underburden the input data have a higher degree of uncertainty. Improving the input data and further development of the uncertainty quantification might improve the model results.

The authors wish to thank their research partners at the Norwegian Geotechnical Institute (NGI), the Norwegian Research Centre (NORCE), the University of Oslo and the University of Bergen. The authors are also grateful for the many helpful comments the reviewers provided to improve the manuscript.

TIB: conceptualization (lead), data curation (equal), formal analysis (lead), investigation (supporting), methodology (lead), software (lead), validation (lead), visualization (lead), writing – original draft (lead), writing – review & editing (lead); ES: conceptualization (equal), data curation (lead), funding acquisition (lead), investigation (lead), project administration (lead), writing – original draft (equal), writing – review & editing (equal); EAHM: data curation (equal), investigation (supporting), visualization (equal), writing – original draft (supporting), writing – review & editing (equal); SAS: formal analysis (equal), methodology (equal), software (equal), validation (equal), writing – review & editing (supporting).

This work was funded by the Norges Forskningsråd as part of the FRISK project ‘Quantification of fault-related leakage risk (grant No. 294719) and the Norges Forskningsråd Norwegian CCS Centre (grant No. 257579/E20).

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

The datasets generated during and/or analysed during the current study are available from the authors upon reasonable request.