Capturing secondary permeability and fracture connectivity through geological time by fracture back-stripping

doi:10.21203/rs.3.rs-4844693/v1

Download PDF

Article

Capturing secondary permeability and fracture connectivity through geological time by fracture back-stripping

https://doi.org/10.21203/rs.3.rs-4844693/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Fractured rock petrophysical studies rarely use temporal constraints, thus hindering fracture-related permeability (K) and connectivity (C_B) estimation during past geofluid migration. Presenting a new conceptual approach, we use a stochastic method incorporating absolute ages to reconstruct fracture arrays back in geological time. Generating ‘grown’ discrete fracture network (DFN) models, we simulate the hydraulic behaviour of fractured rock volumes from the late-Cretaceous/Palaeocene to the Devonian, via progressive fracture back-stripping. Our study reveals that maximum bulk permeability (K₁) increased through time from 9.47e-14 m² to 4.44e-13 m² (~ 3 orders of magnitude) along with the maximum horizontal permeability orientation (K_Hmax) shifting from NE-SW to NW-SE. Similarly, C_B increased (1.01 to 1.79) with fracture saturation, peaking in the mid-Cretaceous. Validating our method, timing of dated offshore hydrocarbon activity fit our results. Back-stripping of time-constrained fracture sets is therefore a powerful method to investigate and quantify the dynamic evolution of petrophysical properties through geological time.

Earth and environmental sciences/Solid Earth sciences/Tectonics

Earth and environmental sciences/Hydrology

Time constraints, such as the age of fracturing or veining, are seldom utilised in fractured media petrophysical studies due to the general difficulty in deriving such constraints. In turn, the lack of temporal control makes it challenging to establish which fracture-related petrophysical characteristics, such as permeability and fracture connectivity, played a role during past reservoir charging events, and their overall evolution through geological time. Applying time constraints to a well-studied fracture network may resolve which fracture sets, or combination thereof, afforded preferential fluid pathways at specific times. Furthermore, as geofluids may form and migrate episodically^1,2, the transient increase of fracture-related permeability (particularly in rocks with low matrix permeabilities) may be a key control on when and how fluid flow occurs. Appreciating when higher secondary (‘paleo’-) permeability formed in a rock volume, and when significant fracture connectivity was established, can, in turn, provide crucial insights into the dynamics of reservoir charging and geofluid flow histories, significantly aiding geofluid exploration efforts, reservoir production or storage strategies.

Old crystalline basement rocks become progressively saturated with fractures over time when experiencing extended polyphase deformation (Fig. 1), with new fractures forming, or inherited fractures accommodating different deformational episodes³. The addition of each new fracture set to a pre-existing fracture network changes how geofluids flow through a host rock⁴. Accordingly, a fractured rock volume today represents the sum of past fracture, petrophysical, and fluid flow configurations through geological time (Fig. 1). However, the study of selected portions of this history typically requires unravelling a seemingly chaotic array of cross-cutting fractures by relative timing, which is difficult to do⁵. Absolute dating techniques (e.g. Ref^6,7), however, can assist in establishing the formation ages of individual fracture sets, thus permitting us to separate out a specific fracture event in absolute time. Even so, assessing the secondary structural permeability associated with these fractures at the time of their formation remains challenging. Indeed, secondary permeability is both transient and dynamic, being first created by fracturing and then rapidly destroyed by the precipitation of different vein-infilling minerals within brief spans of geological time^4,8,9. As a result, evidence of paleo-permeability in fractured media today is typically only preserved as sealed fractures¹⁰.

To address these challenges, we present a unique case study using a novel stochastic approach to reconstruct complex fracture networks in absolute time, replicating the transient permeability and fracture connectivity evolution of a selected rock volume. Fracture networks are routinely modelled using the discrete fracture network (DFN) approach (e.g. Ref^11–14), which relies on inputs such as fracture size, intensity, shape, orientation, and aperture¹⁵ to generate statistically sound three-dimensional (3D) geometric fracture representations. These ‘static’ (representing the geometric description of a one single geological moment) DFN models are then used to simulate a range of bulk rock petrophysical properties, such as permeability and fracture connectivity¹⁴. To make our stochastic approach time-integrated, this study uses the 'grown' FracMan¹⁶ DFN modelling method^17–19, whereby modelled fractures ‘grow’ in sequential time steps, with dynamic interactions between fractures¹⁸. Crucially, this approach can replicate complex fracture networks involving multiple fracture sets of different relative ages, with younger fractures being affected by older pre-existing fractures¹⁸.

As DFN modelling inputs we used a combination of multiscalar deterministic datasets and temporal constraints provided by Ref²⁰ from Smøla island, mid-Norwegian passive margin (Fig. 2a). Smøla is an ideal setting for this study, it being a crystalline basement exposure that has experienced at least five deformation episodes from the Devonian to Early Palaeogene, which formed complex fracture and fault arrays that have been comprehensively structurally characterised and dated²⁰. Using these inputs, we could therefore generate multiple DFN models representing fracture networks from the present-day back to the Palaeozoic, quantifying for the first time the progressive evolution of key petrophysical and network properties through geological time.

Our modelling shows that bulk permeability and fracture connectivity generally increased through geological time with progressive fracture saturation. To test the approach and validate our results, we compared them against the example of hydrocarbon secondary migration occurrence and timing offshore Norway, with our modelled enhanced permeability and connectivity values coinciding with the onset of most migration activity in the Mid-Cretaceous. The comparison with hydrocarbon activity offshore showed how past reservoir charging could be related to progressive secondary permeability and connectivity development. Furthermore, our determined maximum permeability orientations (anisotropy) highlighted which fracture sets may have represented important fluid pathways at that time (which, in this study, were N-S to NW-SE oriented Norway-Greenland rifting- and break-up-related structures).

Smøla island (Figure 2a) is comprised of crystalline basement rocks (including monzogranite, diorite, gabbro)²⁰, which typically have low matrix permeability (between 10^-18 – 10^-19m²)²¹. To generate our DFN models, we used numerous input parameters, which are described below, and in more detail in the Supplement. Crucial to enabling us to undertake this study, we have made use of multiscalar spatial datasets (lineament trace mapping), high-resolution drill hole and field mapping datasets, and absolute dating of brittle deformation.

Deterministic inputs

Surface lineament traces on Smøla (Figure 2) were mapped at four scales (deterministic trace maps: 1:500, 1:7 500, 1:25 000, and 1:100 000) on orthophoto imagery, digital terrain models (DTM), and aeromagnetic survey data. A total of 9554 fracture traces were picked and the distribution of their trace lengths over four orders of magnitude is best accounted for by a power-law relationship (Figure 3). The power-law distributions, indicating that the number of fractures of relative lengths are consistent across multiple scales (exhibiting a fractal fracture pattern), let us model the fracture networks at our chosen volume dimensions. Lineament trace mapping provided the fracture size distributions and scaling exponents for each of the four identified lineament sets, which are oriented ~ N-S (to NNE-SSW), NE-SW, E-W, and NW-SE (to WNW-ESE). The details of these parameters are available in the Supplement.

In addition to the two-dimensional (2D) multiscalar datasets, downhole televiewer data from four oriented diamond drill holes (by the Geological Survey of Norway-NGU) afforded crucial in-situ geometric information on drill-intersected structures (Figure 4a). These data were sorted by pole to plane clusters, with twelve fracture sets identified over the four drill holes (Figure 4a). Ref²⁰ recognised five deformation episodes (labelled D₁ to D₅) with associated systematic fracture sets on Smøla from field mapping and high-resolution logging of the same four oriented diamond drill holes. Each set is associated with distinct regional geometric trends, synkinematic mineralisations (Figure 4b & c), different deformation episodes, deformation ages derived from K-Ar geochronology of fault gouges²⁰ and U-Pb dating from other studies²². Comparing both the downhole televiewer data to the field studies by Ref²⁰, the following five fracture sets from the televiewer data were thus selected for the DFN modelling (Table 1 & 2) and associated with drill hole linear fracture intensities (P₁₀, refer to Supplement), fracture characteristics (from field work and drill logging), and ages:

The earliest D₁ fractures (mean dip and dip direction: 01/123) are infilled by epidote-prehnite mineralisations (Figure 4c) and formed during a deformation episode dated to >390 Ma²³;
The two D₂ fracture sets (mean dip and dip direction: 00/173 and 74.5/107) are associated with sericite-chlorite-carbonate mineralisations (Figure 4c), and were dated to ~300 Ma and ~204 – 196 Ma²⁰;
D₃ fractures (mean dip and dip direction: 10/075) are decorated by chlorite-hematite mineral coatings (Figure 4c) and were dated to ~128 - 100 Ma²⁰;
The D₄ and D₅ fracture sets were grouped together for simplicity, owing to similar geometries (mean dip and dip direction: 08/242), and are both mostly tensile features, formed post ~75 Ma (70-60 Ma)^20,22 with quartz-calcite-zeolite vein infill (Figure 4c).

DFN modelling

Using the deterministic input values above, we applied the stochastic ‘grown’ DFN modelling approach described by Ref^18,19 to generate five fracture sets within a 500 m³ modelling volume. These fracture sets, which represent the natural fracture trends that accumulated in the Smøla basement during progressive deformation episodes, were used to model the D₁, D₂, D₃, and D₄ and D₅ (labelled as D₄₊₅) fracture sets.

From the drill hole televiewer data, each of the selected fracture sets has a statistical orientation distribution (Table 1), and an associated target volumetric fracture intensity (P₃₂) converted from the relevant downhole linear fracture intensity (P₁₀) as shown on Table 2. Lower-limit P₃₂ values were also calculated for comparison from the areal intensity (P₂₁) values from the 1:500 trace map lineaments (Figure 2d; see Supplement). Each fracture set is attributed a related fracture size (minimum and cut-off size), and a power-law size exponent (α; Figure 3) from the lineament trace maps. Internal fracture properties are also provided for each set, such as aperture values (appraised from extensive drill core vein and fracture measurements) and fracture permeability. Termination relationships (Table 3) between the different lineament traces, which have been correlated to the fracture sets, were also included in the modelling. Further information is provided in the Supplement regarding all the input parameters.

Table 1. Mean orientations of selected fracture sets from the downhole televiewer data. Fracture sets are named by general strike orientation, except for ‘Shallow’ which represents the shallowly dipping D₂ fracture set.

Drill Hole	Fracture Set	Fracture Count	Mean Pole Trend (°)	Mean Pole Plunge (°)	Def. Episode
Drill Hole	Fracture Set	Fracture Count	Mean Pole Trend (°)	Mean Pole Plunge (°)	Def. Episode
BH2	NE-SW	53	122.9	0.66	D₁
BH1	E-W	122	173.2	0.19	D₂
BH1	Shallow	133	106.7	74.53	D₂
BH3	N-S	31	74.65	10.04	D₃
BH4	NW-SE	31	242.04	8.02	D₄₊₅

Four ‘grown’ DFN models were stochastically generated as single realisations (Figure 5), starting with the D₄₊₅ model representing the current fractured rock volume on Smøla (post ~75 Ma). For each temporal iteration, fracture sets of younger sequential ages were progressively removed (back-stripped). The first model (D₄₊₅) thus contains five fractures sets (8 244 732 individual fractures), with an average P₃₂ value of 1.56 m^-1, and contains fractures formed during all the deformation episodes (D₁ to D₅). The second model, with four fracture sets (5 206 985 individual fractures) and an average P₃₂ value of 1.3 m^-1, represents the rock volume post-D₃ deformation (~128 - 100 Ma). The third model, with three fracture sets (5 193 773 individual fractures), and an average P₃₂value of 1.24 m^-1, represents the volume after D₂ (~300 Ma and ~204 – 196 Ma). The fourth model, with one fracture set (10 107 858 individual fractures) and an average 0.67 m^-1 P₃₂ value, represents the volume post-D₁ deformation (>390 Ma).

Table 2. Key model input parameters by fracture sets and ordered by deformation episode. The fracture intensity P₃₂ values are estimated from P10 using the *Ref²⁴ conversion method. Refer to the Supplement for further information on the target P₃₂ values used. Truncation sizes (min=minimum, max=maximum) relates to equivalent radius of modelled fractures. K stands for permeability value of the fracture infill.

Fracture Set Name	P₁₀ (1 m^-1)	P₃₂* (m².m^-3)	Target P₃₂ used				Size (min) (m)	Size (max) (m)	Scaling exponent (α)	Aperture (mm)	K (m²)
Fracture Set Name	P₁₀ (1 m^-1)	P₃₂* (m².m^-3)	D₁	D₂	D₃	D₄₊₅	Size (min) (m)	Size (max) (m)	Scaling exponent (α)	Aperture (mm)	K (m²)
NE-SW	0.48	0.76	0.6	0.3	0.3	0.3	1	20	1.984	0.03	5.1e-12
E-W	1.07	1.68		0.5	0.3	0.3	10	80	1.619	0.04	5.1e-12
Shallow	1.17	1.84		0.4	0.2	0.2	5	60	1.619	0.04	5.1e-12
N-S	0.44	0.66			0.6	0.3	10	80	2.273	0.09	5.1e-12
NW-SE	0.29	0.47				0.5	5	45	2.005	0.1	5.1e-12

Table 3. Probability of termination (%) of a younger mapped lineament trace set against an older ‘passive’ set. Based on geometric similarity to the fracture sets (Table 1), the NE-SW traces are associated with D₁, the E-W traces with D₂, the N-S traces with D₃, and the NW-SE traces with D₄₊₅for the modelling. The ‘Shallow’ fracture set has no corresponding lineament set, as there is no clear 2D map view trace. The adjusted termination probability values used (including for the ‘Shallow’ set), are detailed in the Supplement.

		Younger (lineament trace sets)
		NE-SW	E-W	N-S	NW-SE
Older	NE-SW	19%	21%	15%	27%
	E-W	39%	36%	67%	31%
	N-S	16%	27%	6%	30%
	NW-SE	26%	16%	12%	12%

Within the 500 m³ modelling volumes, we nested a second 200 m³ grid volume with 1000 grid cells (Figure 5). The 200 m³ grid was then used to estimate P₃₂ volumetric intensities and, using the ‘Oda crack tensor’ from Ref²⁵, the equivalent permeability tensor for each grid cell²⁵. From the permeability tensor output (magnitude and orientation), the principal permeability components (K_i), subdivided into maximum (K₁), intermediate (K₂), and minimum (K₃) eigenvectors, could be estimated. For each deformation episode, the three eigenvectors for each of the 1000 grid cells are plotted as poles on stereoplots (Figure 5). An overall mean eigenvector for each of the K₁, K₂, and K₃ poles was then plotted, providing the average 3D orientation of the principal permeability tensor components (Figure 5).

The principal permeability tensor for the D₄₊₅DFN model represents the current permeability anisotropy of the fractured Smøla rock volume; the components exhibit a sub-vertical K₁, a sub-horizontal K₂ oriented NNW, and a K₃ oriented WSW. The average permeability values are 4.44e-13 m² for K₁, 4.01e-13 m² for K₂, and 1.77e-14 m² for K₃ (K₁ > K₂ > K₃). The D₃ tensor components show a K₁ oriented sub-vertically, K₂ oriented ~S, and K₃ oriented ~W to WNW. The permeability values are ~1.6 orders of magnitude lower than those for D₄₊₅, with K₁ at 3.25e-13 m², K₂ at 2.84e-13 m², and K₃ at 1.31e-13 m². The D₂ tensor has a sub-horizontal K₁ oriented ENE-WSW, a sub-vertical intermediate K₂, and a K₃ oriented NNW-SSE. The average permeability values are ~1.7 orders of magnitude lower than D₃, with 2.06e-13 m² for K₁, 1.55e-13 m² for K₂, and 1.10e-13 m² for K₃. For D₁, the earliest permeability tensor, K₁ is sub-vertical to steeply inclined, K₂ is shallowly plunging to sub-horizontal and oriented NNE to NE, and K₃ is oriented NW. The average permeability values are ~2.2 orders of magnitude lower than D₂, with K₁ at 9.47e-14 m², K₂at 9.36e-14 m², and K₃ at 1.55e-14 m² (K₁ ≈ K₂ > K₃). Overall, there is a consistent increase in the principal permeability tensor magnitude through time periods from D₁ to D₄₊₅ (increase of ~2.6 (for K₁), ~2.5 (for K₂), and ~3 orders of magnitude (for K₃)).

The principal permeability tensor component orientations and magnitudes are strongly influenced by the youngest fracture set(s) in each DFN model (most evident in the D₁ model). In the D₃ and D₄₊₅models the sub-vertical intersections between the different steeply dipping fracture sets result in K₁ orientations typically sub-vertical and parallel to these intersections. The D₂ model, however, is a combination of these two situations, with the sub-horizontal intersections of the two D₂ fracture sets and the earlier D₁ fracture set causing K₁ to be sub-horizontal. In general, the principal permeability eigenvector orientations rotate between the models, with K₂ (or K₁ for D₂) changing from towards the NE, to S, then NNW, and K₃ shifting from NW, NNW to W, WSW from D₂ to D₃ (Figure 5).

Topology and fracture connectivity analysis

Both the 3D modelled fractures (Figure 6a) and the 2D mapped (deterministic) lineament traces (Figure 6b) have been quantitatively characterised, following the methodology of Ref²⁶, for network topology and spatial connectivity (Figure 6c). To make comparisons with the final DFN (D₄₊₅) model (representing the same time interval as the deterministic trace maps), and to confirm scale-independent self-similarity between the different deterministic trace maps (mapped at different scales from different parts of Smøla; Figure 2), all the trace maps were analysed using a branch and node counting method²⁷. The proportional counts of branches and nodes was then used to confirm whether the maps form part of the same fracture system, and to provide network topology and connectivity characteristics.

Four modelled trace maps were produced with the 3D fractures represented as 2D lines on horizontal slices through each of the four DFN models (traces being intersection lineations on the planar surface). Three circular representative scan areas (20 m diameters) were selected within each model trace map for network topology analysis (selected scan areas shown on Figure 6a). The branch types include II (isolated), C-I (semi-isolated), and C-C (connected)-type branches. The nodes counted included I (isolated), Y (terminating), and X (cross-cutting), and E (edge)-type nodes²⁶. In total, 2001 nodes and 2530 branches were counted across the maps. Subsequently, each of the four deterministic lineament trace maps were similarly analysed following the same procedure. For each trace map at the different scales, we selected three circular representative scan areas from 4 km (1:100 000) to 40 m (1:500) in diameter. Within each scan area, all the network branches and nodes were similarly counted (selected scan areas shown on Figure 6b). Overall, 2685 branches and 2154 nodes were counted for all the different scales.

The topology ternary diagrams of Figure 6c plot the proportions of nodes and branches from the modelled and deterministic trace maps. Out of the modelled trace maps, the D₁ trace map yielded the highest proportion of I-nodes (~76 %), with D₂ and D₃ having the highest number of X-nodes (~ 41% and ~45 % respectively), and D₄₊₅ the highest proportion of Y-nodes (~32 %). In contrast, the deterministic trace maps show on the node ternary plot that the 1:500 trace map has the highest proportion of Y-nodes (~64 %), while the 1:7500 and the 1:100 000 trace maps exhibit the highest proportion of I-nodes (~41 % and 38 %, respectively). Overall, X-nodes become more abundant with the increasing scale (~9 % to 24 %), suggesting the resolution of the lineament mapping may affect the detection of node type. For all the maps, however, X-nodes are less common than other node types. To compare the two types of trace maps, the final modelled D₄₊₅trace map has fewer Y-nodes, but a similar abundance of C-C branches as the mapped lineament trace maps across the scales. Overall, the modelled trace maps back through time show a trend of increasing number of isolated branches and nodes. Finally, the similar branch and node proportions in the deterministic trace maps suggest they belong to the same fracture system, supporting the use of lineament trace-derived fracture length distributions (Figure 3) in our DFN modelling.

To assess fracture connectivity for the modelled and deterministic maps (Figure 6c), we used the connections per branch (C_B) average value (combined for each scale’s scan areas). The C_B parameter is defined as C_B= (3N_Y + 4N_X) / N_B, where N_y is the number of Y-nodes, N_x is the number of X-nodes, and N_B is the total number of branches expressed as 0.0<C_B<2.0²⁶. For the modelled trace maps, the calculated C_B values are 1.01 for D₁, 1.74 for D₂, 1.82 for D₃, and 1.79 for D₄₊₅. In time, the networks become therefore well-connected from D₂ onwards, with D₄₊₅slightly less well-connected than D₃ owing to the more isolated fractures (now represented by quartz-calcite veins that have shorter fracture equivalent radius sizes). However, the change in connectivity from D₁ to D₄₊₅indicates an evolution in the inter-connectivity of the fracture system through time, with a trend towards higher connectivity with progressive fracture saturation. For comparison, the deterministic lineament traces have calculated C_B values of 1.5 for the 1:100 000 map, 1.75 for 1:25 000, 1.65 for 1:7 500, and 1.8 for the 1:500 map. In general, all the maps have C_B values >1.5, suggesting a good degree of interconnectivity on average. Importantly, the C_B values, and the proportions of branches and nodes for the D₄₊₅modelled trace maps, overlap well and show similarity to the lower-resolution deterministic trace maps (1:100 000 and 1:25 000).

The variety of our deterministic inputs and our novel stochastic approach allowed us to estimate the evolution of secondary permeability and fracture connectivity within a fractured basement rock volume over geological time. The petrophysical conditions allocated to the recognised, characterised, and dated deformation episodes are assumed to have existed at distinct times corresponding to our D₁ to D_4 + 5 models and have evolved with time.

Fractures forming by brittle failure may result from transient fault-valve behaviour, controlled by critical combinations of pore-fluid pressure and differential stress conditions^8,9. Assuming this, fractures can form and progressively evolve during each deformation episode, as also shown by crack-and-seal vein textures in some of the veins on Smøla²⁰. During the episodes, dynamic conditions existed where paleo-apertures, represented by transiently open fractures, repeatedly formed allowing fluid ingress and flow⁹. These open fractures were also repeatedly sealed in brief geological time spans by the precipitation of vein minerals²⁸. Maximum caution is thus necessary when discussing bulk petrophysical characteristics, as they represent not only the effects of transient conditions but may also be the result of multiple events during a single deformation episode that we are now grouping together.

Evolution in secondary permeability and fracture connectivity through time

Remarkably our modelling results make it possible to assess how permeability anisotropy and fracture connectivity have varied through time during the evolution of a tract of the mid-Norwegian passive margin. Insights into how the actual geology of Smøla impacted the results are discussed below. For ease of comparison between the different DFN models, we project our computed maximum and intermediate permeability eigenvectors (K₁ and K₂) of the 3D permeability tensor components onto the 2D horizontal plane as the maximum permeability direction (K_Hmax; Fig. 7).

Model D_4 + 5 (Fig. 5) represents the current fractured rock volume, which is the sum of the accumulated fractures from the Devonian to at least the Late-Cretaceous-Palaeogene (< 75 − 60 Ma). The model predicts a bulk NNW-SSE oriented K_Hmax (Fig. 7), which reflects all the earlier and newly formed fractures. Compared to the earlier models (D₁ to D₃), this model has the highest principal permeability tensor (PPT) component magnitudes (4.44e-13 m² for K₁ to 1.77e-14 m² for K₃), implying that PPT magnitude increases with the number of (interconnected) fracture sets (fracture saturation). During late-Cretaceous to Tertiary times, fluid migration would have therefore focussed along structures oriented ~ NNW-SSE (sub-parallel to our K_Hmax) with steep dips (K₁ sub-vertical), before later mineral precipitation (calcite-quartz-zeolite in our case²⁰) and infill (e.g. Ref⁸). Despite the relative high permeability, the D_4 + 5 model has lower connectivity (C_B = 1.79) compared to the earlier D₃ model (C_B = 1.82) (Fig. 7). Influencing the modelled connectivity, the D_4 + 5 veins are commonly isolated features, rarely cross-cutting older features (contributing few X-nodes), and commonly utilise pre-existing fractures (reactivating older fractures as tensile veins)²⁰, which would indeed result in lower fracture connectivity.

Back-stripping out the D_4 + 5 fractures, model D₃ represents the rock volume during the mid-Cretaceous (~ 128 − 100 Ma). It constrains a N-S to NNW-SSE K_Hmax (Fig. 7), and lower PPT magnitudes (3.25e-13 m² for K₁ to 1.31e-13 m² for K₃). Potential fluid-flow during the mid-Cretaceous would thus have occurred along fractures oriented sub-parallel to this K_Hmax orientation. In terms of connectivity, the fracture termination relationships (Table 3) promoted abundant cross-cutting traces for the D₃ model, which, in turn, resulted in the highest obtained fracture connectivity (C_B = 1.82; Fig. 7). The fracture orientation trends used here shift dramatically from D₂ to D₃ (~ E-W to N-S; Fig. 5), which correlates with a regional stress-field shift in mid-Norway in the mid-Cretaceous³⁰. Whether this is the cause for new cross-cutting fractures rather than reactivating earlier D₂ or D₁ fractures is not however clear. Notably though, the D₃ fractures are indurated and infilled with cataclasite (to gouge; chlorite-hematite mineralised)²⁰. Possible fluid ingress into these fractures may have been thus aided by host rock clasts potential ‘propping’ the features open during fluid flow events³¹. Overall, the D₃ model supports significant fluid flow during the mid-Cretaceous, with relatively high permeability magnitudes (aided by ‘propped’ fractures), and the highest fracture connectivity.

After stripping out the D₃ fracture set, the D₂ model reproduces the rock volume during the Late Carboniferous to Permian and the late Triassic to early Jurassic (~ 300 Ma to ~ 204 − 196 Ma). For that time, modelling yielded a permeability anisotropy K_Hmax oriented NE-SW to ENE-WSW, with a sub-horizontal K₁ orientation (K₁ is sub-vertical in the other models), with the second-lowest PPT magnitudes (2.06e-13 m² for K₁ to 1.10e-13 m² for K₃). The D₂ model represents a prolonged episode, which we modelled by means of three fracture sets (two for D₂ and one for D₁). Both D₂ sets contain phyllosilicate-rich (chlorite-sericite) phyllite shears to shear fractures, with fault gouge²⁰. In fact, the D₂ fractures type and the orientation of the permeability anisotropy indicate complex fluid-flow regimes along either the sub-vertical ~ E-W structures or laterally along shallowly dipping structures (as suggested by the K₁ orientation). Further complicating this picture, authigenic phyllosilicates may have impeded potential fluid-flow (which our modelling ignores), with the D₂ fractures later acting as impermeable barriers³². Irrespectively, the fracture connectivity of this model, being nearly equal to the D₃ and D_4 + 5 models, does demonstrate how a high degree of connectivity (C_B >1.5) was achieved by the late Triassic. The sudden increase in C_B indicates that, at this resolution, fracture networks rapidly become interconnected during deformation. (e.g. Ref⁴).

Finally, back-stripping out the D₂ fracture set, only the earliest initial fractures remain, recreating the basement volume during the Devonian (> 390 Ma). The D₁ model, with early tensile epidote-prehnite veins²⁰ infilling those fractures, provides a NE-SW oriented K_Hmax and the lowest PPT magnitudes (9.47e-14 m² for K₁ to 1.55e-14 m² for K₃) (Fig. 7), ~ 3 orders of magnitude lower than the D_4 + 5 model. Overall, the D₁ model reflects the mostly isolated modelled fractures, and is the least connected of all DFN models, with a C_B value of 1.01. The low C_B value notwithstanding, owing to the low permeability of the host-rock, these initial fractures may still have imparted a significant bulk permeability anisotropy, with fluid flow during the Devonian focused along hydraulically connected fractures sub-parallel to the K_Hmax where available.

Reactivation

As the host rock matrix permeability is intrinsically very low, the fracture sets modelled are the basement’s most permeable features. However, without clear control on the reactivation of pre-existing features, it is challenging to fully constrain the 3D principal permeability tensor. Our modelling addresses this by incorporating pre-existing features with lower volumetric intensity, P₃₂ values (typically decreased by ~ 50% from their initial intensity value) to try and simulate selective reactivation. A significant limitation in our modelling, however, concerns fracture reactivation between different DFN model generations.

Assuming adequate pore fluid pressures, and necessary cohesion and shear strength values ^33,34, fractures optimally oriented 30° to the maximum principal stress (σ₁) and containing the intermediate principal stress (σ₂) should reactivate as shear fractures, while those oriented ~ 0° to σ₁ (and containing σ₂) should reactivate as tensile fractures (with fractures between these two endmember situations reactivating as hybrid fractures)^34,35.

For Smøla, stress fields during the Triassic-Jurassic (D₂) and Early to Late Cretaceous (D₃ and potentially D_4 + 5) indicate E-W and NW-SE to WNW-ESE crustal extension (≈ σ₃) orientations, respectively³⁰. The stress field orientations suggests that, assuming a sub-vertical σ₁ (maximum horizontal stress ≈ σ₂), a portion of favourably oriented fractures from D₁ to D₃ would have been progressively reactivated. Examples of this are late D_4 + 5 quartz-calcite tensile veins in both D₂ and D₃ fractures²⁰. Considering mid-Norway's more recent glacial history, isostatic uplift from ice unloading may have also reactivated many low-angle fractures (sub-vertical σ₃)³⁶.

Quantifying the proportion of reactivating fractures in each DFN model remains unresolved. Dilation and shear tendency studies using assumed paleostress fields may define statistical probabilities of reactivation, providing a basis for intensity estimates, and confirming bulk permeability trends, but this remains presently unanswered in this study.

Novelty of this approach and wider implications

Studying the evolution of various petrophysical characteristics through selective back-stripping of components formed at different times at the scale of a large tract of a passive margin represents an absolute novelty. A similar approach has previously been proposed only in a 3D tomography study of paleo-porosity evolution (e.g. Ref³⁷), and, in an opposite sense, in stepwise forward modelling (restoration) to understand progressive deformation and evolution of permeability (e.g. Ref³⁸). However, those studies lack absolute time constraints, making it difficult to place any results within geologically meaningful time frames. In contrast, our modelling benefits from an absolute time-dimension and an unprecedented well-resolved structural characterisation of the fracture sets on Smøla. These constraints enabled us to precisely control the structure and timing of fracture formation in our DFN reconstructions. Additionally, the ‘grown’ DFN modelling approach includes temporal aspects related to fracture propagation and interaction^18,39,40, simulating pre-existing fractures, during each generation of DFN. Also, considering that sub-seismic-resolution faults and fractures (< 10–30 m in size) in impermeable reservoirs may significantly affect fluid flow⁴¹, using the presented stochastic approach our modelling can model fractures at that resolution. The flexibility in modelling resolution is an advantage when comparing to 3D tomography (fine resolution on the mm to cm-scale) and forward modelling techniques (large resolutions on the km-scale).

To further test and validate our approach, we now use hydrocarbon dynamics from offshore Norway as an example owing to available absolute time-constraints. Hydrocarbons are generally dated using geochronology and biostratigraphic techniques such as Re-Os radiometric dating and biomarkers^42,43. Our DFN models with the associated ages can, therefore, be compared against known reservoir fluid generation or migration events (assuming our 500 m³ modelling volume and the 200 m³ fracture upscaling volume represent a typical reservoir, see Supplement for further details), and broader tectonic events through geological time.

The D₂ DFN model, for example, coincides with the earliest known episode of hydrocarbon generation and secondary migration (Triassic-Jurassic) in the Embla field (Late Triassic fluid migration) and the Halten Terrace, Rås Basin, and Grip High (Jurassic oil generation)^44,45 (Fig. 7). The Triassic and Jurassic witnessed crustal stretching and thinning due to the drifting of Norway and Greenland²⁹ and strike-slip faulting along the Møre-Trøndelag Fault Complex (MTFC), a major NE-SW feature extending from Central Norway into the northern North Sea^20,46. Transtensional kinematics on the MTFC, splay development, and related extension influenced the orientation of the fracture sets that we have used for our modelling, with many being oblique to the strike of the MTFC. Despite hydrocarbon generation in the Halten Terrace during the Late Triassic to Jurassic, major reservoir-scale secondary migration may have been limited owing to the lack of permeable or communicable E-W faults and reactivated NE-SW fractures.

During the Cretaceous to Paleogene (D_4 + 5 and D₃ models), most hydrocarbon generation and secondary migration events initiated, extending into the Neogene and later^44,47–50 (Fig. 7). Importantly, while hydrocarbon generation is linked to the burial and thermal history (maturation) of source rocks⁵¹, our D₃ and D_4 + 5 models do, however, temporally correspond to secondary migration events in the Halten Terrace (mid-Norwegian margin) and Viking Graben (SW Norway). From the Cretaceous onwards, crustal extension of the newly forming passive margin led to major rift basins, basement structural highs (representing ‘buried hill-type’ hydrocarbon traps³¹) and NE-SW to NW-SE striking extensional faults, such as the Klakk fault complex oriented ~ NNW-SSE (similar geometries as our D₃ and D_4 + 5 fracture sets). Therefore, structures oriented ~ N-S to NW-SE forming during the Norway-Greenland rifting and break-up, being parallel to our K_Hmax orientations at these times, may well have been important fluid pathways during the Mid-Cretaceous (D₃) onwards. These conduits would enable communication between basin source rocks and potential trap sites, facilitating fluid percolation and reservoir charging of, for example, fractured basement plays with similar lithologies as Smøla.

As demonstrated by this example, our approach can be adapted to quantify similar geofluid migration histories and is applicable in various contexts where understanding the evolution of bulk rock permeability or fracture interconnectivity is crucial, such as in fracture-hosted groundwater systems, nuclear waste repositories, and CO₂ storage studies.

Our study offers a novel conceptual approach to constrain fracture-related petrophysical properties through time. Using temporally constrained basement-hosted fractures from Smøla, mid-Norwegian passive margin, stochastic modelling reconstructs the evolution of fractured rock volumes as ‘grown’ DFN models, from the present back to the Devonian. To improve realism, we integrated various deterministic remote sensing, geophysical, drill core, and field studies datasets, combined with absolute dating results, to generate stochastic DFN models with petrophysical estimates.

Our well-constrained orientation trends and sorting criteria of different fracture sets, along with absolute ages, allowed us to recreate the fractured rock volume during each deformation episode. By applying a back-stripping methodology to the fracture sets, we deconvoluted a complex present-day fracture network into four time-specific models representing the Late Cretaceous-Palaeocene (and later), the mid-Cretaceous, Late Triassic-Jurassic to Carboniferous-Permian, and finally the Devonian (spanning over ~ 300 Ma). Each of these models are stochastic representations of fracture networks back through time and provide estimates of the bulk principal permeability tensor (component orientations and magnitudes), and fracture connectivity at specific episodes in geological time.

Comparison of our results with an example of absolute-dated geofluid activity, namely petroleum generation and secondary migration offshore Norway, shows that our modelled time-specific petrophysical characteristics fit well with known episodes of hydrocarbon activity and reservoir charging. The ability to estimate permeability and connectivity properties during past fluid migration events represents a significant advancement in understanding fluid flow histories. Understanding how the orientation of permeability anisotropy has changed through geological time can indicate which structures were subparallel to the maximum permeability anisotropy, and thus served as preferential fluid pathways during specific episodes. Additionally, identifying when peak fracture connectivity occurred, which in our example of offshore hydrocarbons activity was during the mid-Cretaceous, provides a basis for understanding why geofluid activity may have commenced at that time (notwithstanding the burial history and source rock maturation in our example).

While incorporating more geological characteristics that affect permeability (such as fracture infill types) and structural reactivation into our models remains beyond the scope of this study (as well as beyond the current resolution possibilities), those inputs could further refine similar approaches in the future. Nonetheless, the framework highlighted in this study remains statistically sound and can be adapted and applied broadly to other geological terranes where absolute time constraints on fractures are available, providing a novel perspective on the history of fluid flow in fractured media.

Acknowledgements

This research has been supported by the ongoing BASE 3 project (BASE – Basement fracturing and weathering on- and offshore Norway) and the Research Council of Norway (Norges Forskningsråd), grant no. 319849. BASE 3 is a joint research project launched and steered by the Geological Survey of Norway (NGU), with the industry partners Equinor ASA, Lundin Norway AS, AkerBP ASA, Spirit Energy Norway AS, and Wintershall Dea Norway AS, and the Norges Forskningsråd (NFR). The authors would like to thank the WSP UK Edinburgh office for the assistance and guidance while completing this study. We thank the BASE 3 project colleagues for discussions and providing constructive input. All the colleagues from the University of Bologna and, in particular, the Deformation, Fluids and Tectonics (DFT) research team, are specially thanked for their continual support and encouragement.

Author contributions

MSH contributed to conceiving the study, collected the data, undertook the statistical analyses, modelling, and results validation, produced the figures and wrote the paper (original draft and editing). MC guided model production and contributed to the editing of the manuscript. JK performed funding acquisition and contributed to editing the manuscript. GV conceptualised the study, obtained funding, acted as PhD project supervisor, and contributed to the writing of the manuscript (original draft and editing).

Competing interests

The authors declare no competing interests.

Materials and correspondence

Correspondence and requests for materials should be addressed to Matthew Sotheron Hodge or Giulio Viola.

Data availability

The key datasets that support the findings of this study are available in the Mendeley Data repository (https://data.mendeley.com/datasets/ctpphnptc5/2) with the identifier DOI: 10.17632/ctpphnptc5.2

Viola, G., Venvik Ganerød, G. & Wahlgren, C.-H. Unraveling 1.5 Ga of brittle deformation history in the Laxemar-Simpevarp area, southeast Sweden: A contribution to the Swedish site investigation study for the disposal of highly radioactive nuclear waste. Tectonics 28, (2009).
Drake, H., Tullborg, E. L. & Page, L. Distinguished multiple events of fracture mineralisation related to far-field orogenic effects in Paleoproterozoic crystalline rocks, Simpevarp area, SE Sweden. Lithos 110, 37–49 (2009).
Viola, G., Kounov, A., Andreoli, M. A. G. & Mattila, J. Brittle tectonic evolution along the western margin of South Africa: More than 500Myr of continued reactivation. Tectonophysics 514–517, 93–114 (2012).
Cox, S. F. Coupling between deformation, fluid pressures, and fluid flow in ore-producing hydrothermal systems at depth in the crust. Economic Geology 100th Anniversary Volume 1905–2005 39 (2005).
Scheiber, T. & Viola, G. Complex Bedrock Fracture Patterns: A Multipronged Approach to Resolve Their Evolution in Space and Time. Tectonics 37, 1030–1062 (2018).
Viola, G. et al. Deconvoluting complex structural histories archived in brittle fault zones. Nat Commun 7, (2016).
Torgersen, E., Viola, G., Zwingmann, H. & Harris, C. Erratum to ‘Structural and temporal evolution of a reactivated brittle-ductile fault - Part II: Timing of fault initiation and reactivation by K-Ar dating of synkinematic illite/muscovite’ [Earth Planet. Sci. Lett. 407 (2014) 221–233]. Earth and Planetary Science Letters vol. 410 211 Preprint at https://doi.org/10.1016/j.epsl.2014.11.013 (2015).
Sibson, R. Crustal stress, faulting and fluid flow. Geological Society, London, Special Publications 78, 69–84 (1994).
Cox, S. F. The application of failure mode diagrams for exploring the roles of fluid pressure and stress states in controlling styles of fracture-controlled permeability enhancement in faults and shear zones. Geofluids 10, 217–233 (2010).
Hosono, H., Takemura, T., Asahina, D. & Otsubo, M. Estimation of paleo-permeability around a seismogenic fault based on permeability tensor from observable geometric information of quartz veins. Earth, Planets and Space 74, 141 (2022).
Long, J. C. S., Gilmour, P. & Witherspoon, P. A. A Model for Steady Fluid Flow in Random Three-Dimensional Networks of Disc-Shaped Fractures. Water Resour Res 21, 1105–1115 (1985).
Dershowitz, W. S. Rock joint systems. (Massachusetts Institute of Technology, 1984).
Dershowitz, W. S. & Einstein, H. H. Characterizing rock joint geometry with joint system models. Rock Mech Rock Eng 21, 21–51 (1988).
Welch, M. J., Lüthje, M. & Oldfield, S. J. Modelling the Evolution of Natural Fracture Networks: Methods for Simulating the Nucleation, Propagation and Interaction of Layer-Bound Fractures. (Springer, 2020). doi:https://doi.org/10.1007/978-3-030-52414-2.
Davy, P. et al. DFN, why, how and what for, concepts, theories and issues. in ARMA International Discrete Fracture Network Engineering Conference D013S000R001 (ARMA, 2018).
WSP UK Limited. FracMan Energy Edition, Version 8.3. Preprint at https://www.wsp.com/en-gb/services/fracman (2024).
Bonneau, F., Henrion, V., Caumon, G., Renard, P. & Sausse, J. A methodology for pseudo-genetic stochastic modeling of discrete fracture networks. Comput Geosci 56, 12–22 (2013).
Libby, S. et al. Grown Discrete Fracture Networks: a new method for generating fractures according to their deformation history. in Proceedings of the 53rd US Rock Mechanics/Geomechanics Symposium 23–26 (2019).
Libby, S. et al. Exploring the impact of fracture interaction on connectivity and flow channelling using grown fracture networks. Quarterly Journal of Engineering Geology and Hydrogeology 57, qjegh2023-010 (2024).
Hodge, M. S. et al. Multiscalar 3D temporal structural characterisation of Smøla island, mid-Norwegian passive margin: an analogue for unravelling the tectonic history of offshore basement highs. Solid Earth 15, 589–615 (2024).
Selvadurai, A. P. S., Boulon, M. J. & Nguyen, T. S. The Permeability of an Intact Granite. Pure Appl Geophys 162, 373–407 (2005).
Hestnes, Å. et al. Constraining the tectonic evolution of rifted continental margins by U-Pb calcite dating. (2023) doi:10.1038/s41598-023-34649-z.
Fossen, H., Khani, H. F., Faleide, J. I., Ksienzyk, A. K. & Dunlap, W. J. Post-Caledonian extension in the West Norway-northern North Sea region: The role of structural inheritance. in Geological Society Special Publication vol. 439 465–486 (Geological Society of London, 2017).
Wang, X., Dove, C., Dunne, J. E., Gutierrez, W. M. & Westman, M. S. Stereological Interpretation of Rock Fracture Traces on Borehole Walls and Other Cylindrical Surfaces. (Virginia Polytechnic Institute and State University, 2005).
Oda, M. Permeability tensor for discontinuous rock masses. Géotechnique 35, 483–495 (1985).
Sanderson, D. J. & Nixon, C. W. The use of topology in fracture network characterization. J Struct Geol 72, 55–66 (2015).
Nyberg, B., Nixon, C. W. & Sanderson, D. J. NetworkGT: A GIS tool for geometric and topological analysis of two-dimensional fracture networks. Geosphere 14, 1618–1634 (2018).
Sibson, R. H. Structural permeability of fluid-driven fault-fracture meshes. J Struct Geol 18, 1031–1042 (1996).
Peron-Pinvidic, G. & Osmundsen, P. T. The Mid Norwegian - NE Greenland conjugate margins: Rifting evolution, margin segmentation, and breakup. Marine and Petroleum Geology vol. 98 162–184 Preprint at https://doi.org/10.1016/j.marpetgeo.2018.08.011 (2018).
Tartaglia, G. et al. Time-constrained multiphase brittle tectonic evolution of the onshore mid-Norwegian passive margin. GSA Bulletin 135, 621–642 (2023).
Holdsworth, R. E. et al. Natural fracture propping and earthquake-induced oil migration in fractured basement reservoirs. Geology 47, 700–704 (2019).
Faulkner, D. R. et al. A review of recent developments concerning the structure, mechanics and fluid flow properties of fault zones. J Struct Geol 32, 1557–1575 (2010).
Nortje, G. S. et al. New faults v. Fault reactivation: Implications for fault cohesion, fluid flow and copper mineralization, Mount Gordon Fault Zone, Mount Isa District, Australia. Geol Soc Spec Publ 359, 287–311 (2011).
Sibson, R. H. A note on fault reactivation. J Struct Geol 7, 751–754 (1985).
Blenkinsop, T. G. Relationships between faults, extension fractures and veins, and stress. J Struct Geol 30, 622–632 (2008).
Stephens, M. B. et al. Review of the deterministic modelling of deformation zones and fracture domains at the site proposed for a spent nuclear fuel repository, Sweden, and consequences of structural anisotropy. Tectonophysics 653, 68–94 (2015).
Macente, A., Dobson, K. J., MacDonald, J., Wadsworth, F. B. & Vasseur, J. The Evolution of Paleo-Porosity in Basalts: Reversing Pore-Filling Mechanisms Using X-Ray Computed Tomography. Transp Porous Media 145, 697–717 (2022).
Roberts, D., Luneburg, C., Jin, D. & Kato, J. Utilizing a kinematic and mechanical modelling workflow to constrain fracture network characteristics: Application to the Teapot Dome, Wyoming, USA. J Struct Geol 159, 104596 (2022).
Maillot, J., Davy, P., Le Goc, R., Darcel, C. & de Dreuzy, J. R. Connectivity, permeability, and channeling in randomly distributed and kinematically defined discrete fracture network models. Water Resour Res 52, 8526–8545 (2016).
Davy, P., Le Goc, R. & Darcel, C. A model of fracture nucleation, growth and arrest, and consequences for fracture density and scaling. J Geophys Res Solid Earth 118, 1393–1407 (2013).
Walsh, J. J., Watterson, J., Heath, A., Gillespie, P. A. & Childs, C. Assessment of the effects of sub-seismic faults on bulk permeabilities of reservoir sequences. Geological Society, London, Special Publications 127, 99–114 (1998).
Selby, D. & Creaser, R. A. Direct Radiometric Dating of Hydrocarbon Deposits Using Rhenium-Osmium Isotopes. Science (1979) 308, 1293–1295 (2005).
Bhullar, A. G., Karlsen, D. A., Backer-Owe, K., Seland, R. T. & Le Tran, K. Dating reservoir filling — a case history from the North Sea. Mar Pet Geol 16, 581–603 (1999).
Ohm, S. E., Karlsen, D. A., Phan, N. T., Strand, T. & Iversen, G. Present Jurassic petroleum charge facing Paleozoic biodegraded oil: Geochemical challenges and potential upsides, Embla field, North Sea. Am Assoc Pet Geol Bull 96, 1523–1552 (2012).
Matapour, Z. & Karlsen, D. A. Ages of Norwegian oils and bitumen based on age-specific biomarkers. Petroleum Geoscience 24, 92–101 (2018).
Watts, L. M., Holdsworth, R. E., Roberts, D., Sleight, J. M. & Walker, R. J. Structural evolution of the reactivated Møre-Trøndelag Fault Complex, Fosen Peninsula, Norway. J Geol Soc London (2023) doi:10.1144/jgs2022-139.
Goff, J. C. Hydrocarbon generation and migration from Jurassic source rocks in the E Shetland Basin and Viking Graben of the northern North Sea. J Geol Soc London 140, 445–474 (1983).
Georgiev, S. V et al. Comprehensive evolution of a petroleum system in absolute time: The example of Brynhild, Norwegian North Sea. Chem Geol 522, 260–282 (2019).
Khorasani, G. K. Factors Controlling Source Rock Potential of the Mesozoic Coal-bearing Strata from Offshore Central Norway: Application to Petroleum Exploration. Bulletin of Canadian Petroleum Geology 37, 417–427 (1989).
Ali, A. Basin Modelling: HC Generation Modelling of the Åre, Melke and Spekk Formations, Haltenbanken Area. (University of Oslo, 2012).
Palumbo, F., Main, I. G. & Zito, G. The thermal evolution of sedimentary basins and its effect on the maturation of hydrocarbons. Geophys J Int 139, 248–260 (1999).

There is NO Competing Interest.

HodgeetalSupplementarymaterial.pdf

Download PDF

Version 1

posted

You are reading this latest preprint version

Capturing secondary permeability and fracture connectivity through geological time by fracture back-stripping

Status:

Version 1

Abstract

Figures

Introduction

Results

Deterministic inputs

DFN modelling

Topology and fracture connectivity analysis

Discussion

Evolution in secondary permeability and fracture connectivity through time

Reactivation

Novelty of this approach and wider implications

Conclusions

Declarations

Acknowledgements

Author contributions

Competing interests

Materials and correspondence

Data availability

References

Additional Declarations

Supplementary Files

Status:

Version 1