Assessing fault zones as fluid-migration pathways requires the characterization of permeability both across and along faults, as well as the adjacent volume. The hydraulic properties of the Vette Fault Zone, North Sea, are described by modelling the mixing of host-rock lithologies into the fault zone, and the fault width is derived from empirical relationships as a function of throw and clay content. To better understand the sensitivity related to the uncertainties in overburden lithologies and fault-width correlations, a parametric study with 1125 model realizations were solved in a 2D steady-state, single-phase, subsurface flow model. The fault zone, included as a discrete permeable structure, significantly alters the flow field compared to a model that only considers lithological juxtaposition. The most prominent hydraulic communication in the Vette Fault Zone is downwards from the storage reservoir where sand is mixed into the fault zone. Increasing the host-rock permeability in the overburden also increases the fault permeability and shifts the inflection point for down-fault flow, causing the pressurized reservoir to drain towards the overburden and the top surface. For CO2 storage application, the models highlighted the potential for downward communication along the fault for brine, and the CO2 capillary sealing towards the overburden.

Thematic collection: This article is part of the Fault and top seals 2022 collection available at:

Implementation of large-scale CO2 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 CO2 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 CO2 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 CO2 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 CO2 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 CO2 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 CO2 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 (Vshale) (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 (Vshale) 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 Vshale 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.

Fault flow model

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 CO2 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 Vclay (%) and permeability of the host rock, as described below.

Fault width typically increases with throw (Torabi and Berg 2011; Alaei and Torabi 2017). However, variations of width with throw also depend on the clay content of the host rock. Based on data from Sperrevik et al. (2002), Bense and Person (2006) derived a correlation between the fault-width growth coefficients per metre of fault throw, dw (m/m), and Vclay (%):
where the parameters a and b are fitting parameters. The data from Sperrevik et al. (2002) give a = 0.07 and b = −0.02, so for a fault section with Vclay=0, the width increases with the throw by dw = a = 0.07; and for a clay-rich fault section with Vclay=100%;, the width increases with the throw by dw = 0.0095. The total width, W (m), of the fault is the sum of the fault-width growth contribution from the downthrown and upthrown host rock:
where subscripts d and u refer to the downthrown side and the upthrown side, respectively. Dd (m) is the fault throw on the downthrown side relative to the upthrown side of the fault, and Du (m) is the fault throw on the upthrown side relative to the downthrown side of the fault. In the case where one side of the fault is considered constant/stationary relative to the other side (e.g. in the case of a listric fault zone), the hanging wall is the downthrown side such that Du is zero. It should be noted that the fault thickness models might overestimate fault thickness in deeper parts of VFZ due to the listric fault growth, whereas in the overburden the models should work well.
A useful metric when evaluating the effect of faults is the massiveness, or cross-sectional area A (m2), of the fault. The massiveness can be found by integrating the fault width along the length, l (m), of the fault:
Bense and Person (2006) also included algorithms to calculate the anisotropic permeability components. Due to the shearing of the host rock, the hydraulic structure in the fault becomes layered, creating an anisotropy. The reason for the anisotropy is nuanced and can be due to several factors such as stress orientation, porosity, cataclasis, shear structures and veining (Farrell et al. 2014). The anisotropic components of the permeability can be derived using averaging techniques for the estimation of the permeability of a layered porous media. The highest-permeability component is parallel to the fault, kf,max=kf, (can be either parallel to or perpendicular to the slip direction), and the minimum permeability component is generally perpendicular to the fault plane (Farrell et al. 2014), kf,min=kf,. The typical averaging technique for the permeability perpendicular to the fault, kf, (m2), is calculated by the harmonic mean:
where dwd and dwu (no units) are the fault-width growth coefficients from the downthrown and upthrown host rock, respectively; and kd (m2) and ku (m2) are the permeability of the host rock on the downthrown and upthrown side of the fault, respectively. The permeability profiles of the host rock are shown in Figure 3.
The typical averaging technique for the permeability parallel to the fault, kf, (m2), is calculated using the arithmetic mean:
It can be seen from the expressions for permeability in equations (4) and (5) that the values are weighted with the fault-width growth coefficients that depend on Vclay. The expressions in equation (4) for the across-fault direction are similar to the expressions used to estimate the transmissibility and transmissibility multipliers (Manzocchi et al. 1999) for faults to evaluate their sealing properties.

Fault characteristics example

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 Du = 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.

Regional flow model

Model assumptions

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 CO2 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 CO2 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 CO2 storage, the CO2 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 CO2 phase in addition would require transient two-phase fluid-flow simulations to capture the migration and buoyancy of CO2 in the formations, and the need to consider the CO2 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.

Governing flow equations

The pressurization of the storage reservoir will result in an outward fluid flow that can be described by the mass conservation equation for the fluid (equation 6, here at steady state):
where ρ = 1000 kg m−3 is the density of the fluid and the specific discharge, q (m s−1), can be described by Darcy's law (equation 7):
where μ = 1 mPa s is the viscosity of the fluid, k (m2) is the permeability and p (Pa) is the pore pressure.

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, pf). 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).

Fluid flow along (longitudinal) and across (transversal) a fault is expressed by a dimensionally reduced formulation of the mass conservation equation for the fluid where both equations (6) and (7) are integrated across the width (m) of the fault. When the properties in the fault are constant, the upscaled governing equation becomes:
where the subscript t indicates the tangential component, kf, (m2) is the permeability parallel to the fault and pf (Pa) is the pore pressure inside the fault. The tangential derivative of p, tpf, is the tangential projection of a gradient on an edge (or surface in 3D), tpf=(InnT)pf, where I is the unit tensor and n is the unit normal vector of the surface. The source term qΔ (kg m−2 s−1) in equation (8) describes the fluid exchange across the upscaled fault. If the fault is surrounded by impermeable rocks or in the case of continuity in pore pressure across the fault, then qΔ=0. However, when there is a discontinuity, a jump in pore pressure, Δp, across the fault occurs, and can be described using a thin-layer approximation (here across an interface with thickness W, which is also the width of the fault):
where kf, (m2) is the permeability perpendicular to the fault.
The boundary and point conditions are shown in Figure 5. The scenario that is solved for is a pressurized reservoir and this is accomplished by constraining a point inside the reservoir (yellow domain in Fig. 5) to a fixed elevated pore pressure or injection pressure, pinj = 1 MPa (red point in Fig. 5), relative to the initial pore pressure. The bottom boundary (thick black line in Fig. 5) in the model has no-flow boundary conditions, and the lateral sides and top surface (seafloor) are described as open boundaries with constant hydrostatic pore pressure equivalent to the initial pore pressure (red, green and blue boundaries in Fig. 5, respectively). With these constraints in the model the fluid will flow from the reservoir (with the elevated pore pressure), and the flow rate out of the reservoir, QR (kg s−1), is in balance with the flow rate out of the open boundaries to the sides, QL (kg s−1) and QR (kg s−1), and to the top (seafloor), QS (kg s−1):

To evaluate the leakage potential associated with the VFZ, we further define a subsection of the seafloor where the calculated flow rate, QF (kg s−1), 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 (69) 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.

Parametric study for the Vette Fault Zone

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 53 (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 32 (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.

Width and permeability profiles of Vette Fault Zone

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.

Flow pattern in and around the Vette Fault Zone

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.

Main learnings and model limitations

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 911) 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.

Application for the Vette Fault Zone and CO2 leakage

The total flow rate of brine out of the surface in the static calculations here range from 5.4 to 24 kg a−1 m−1 fault length, and for the no-fault model the range is from 5.6 to 93 kg a−1 m−1 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 CO2, these rates would not be representative for CO2 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 CO2 plume would be accumulating. Secondly, the source of brine can in practice be considered infinite, while the source of CO2 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 CO2 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 CO2 becomes increasingly stabilized due to trapping mechanisms (residual, solubility and mineralization). Thirdly, the entry pressure for CO2 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−18 m2 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 CO2 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 CO2 migration and the likelihood of CO2 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 CO2 is important to accelerate the deployment of CO2 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 CO2 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−1 m−1 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−1 m−1 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 CO2 flow rates is beyond the scope of this work as further understanding of the CO2 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.

Using constant fault characteristics for the width and permeability is only applicable over a short fault segment. When evaluating basin-scale effects, such as leakage to the surface from a storage reservoir, the overburden plays an important role and a more detailed characterization of the fault, as demonstrated here, is required. This is motivated by the results shown in this study that the presence of a fault significantly alters the flow field around the fault. Although the fault does not have a large effect on the flow rates towards the surface in this study, this may not always be the case, and, thus, ignoring or overly simplifying the fault is a significant alteration of the geomodel.

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.

This is an Open Access article distributed under the terms of the Creative Commons Attribution 4.0 License (