## Abstract

The chalk reservoir of the Gorm Field, southern North Sea is dome-shaped and faulted owing to a combination of salt diapirism and regional east–west extension. Fractures developed in the structure considerably enhance permeability. The dataset discussed here records fractures in horizontal wells from more than 10 km of image logs and provides a special opportunity to test theoretical models of fracture development with quantitative observations. In an effort to forecast fracture density and fracture orientation, we have estimated the strains in the structure using an elastic dislocation model that incorporates mechanical boundaries in the form of the tectono-stratigraphic interface with salt and tectonic faults. More than 50% of the angular differences between poles to the planes of simulated and observed fractures are less than 30°; 75% are less than 45°. Relative strain magnitude appears to be a useful indicator of fracture density. At the field scale, small strain magnitudes correspond with small non-zero fracture densities and relatively large strain magnitudes correspond with high fracture densities.

The Gorm Field is situated in the southern part of the Danish Central Graben (Fig. 1). It is operated by Maersk Oil. With hydrocarbon production commencing in 1981, to date it has produced over 372.5 million barrels of oil and 580.5 billion scf of gas. The Gorm platform is central to the infrastructure of the Danish Central Graben pipeline system and with estimated remaining reserves of around 25 million stock tank barrels (stb) at current production of 13 000 stb/day, accuracy of infill wells and optimizing interventions is vital for continued operation of the field. The Danian and Maastrichtian chalk reservoirs record small matrix permeabilities less than 8 mD with faults and fractures contributing significantly to enhance permeability (Megson 1992). The Gorm structure, located in the Chalk above a salt diapir, is a faulted dome where at the depth of present hydrocarbon exploitation it has a diameter of approximately 4 km. The fault pattern is dominated by a NNE-striking graben transecting the crest of the dome, bounded by a pair of normal faults (maximum throw of 170 m) that dip steeply (up to 80°) at the very top of the structure (approximately 1900 m), shallowing to 50° on the flanks and becoming listric with depth. The remaining (several tens) faults have dips in the range 40–60°. The structure has been penetrated by at least 111 wells, of which 13 have fracture logs available in the economic interval and, from these, it is possible to determine densities and orientations of faults and fractures that are too small to image seismically.

The ability to make reliable quantitative or semi-quantitative predictions about the spatial and statistical properties of natural factures is becoming increasingly important to the hydrocarbon industry. From a statistical point of view the scaling relationships of fault-size frequencies has been extensively documented since the work of Scholz & Cowie (1990). Providing there are observations at several orders of magnitude of size, these relationships help with the question of the numbers of fractures, but they are unable to address the related problem of location and orientation. Combining information from size, orientation and intensity distributions by way of a Monte Carlo scheme leads to discrete fracture networks (e.g. Dershowitz & Einstein 1988) of connected fracture surfaces. This powerful approach yields realizations that are faithful to the size, orientation and intensity statistics but it requires extensive detailed observations in order to be predictive. An alternative approach that we discuss in this paper is based on the hypothesis that the numbers, size and orientations of small-scale fractures might be related to the amount of deformation of the host rock. It is not usually possible to measure three-dimensional strains in the subsurface. Instead, we have to formulate a mathematical model of the deformation in the volume, generating, for example, deformed models of horizon surfaces and/or other reference data. If these modelled surfaces are sufficiently similar to the observed geometry of the real reference surfaces, we make the conceptual step that the modelled strains both at and between those surfaces are a plausible, first-order approximation of the actual rock strain.

The principal objective of this paper is to determine the feasibility of using routinely available seismic interpretation with a continuum model of the subsurface to predict, in a qualified sense, the relative intensity and orientation of subseismic fractures with only structural interpretation as the input. Here we consider the total deformation associated with emplacement of the salt diapir and the coeval/subsequent faulting. Specifically, we consider the case where the deformation in the volume of interest is driven by displacements at surfaces that form internal mechanical boundaries such as (a) fault planes and (b) the top of the pre-mobilized salt unit within a homogeneous elastic solid. We refer to this as elastic dislocation (ED) modelling (see also Dee *et al.* 2007).

The use of ED theory has been documented in both the structural geology and petroleum geology literature, but its routine application to reservoir fracture prediction has received little attention in peer-reviewed journals. However, Bourne & Willemse (2001) have shown that the displacement–discontinuity boundary element method (Crouch & Starfield 1990) provides an explanation of fault-related opening-mode fractures at Nash Point, UK. Using the same underlying theory, Maerten *et al.* (2002) demonstrated in the Oseberg Sør field of the northern North Sea that modelled stress perturbations associated with major seismically mapped faults are effective predictors of the orientation and intensity of minor seismically mapped faults. Dee *et al.* (2007) have described similar results for a contractional structure in Columbia and extensional structures in the North Sea and Algeria. Applying a similar methodology but based on a kinematic model, Lohr *et al.* (2008) explore the relationship between strain and fracturing in the Lower Permian sandstones of the NW German basin. Unlike the ED method, their three-dimensional deformation pattern is derived from structural restoration.

## Geological background

### The Gorm structure

The prominent dome structure of Gorm has been caused primarily by mobility of underlying Zechstein salt. Regional halokinesis started to develop in the Mid to Late Jurassic, coinciding with extension associated with the formation of the Central Graben (Megson 1992).

The base of the Gorm diapir is interpreted to be seated at the top of the footwall (*c.* 5500 m depth) of a major Triassic fault with a throw of approximately 2000 m. The diapir is approximately elliptical in map projection with a slightly sinuous north–south trend that also reflects the shape of the Triassic fault scarp. The diapir penetrates rocks of Jurassic age with the present day top of salt at an approximate depth of 2750 m (see Fig. 2), 250 m beneath the Base Cretaceous and 750 m beneath the top of the Top Chalk.

The Chalk is significantly thinner above the diapir than on its flanks. We interpret a component of this thinning to be due to strains associated with the mechanical emplacement of the salt, but given the time-frame outlined by Megson (1992), it would seem reasonable that some of the thinning is also sedimentary, that is, the Chalk was being deposited on an evolving structural high.

The domed structure of the Chalk forms a four-way dip closure beneath a cap rock of Lower Tertiary shales and above a thin Lower Cretaceous sequence, creating an effective trap for oil originating from Jurassic source rocks. The reservoir of Gorm constitutes Danian and Maastrichtian chalk layers at a shallowest depth of approximately 1900 m with a thickness of over 500 m combined, 300 m of which lies within the economic zone.

The dome is transected by a broadly NNE–SSW fault system. These faults have upper terminations beneath rocks of Early Tertiary age. The largest fault in this system (F1 in Fig. 2), having a maximum throw of 170 m at the Top Chalk, divides the field into the western A-block and the eastern B-block (Fig. 2).

### Chalk reservoir

Matrix porosity values measured in wells across the field range from 15 to 45%, with matrix permeability ranging from 0.5 to 8 mD. Faults and fractures contribute significantly to enhance the permeability. The maximum oil column height is 71 m with negligible gas columns. The reservoir overlies a thin Lower Cretaceous to Jurassic age succession (and the core of Zechstein salt).

### Seismic data

Structural interpretation (see below) is derived from two seismic data sets of vintage 1988 and 2005. Both surveys were time-migrated then depth converted using a three-dimensional velocity model. We estimate that the smallest faulted offsets imaged in the seismic correspond to approximately 20 m.

### Well data

Later in the text we describe how the principal tests of the ED model depend on direct comparisons with observed fractures in wells. This section describes the subseismic fractures recorded along the well paths.

Including side tracks, 111 wells had been drilled in the Gorm Field by the end of 2010, with borehole imaging (including FMS/FMI/STAR) available for 13 deviated wells (*c.* 17 000 m in total). Within the reservoir interval sedimentological and structural surfaces are identified from the image logs using dip analysis.

Despite the ostensibly wide coverage of the image logs, their interpretation has some important limitations. First, there is an operational sample bias in as much as most of the 13 wells have approximately radial paths with respect to the dome and do not penetrate the top of the reservoir in the crestal graben; they are concentrated in the relatively uplifted flanks. Second, image log data themselves suffers from orientation bias; it is commonly accepted that fractures oriented parallel or subparallel to the well track will be undersampled. Image artefacts can also mask real geological features that are perpendicular to the well bore, resulting in poor recognition of hole-normal fractures.

Collectively, these limitations militate against reliable observations of fractures that are radial or tangential to the structure, so we have to acknowledge that measured densities and orientation patterns are compromised to some degree.

#### Fracture classification

More than 16 000 faults/fractures with variable expression have been identified from the image logs. The vast majority of these (more than 15 500) are small. They do not clearly transect the entire well bore and they are classified as joints, minor natural fractures and drilling-induced fractures. In this study we have limited the analysis to those that are visible in all four quadrants of the image log. Some but not all of these also show lithological variation across the fracture surface, indicating shear displacement. Some have finer-grained material developed in their cores, which again we take as an indication of shear. This subset of major fractures comprises a total of 418 subseismic faults that, when counted along the well path at a scale of 100 m, range in density (*D*_{f}) from 0.0 to 0.22 m^{−1}.

Maerten *et al.* (2002) suggest that the perturbations in stress associated with faults in an elastic model are dependent on the size of the fault boundaries and magnitude of offset. In general, perturbations associated with a set of faults with hundreds of metres of displacement predict orientation and intensity of faults with displacements on the scale of tens of metres. In this work our smallest seismically imaged faults have displacements of approximately 20 m, so we anticipate that we might forecast structures with offsets up to a few metres. For this reason we limit the comparison and analysis to the major fractures picked in the image logs.

#### Structural domains

Since the wells are strongly deviated, each one may cross one or more of the framework faults (see below). In an effort to isolate the fracture density characteristics to specific structural domains, we have broken down the original 13 wells into 25 subsets (Fig. 3). The boundary of a subset or domain corresponds to a framework fault, thus we avoid mixing information from adjacent footwalls and hanging walls. Although this has no implications for comparison of fracture orientations, it is an important consideration for fracture density.

## Methodology

### Structural interpretation

For the purposes of this work the field was interpreted using both the original 3D seismic survey from 1988 and the more recent 2005 survey. Dependable horizon picks come largely from the top of the Chalk sequence, namely the Top Chalk and the top of the M2 zone of the Maastrichtian. Where possible we have included picks from deeper horizons in order to constrain the displacement patterns since these provide further constraints for the ED modelling. Faults were initially picked as segments on in-lines and cross-lines and then correlated to create fault surfaces. Each surface is represented as a triangulated grid derived from the fault segments. The gridding process (grid dimension of 50 m) produces smooth fault planes. Where the interpretation requires two faults to be linked, the branch line is a fault segment common to both surfaces. The structural model incorporates all of the faults that we observe in the seismic data, picked over their entire imaged depth range. However, we suspect that the imaged range is incomplete since most of the faults are not visible below the Base Cretaceous and all of the faults are truncated a short distance (less than 100 m) above the Top Chalk. The latter implies that the Top Chalk was close to the Earth's surface at the time of faulting.

We distinguish (arbitrarily) two groups of faults based on how significantly they impact the mapping of the structure. Faults that are mapped as deep as the base of the Chalk fall into our ‘framework’ set; faults that cannot be traced so deep form our ‘infill’ set. On the east of Gorm (B-block) the framework faults dip predominantly towards the west. The steepest dips, directly above the top of salt (on fault F1), reach 80°. Moving north and south away from the diapir, the F1 fault-dips decrease to between 30° and 60°. F1 and a number of other B-block Framework faults have dips that decrease with depth to values as low as 20°. The remainder of the Framework faults in B-block are relatively less listric with dips between 30° and 60° (in the same fault plane). Conversely, most of the A-block framework faults dip towards the east. Their dips range between 40 and 60° and there is no suggestion of dips decreasing with depth. The infill fault set also has dips in the range 40–70°. The surface dips and orientations of the framework faults are summarized as colour-coded surfaces in Figure 4a and on the equal-area stereographic projection in Figure 4c.

Examination of the fault branch lines shows a relationship between orientation and location above the salt diapir (Fig. 4c). On the topmost salt and northern ridge, the branch lines plunge SE and to the south they plunge north to NW, broadly parallel with the axis of the ridge. On the east and west flanks, the branch lines plunge inwards. This convergent pattern implies that salt emplacement exerted an important secondary influence on the geometry of the extensional faulting. In the stereographic plot of fault orientation (Fig. 4c) we also show the eigenvectors of the poles to the fault surfaces. The third eigenvector (or anti-cluster) clearly plunges at a low angle (<15°) to the north, which is consistent with the normal fault system being primarily Andersonian in nature. The azimuth of the first eigenvector (*c.* 100°) also turns out to be the extension direction in our preferred model.

While a structural model derived from seismic interpretation is always non-unique, we have been careful to ensure wherever possible that our combined fault and horizon interpretation is qualified by structural analysis. It is always important to validate the structural model arising from seismic interpretation but doubly so when that model is to be used as the basis for derivative modelling, as is the case here. We have performed validation using the displacement analysis procedure of Barnett *et al.* (1987) and extended by Freeman *et al.* (2010) to limit interpretation errors by considering the inferred wall rock strain. Figure 4b shows the framework set colour-coded in throw. Note that the colour range is specific to each fault since it is normalized by the maximum displacement on that particular fault. This permits the shape of the throw patterns to be viewed more clearly. There are few examples of the ideal elliptical throw pattern of Barnett *et al.* (1987). This is easily justified by the fact that the maximum displacements are slightly below the imaged tops of the faults that mark the Earth's surface at the time of faulting. Remaining details of the patterns are largely smooth and continuous, but they display abrupt changes where they cross a branch line. As a final step in the validation we have plotted mapped dimensions of the faults against their maximum displacement. The interpreted faults from this paper compare favourably with a database of published length–displacement data (e.g. Schultz *et al.* 2008).

Part of the structural analysis involves generating the intersection polygons that form the horizon boundaries at each of the fault surfaces and form the primary information for calculating offset. By aggregating the offsets (heaves) of all the faults at a particular horizon and in a direction perpendicular to average strike, we are able to estimate the total extension due to faulting. Dividing the extension by the length over which each transect has influence gives an estimate of the extensional strain. Figure 5 shows that for the Top M2 horizon the interpreted structural model indicates an extensional strain of between 0.03 and 0.04. Repeating the exercise at deeper horizons produces substantially similar results. We acknowledge, however, that strain estimates from seismic-scale faults may underestimate the real extension by a factor of 2 (Walsh *et al.* 1996). Because the fault system extends north and south beyond the immediate influence of the salt dome, the pattern is consistent with a response to real tectonic extension rather than a manifestation of outer-arc extension.

### Elastic dislocation modelling

In this section we describe the ED methodology used to estimate the strain in reservoir units of the Gorm Field. Note that these estimates are first order and are derived from governing equations defined for a homogeneous half-space with isotropic, linear elastic properties. Note further that the common justification for the applicability of elastic models to geological situations involving finite but small strains (as opposed to engineering strains) is simply that strains accrued by elastic processes become permanent as stress relaxes by some (unspecified) process over long periods of geological time. Unless stated otherwise, we use the term ‘ED model’ to refer to a forward model.

#### Essential considerations

Expanding on the description given by Dee *et al.* (2007), the basic unit driving deformation in an ED model is a planar element. Solutions for the displacement and strains in the volume associated with a uniform planar dislocation element in a homogeneous elastic medium are given by Comninou & Dunders 1975 (angular), Jeyakumaran *et al.* 1992 (triangular) and Okada 1992 (rectangular). Given a defined displacement (shear or dilation) on a single element, the elastic solutions provide both the displacement (vector) and strain (tensor) at any desired point in the half-space. Given many elements and an arbitrary point in the volume, the sum of the contributions from each element at that point provides the resultant displacement vector and strain tensor in accordance with the rule of superposition for linear elastic solids.

In this paper, all faults and, in addition, the geological interface between the salt and the overburden are considered to be mechanical boundaries. The fault surfaces have each been decomposed into a set of rectangular elements with horizontal upper edges. For each fault the observed dip separation over the region of an element is mapped on to the element and itself, decomposed into horizontal- and dip-shear components (in most cases the horizontal shear components are zero). Because the area of the observed fault surface and the summed area of the elements may differ slightly, the displacement is scaled in such a way that the total moment of the interpreted fault and the total moment for all the elements on a single fault is the same. In keeping with the conclusion that faults are approximately Andersonian and normal in sense, the slip vector on each element of a fault has been chosen to correspond to the average dip vector for the same real fault plane.

The interface between the salt and the overburden is represented by a set of triangular elements. At each of these elements we have estimated the amount of vertical motion relative to an assumed, pre-mobilized top of salt. There are two end-member cases for the initial position of the top of salt: (a) approximately uniform thickness with a volume (in the region of interest) equal to that in the diapir; and (b) zero thickness with all salt in the diapir flowing from outside the region of interest. In the first case the relative motion may be uplift or subsidence while in the second case it will be uplift only. The measure of uplift or subsidence is used as the dilation displacement on the ED element. Collectively the elements defining the salt/overburden interface effectively model the shape change of the interface from the pre-mobilized boundary to the present day configuration. The displacements vectors and strains tensors at a pre-specified point are determined in exactly the same way as described above for faults. Total displacement and strain at a point are then the sum of the contribution from the salt and the faults.

#### Observation points

We refer the term ‘observation point’ to any point in the volume where we require an ED solution. In this work we seek solutions at points on horizon surfaces, arbitrary grids and along the well path. In the latter case the observation points may be at a pre-specified spacing along a well path or they may coincide exactly with the location of a logged fracture.

#### Retro-deformation

An ED model is usually described with all mechanical boundaries defined in terms of the present day, deformed, state. Faults will have undergone relative rotations, horizons have offsets at faults and well paths cross faults without interruption. A forward model from such a starting configuration will misplace horizons and introduce offsets in wells that cross faults; all observation points will be moved from their present-day configuration. To counter this, our models are retro-deformed, entirely elastically, using a numerical inversion process prior to being forward modelled. This process ensures that (a) faults are rotated to a pre-deformation configuration, (b) the forward model places horizons in their correct present-day positions and (c) wells have offsets at faults in the pre-deformed state but these offsets disappear during forward deformation. A useful byproduct of the retro-deformation is a restored model that might, in some circumstances, also be used for structural validation.

#### Modelling parameters, strain, stress and failure

All fault surfaces are represented by 20 m square elements each with a dip slip or slightly oblique slip displacement. The interface between the salt and the overburden is modelled with equilateral triangular elements (500 m sides), each with a unique dilation displacement.

For any given set of boundary conditions and observation points, strains in the model are influenced, additionally, by Poisson's ratio and distance of an observation point from the top surface of the half space, that is, the Earth's surface. Maintaining consistency with the geological observations, the modelling was performed with the top of the shallowest observed fault at 0 m, that is, the model coordinates have been uplifted by 1800 m from the their present-day locations. Poisson's ratio was set at 0.2, which is appropriate for normally compacted chalk.

Given the three-dimensional strains and Young's Modulus, the stress tensor at a point is calculated from the standard equations of linear elasticity. Following the Δχ method of Bourne & Willemse (2001), the combination of the stress tensor and a Mohr–Coulomb failure criterion distinguishes between three possible states of fracturing: (a) no fracturing; (b) tensile fracturing; and (c) shear fracturing. Here we use a cohesive strength of 2 MPa and a coefficient of internal friction of 0.55 (see Discussion).

#### Model configurations

In an attempt to address the sensitivity of model outcomes to uncertainty in both geological input and geological control, we have run nine basic configurations. Displacements on faults are the same for all models and vary on a per-fault basis, as shown in Figure 4b. The following primary configurations are defined:

faults:

include;

exclude;

uplift/subsidence from salt:

use the full extent inferred from the seismic interpretation;

use half of the observed uplift/subsidence;

subsidence in the exterior of the diapir:

subsidence owing to withdrawal of locally sourced salt;

no subsidence, reflecting the case where all salt has moved in from outside the area of interest.

The nine cases and their mnemonics are formalized in Table 1. Each of these cases has been run with a range of three different, constant volume, remote strains (see Table 2) and with the maximum principal extension in five different orientations (80, 90, 100, 110 and 120°). In what follows, a set of mechanical boundaries, remote strains and orientations combine to form what is referred to as a ‘model-run’.

#### Orientation comparison strategy between simulated and observed fractures

Comparison of simulated and observed fracture orientations is central to the validation and refinement of the modelling process. In a model run, either a single tensile fracture or a conjugate pair of shear fractures is simulated at the locus of every fracture interpreted from the image logs. For a simulated tensile fracture the difference in dihedral angle between its pole and the pole of the observed fracture is recorded as Δ*p*. The situation is slightly more complex for simulated shear fractures (overwhelmingly dominant in our model-runs) since there are two conjugate surfaces for each single observation. In this latter case, it is the minimum of the two possible dihedral angles that is recorded (also Δ*p*). The results of each of the individual comparisons are then represented in the form of histograms. This comparison strategy is expected to be robust because it cannot be influenced by the sample bias from the image logs discussed earlier.

We have chosen to use the comparison between poles rather than strikes because it is robust in measure in three dimensions. Difference in strike is not. For example, trivially, two mutually perpendicular surfaces with the same dip direction have the same strike. Consider further a pair of conjugate surfaces plunging in such a way that the difference in strike is 45°. The minimum difference in strike between a random surface and the conjugate pair has a 75% chance of being less than 45°. In general, comparison of strike angle between an observed fracture and a conjugate pair yields a distribution that is strongly skewed towards zero, giving an entirely false impression of goodness of fit. Similarly, the distribution of minimum differences in orientation of poles to a conjugate pair and a random surface is also not uniform. With this in mind we have to be careful when we interpret a Δ*p* distribution from one of the model runs.

First we state the null hypothesis that the distribution of Δ*p* is due to a random process. The test population for this case is readily approximated in the form of a cumulative density function (cdf). Given a coefficient of internal friction of 0.55, the distribution is bell shaped with a mean of 43.7° and slightly skewed towards 90° (Fig. 6a). We assess the goodness of fit between our modelled/empirical cdf and the test cdf by reference to the Kolmogorov–Smirnov test (see Davis 1986, p. 99). In particular we choose a level of significance of 0.05 meaning that the chance of incorrectly rejecting the null hypothesis is 5%.

#### Deformation intensity as an indicator of fracture density

In previous literature maximum Mohr–Coulomb shear stress has been proposed as a possible indicator for relative fracture density because of its direct interpretation in terms of the Mohr–Coulomb failure envelope (e.g. Maerten *et al.* 2002; Bourne & Willemse 2001). This quantity embodies the strains, their conversion to stress and the geometric properties of the failure envelope. However, the physical meaning of Mohr–Coulomb shear stress is quite remote from the immediate solutions of an ED model, which are displacements and strains. Accordingly, we have chosen to use displacement and strain as the quantities against which we compare measured fracture densities. Specifically, we employ an invariant measure of total strain,
coined as strain magnitude (see Brandon 1995), where the ε_{i} are the three principal strains at an observation point.

## Results

### Fracture orientation

Fracture orientation analysis combines the Δ*p* values from all domains in all wells. From configurations in Tables 1 and 2 and the five orientations of the principal horizontal extension, there is considerable variation in the distributions of Δ*p*. The largest variations are attributable to the bounding cases of orientations and magnitudes of remote strain. In general the configurations without influence of salt (S0F1) are consistently the poorest performers over the given range of remote strain values and orientations. The best fit from these models (least difference in poles) corresponds to an extension direction of 100° and a remote strain of ε_{v}=−0.06, ε_{hmin}=0, ε_{hmax}=+0.06, mean Δ*p*=35° with the median 30.7° and the 75th percentile at 46.2° (Fig. 6b). At the 5% level of significance we are comfortably able to reject the null hypothesis, concluding that the model predicts orientations that are better than random.

Of the model runs that do consider salt emplacement but not faults, all of these also lead us to reject the null hypothesis, leaving us with the alternative hypothesis that salt emplacement alone is also a better than random predictor of Δ*p*. The combination of faults and salt in the same model gives the best fits. All of these models enable us to reject the null hypothesis and they improve the mean Δ*p* and increase the skewness towards zero over the cases without faults or without salt. We note, however, that for the same remote strains the differences in mean and median Δ*p* between the best and worst fits is small, typically less than 5°. The configuration with the smallest Δ*p* is SBF1S1 with remote strains ε_{v}=−0.06, ε_{hmin}=0, ε_{hmax}=+0.06 at 100° having a mean Δ*p=*33.7° with the median 29.7° and the 75th percentile at 44° (Fig. 6c). This is our preferred model. Note that this model is the closest rendering of the actual interpreted geology.

### Fracture density

The relationship between observed fracture density, *D*_{f}, and ε_{m} varies between model runs but here we restrict the description to the preferred model from the orientation analysis, noting that from all runs this also gives us the most informative results for fracture prediction. For all locations of observed major fractures, ε_{m} has a small range between 0.06 and 0.25. The range of observed fracture densities is also small, 0.0–0.22 m^{−1}. Of the 25 structural domains, four have three or fewer data points and are not analysed further. In five domains there is a strong sense of a positive correlation between *D*_{f} and ε_{m} (Fig. 7a). Seven of the domains show small clumps of data with low variance in both *D*_{f} and ε_{m} (Fig. 7b). Four of the domains show large variance in *D*_{f} but little in ε_{m} (Fig. 7c), while the remaining five (Fig. 7d) have distributions that are either equivocal or hint at negative correlations. An important observation from Figure 7 is that, irrespective of the existence of a local correlation, the ranges of *D*_{f} and ε_{m} vary between domains.

Combining the results from all domains (Fig. 8) gives a clearer picture of a useful relationship. It confirms that the highest *D*_{f} are seen at the highest ε_{m} and vice versa. A linear regression has a slope of 0.24 and an intercept at *D*_{f}=−0.0045, endorsing the expectation that *D*_{f} ought to be zero at zero strain. This correlation is statistically significant (at a significance level of 5%) but weak with a coefficient of *r*=0.3. A more pragmatic interpretation is that the envelope of the results forms the upper limit of expected *D*_{f} for a given ε_{m}. We have also contoured the plot at 10, 50 and 90%. The contours follow the trend of the regression line at values of ε_{m}<0.015, tending to flatten at higher ε_{m}. So for an ε_{m} of 0.1 there is a 50% chance of observing densities greater than 0.02 m^{−1} while at an ε_{m} of 0.2 there is a 50% chance that densities will be greater than 0.06 m^{−1}.

An example of the spatial variation in ε_{m} is shown in Figure 9 at the Top M2 horizon for model run SBF1S1. The wells are subhorizontal and transect the horizon surface. The highest values appear in the crestal graben and in the footwall of fault F1. Visual inspection shows that, in general, in the regions coloured cyan/blue/violet, the density of major fractures is lower than that in the regions coloured green/yellow/red.

Subsequent to the modelling and analysis reported here a new well, N40, has been drilled into the footwall of F1 at the location shown in Figure 9. This well reported densities of natural fractures far in excess of what is typically expected in the low-fracture-density B Block. These enhance production to such an extent that the well now accounts for 25% of the daily production in the Gorm Field (Arnhild pers. comm.). While the well was not planned on the basis of this project, caveats notwithstanding, an area of predicted high strain without calibration by well data yielded a production sweet spot.

## Discussion

With trends in hydrocarbon exploration moving towards unconventional resources where exploitation depends on fracture permeability, there is a growing requirement to be able to forecast the characteristics of natural fractures at an early stage of the exploration and production process. Although seismic image processing can, in some cases, indicate small faults that would not ordinarily be picked on seismic sections, in general, fractures in the subsurface are not detected prior to drilling. At the resolution of the seismic method, such fractures are not only invisible but the seismic provides no way to determine if or how much the local rock volume has been deformed. The presence of deformation is inferred from non-planar horizons and offsets at seismically imaged faults. Deformation can be quantified directly at fault surfaces (e.g. Freeman *et al.* 2010) and it can be inferred, quantitatively, by kinematic or mechanical restoration techniques (e.g. Lohr *et al.* 2008) or by kinematic or mechanical forward modelling. In this work we have estimated the strains in the rock volume and their orientations using an ED, forward, mechanical model. We acknowledge that a homogeneous linear elastic rheology is unlikely to give a complete description of the real deformation that evolved over a period of 60 Ma. However, when the primary objective is to forecast the mechanical state of rocks prior to drilling, ED methods are strong contenders as pragmatic solutions.

In our best case model the simulated fracture orientations (poles) fit the observations considerably better than a model based on random variations with more than 50% of the Δ*p* less than 30°. Nevertheless, there are few Δ*p* values lower than 15°. This is of some concern. As noted previously, the nature of the Δ*p* distribution is influenced by the angle of internal friction used in the Mohr–Coulomb failure criteria. We ran sensitivity tests using extreme end-members of the coefficient of internal friction (μ=1, μ=0.2) that bracketed the value of μ=0.55 in the preferred model, and although the details of the Δ*p* distribution differed in each test, the numbers of Δ*p* less than 15° did not change significantly. It is possible that the lack of small Δ*p* might be due in part to spatially varying mechanical properties that we are unable to take account of. It is also possible that, since the Mohr–Coulomb criterion is linear and the orientation of the failure planes is always independent of the stress magnitudes, there is a variation in real conjugate angle that is not catered for by our model. Although beyond the scope of this paper, future work should explore and integrate (a) the effects of three-dimensional failure criteria (see Schöpfer *et al.* 2013), (b) the use of stress-dependent conjugate angles, as is the case in a Mohr failure envelope, and (c) the possibility that strain-controlled orthorhombic fault systems might form a more comprehensive template for the observed fracture orientations.

The cross-plotting of *D*_{f} against ε_{m} and its envelope trend provides a useful insight to the locations on a map of ε_{m} where the density of major fractures is relatively high/low, although there is considerable spread in the data. In comparable studies (e.g. Dee *et al.* 2007; Lohr *et al.* 2008) the field-scale calibrations appear simpler and tighter than we have been able to show. However, those studies refer to far fewer data points (five and four, respectively), where each point represents the average density measured in a near-vertical well over an entire reservoir interval. By contrast, in this work the fractures are nearly all identified from the subhorizontal parts of the well path over a total linear distance of more than 10 km across numerous stratigraphic/lithological divisions of the chalk reservoir. This has required a more general approach to the data collection and spatial analysis. Since between the mechanical boundaries the ED strain field is smooth and continuous, we argue that a strategy that samples *D*_{f} continuously is most appropriate inasmuch as it is best suited to capturing spatial variation. This leads to hundreds of points, each of which is approximately equivalent to a single point in the aforementioned papers. In some part, the spread that we observe is almost certainly influenced by the directional sampling bias, but we are unable to speculate on whether an unbiased sample would improve or degrade the correlation.

Our results show that consideration of the emplacement of salt is a significant factor in influencing the fit between observations and model. Mapping of the top of the salt dome was crudely defined at a linear scale of 500 m. It is possible that, with better seismic imaging, more accurate top salt mapping might improve the model's results but there are other factors that will affect the spread in Figure 9. First, although we have not systematically studied the issue of heterogeneity, we are aware of considerable spatial variation in the lithological and physical characteristics of the chalk reservoir. Clearly this should be considered in follow-up work, as should the application of geostatistics, especially collocated co-kriging of *D*_{f} with ε_{m}, with the aim of optimizing spatial correlations. Second, because the elastic deformation that we model is expected to be only a first-order approximation, there will of course be departures between the magnitude of the real strains compared with those calculated in our model.

Despite the above caveats, there appears to be a genuine correspondence between observed and modelled fracture orientations and densities, and although the correlations are not highly tuned, their general form indicates that prior to drilling it is possible to forecast orientations and relative densities of natural fractures.

## Conclusions

Strain calculated from an ED model that incorporates mechanical boundaries in the form of the tectono-stratigraphic interface with salt and tectonic faults can be correlated with fracture orientation and fracture density.

Given a broad range of tests of sensitivity to geological and tectonic controls, the geometry of our preferred model is simply the structure as interpreted and mapped from the seismic data.

Comparison of fracture strike is not a robust test of the similarity of orientation distributions of fracture surfaces. The difference in orientations between the poles to the planes of observed and modelled fractures is a precise measure.

More than 50% of the angular differences between the simulated fractures in our preferred model and observed fractures are less than 30°; 75% are less than 45°.

Strain magnitude is an effective predictor for fracture density. At the field scale low, non-zero, fracture densities correspond to small strain magnitudes and high fracture densities correspond to large strain magnitudes.

## Acknowledgments

We would like to acknowledge Maersk Oil, as part of the DUC partnership (Maersk Oil, Chevron and Shell), for permission to publish the data. The paper has benefited considerably from the comments of R. Shackleton and an anonymous reviewer.

- © 2015 The Author(s). Published by The Geological Society of London. All rights reserved