The formation of open fractures in cohesive materials, results of scaled analogue and numerical modelling on fault zone porosity development

We compare analogue and numerical models of dilatational fractures at low confining stress. These structures form an effective conduit for fluid flow in the field, but are difficult to model since they form in cohesive materials at low stresses. We use a truly cohesive powder for the analogue models and a Discrete Element Model (DEM) with brittle-elastic bonds for the numerical modelling. We show that despite variations in the model type, small differences in the location of initial fractures and the way these structures link-up to control the evolution of the model, the observed structures are robust. Three structural zones develop where different fault types dominate. In 3D numerical models we show an increase of the porosity on the fault zone with increasing deformation. The progradation direction is shown to be controlled by the position of the fracture. The combination of analogue models with cohesive powder and DEMs with internal cohesion is an excellent tool to study the evolution of open fractures.


Introduction
Brittle deformation of brittle cohesive materials under low confining stress can lead to the formation of dilatant fault sections.This can have profound effects on the fluid flow properties of the system.The diameter of these fault related cavities can easily be 5 orders of magnitude larger then that of pore spaces, and because the flux shows a third order dependence on the channel diameter these cavities can form extremely effective conduits for fluid flow if they are interconnected (e.g.[1][2][3][4][5][6][7]).
Carbonate reservoirs contain a large portion of the world's hydrocarbon reserves [1,[8][9][10], but good quality outcrops of massively dilatant fault zones in carbonates are rare.Outcrops in Tertiary carbonates on Jebel Hafeet, on the border between the United Arab Emirates and Oman, expose some examples in carbonates deformed at shallow depths.Jebel Hafeet is one of a series of foreland anticlines of the Oman Mountains [11].The normal fault zones in the area can be massively dilatant, with openings of several decimetres being common.Presently, these cavities are either filled with different material, such as carbonate veins, crushed wall rock or sediments, or they are still partially open (Fig. 1a and b). a e-mail : h.vangent@ged.rwth-aachen.de

EPJ Web of Conferences
On the mountain crest, wide surface fissures are found.These open structures, which strike parallel to the fold axis of the anticline (Fig. 1c and d), have apertures of more than a meter.Also observed are angular blocks of carbonate, dislodged and rotated between the parallel walls (Fig. 1c).However, despite first class outcrop conditions, the depth of the fissures on the mountain crest, and the connectivity of the cavities along the normal faults, remains difficult to determine.Results from Holland.
deformed cohesive basalts on the Big Island of Hawaii show similar structures and suggest that the dilatant segments of faults connect to a single (non dilatant) fault at depth, and that larger, unobserved cavities are located below rotated blocks [12].
The mechanical background of open fault formation can studied in analogue models using cohesive materials (see for example [5][6][13][14][15][16][17]) as well as with numerical methods, like the Discrete Elements Method [18][19][20].Here we will attempt to integrate these results as well as present some new insight on the progradation direction of these faults.
1 Model descriptions

Analogue model
The set-up chosen in all models is a graben bounded by normal faults with a dip of 60° [12,16].However, in the numerical models, this angle can be varied [18][19].The analogue model is 70 × 15 cm in plan view.The box was filled with 20 cm of modelling material.When a central table is pulled down by a motor, with a constant velocity, two side tables are pushed symmetrically to the sides, forming a master graben.The maximum horizontal elongation is 11% (see also [12,16] for more details).
The material used for the analogue modelling was dry hemihydrate powder (CaSO4 Ɣ ½ H2O), which is a truly cohesive, granular material with a small grain size (in the range of 10 to 400 ȝm).
The material properties of this powder are previously described [12,16].One particularly relevant property of this powder is that due to compaction (i.e. the reduction of porosity under increasing vertical stress) the powder becomes significantly stronger over the depth of the model (see [16]).This, combined with the increasing horizontal and vertical stress over depth, results in three structural zones in the model, a zone of tensile failure in the top 7 cm of the model, pure shear failure in the bottom 7 cm and a zone of hybrid-or mixed mode failure in between (see [16]).It is in this zone that dilational jogs form (sensu: [4,[21][22]).The transitions between these zones observed in the analogue models are fully consistent with predictions based on the material properties and stress calculations.
Deformation was documented with high resolution photographs.These images then served as the input of a Particle Image Velocimetry (PIV) analysis [12,16,23].PIV is a non-intrusive, optical method to determine displacements, and to calculate strains, and rotations.It is an excellent way to quantify analogue models, and to allow a better comparison with data obtained from numerical models [16].

Numerical model
For the numerical modelling the parallel Discrete Element (DEM) Software ESyS-Particle [24] was used.In the Discrete Element Method the brittle-elastic deformation of solid materials is modelled by a large number of spherical particles interacting with their nearest neighbours.The particle interaction can consist of breakable bonded interactions in case of cohesive materials or frictional interactions when there is no cohesion.The elastic and breaking parameters of the interactions can be scaled to match the material properties of the hemihydrate powder.However, the depth dependent increase of the cohesive strength has not been implemented so far.In the numerical models, only one half of the symmetrical set-up of the analogue models, i.e.only a single normal fault, was used in order to save calculation time.Although the models are essentially scale-less, the particles are significantly larger compared to the size of the model then the hemihydrate grains in the analogue model.Initially, 2D models were run, with particle number around 25000 to 70000, depending on the packing.
For a parameter study on the effects of changing cohesion, fault angle and packing, we refer to [18,20,25].

Comparison of the model results
Figure 2 shows a comparison between a numerical and an analogue model.This figure clearly shows that the first order structures are robust.Both models show basically three different trough going structures that are observed in nearly all models.The first (A) consists, in both models, of a near vertical section (i) at the top of the model, which changes dip at a certain depth and connects to the basement fault.The dipping section shows several cavities in the analogue model results (ii) and also in the numerical model.In nearly all models this is the first fracture to form.The middle fault (B) also has a vertical top section, but the deeper section connecting to the basement is a lot steeper than that of fault A. In the analogue model, this deeper section is less developed.The third fault (C) shows a peculiar, curved shape in the upper section (iii).Both faults have overall the same shape, but whereas in the analogue experiment this fault dips towards the centre of the box, the same fault in the numerical models dips towards the other end of the model.Between faults B and, a slightly dipping monocline forms, this, with continued deformation, will form a rotating block.Latter in the experiment, this block starts to disintegrate completely [16].
In both models, the upper section displays massive dilatancy, while the middle section has confined cavities.Toward the base of the model, the fault gets increasingly non-dilatant.Similar faults are inferred on Hawaii [12].Three different structural zones can be clearly identified in both figures, although there are differences in the exact position of the boundaries, as well as in the exact position of the three (in some experiments, two) fractures at the surface (A, B and C in the figure).The fractures A and B are the first to form, both in the analogue and numerical models, and the position of these fractures is variable with in a zone around the upper termination of the basement fault.Their position strongly controls the further evolution of the model, and which of these fracture ends up connecting to the basement.This then has strong influences on the final geometry.The variable position and number of the tensile surface fractures is observed in repeated analogue models and numerical models with different particle packings, but all other parameters constant.This suggests these variable positions are an inherent randomness in the model, which has a significant influence on the final geometry model.

3D Numerical models
Although the models described above show the development of the development of these faults in side-view in a lot of detail, it is the 3D connectivity of the cavities that controls the fluid flow through the fault zone.To get insight into this, the 2D numerical model was extended to a full 3D model with over 1.5 M particles.Because of the heavy computational load, only the first 4% strain was modelled.To visualize the evolution of open spaces in the model, the local porosity of the model was estimated by placing a grid over the model space and calculating the volume of the intersection between the DEM particles and each grid cell.A threshold filter was then used to identify empty spaces in the model.Results (see Fig. 3) show that there is a marked increase in porosity along the fault plane with continued slip.On the surface and on the fault plane it self, fractures are seen to open, and close again under the influence of gravity.
A protocol is currently being developed to better visualize these cavities and to the show the percolating networks.The main issue is that the current method is insensitive to the anisotropy in the 22016-p.5 shape of the empty spaces we are trying to detect, i.e. while the faults are extended in 2 dimensions they are relatively narrow in the dimension perpendicular to the fault plane.The method can therefore not always distinguish between open faults and the "natural" inter-particle porosity of the 22016-p.6 14th International Conference on Experimental Mechanics DEM which, due to resolution limitations, is of a similar magnitude than the aperture of small faults but of isometric shape.

Fracture propagation
Tensile fractures at the surface can form in two mechanism.When they form at the surface due to tensile forces, they propagate down into the hybrid failure zone and change dip to propagate down.An alternative model is that the fractures initiate at the base of the model in the non-diliatant failure zone, where the stresses are high.They then propagate upwards through the mixed mode zone, up to the tensile failure zone, where they change dip and become tensile before they breach the surface.Normal faults in Ethiopian basalts show a vertical, tensile segments at the top [26].It is argued these faults are surface nucleated and propagate down, based on field observations and Mohr-Coulomb calculations.On the other hand, field observations in south-western Iceland and numerical modelling support the opposite model of depth nucleated, upward propagating fractures [27].
Our models show both types of faults (Fig 4 ) Faults A and B are mostly surface nucleated, while fault C propagates up.This can be explained by considering the model as an ideal elastic beam that is flexed by an inclined fault in the middle and (while one side moves down) both ends are kept horizontal.An S-shaped fold will develop in such a beam, with the highest tensional stresses developing on the outer arches of this fold.The upper outer arch is located on the footwall block of the fault, which is the location of faults A (and B), while the lower outer arch is positioned in the footwall, which corresponds to the initiation point of fault C.

Discussion and Conclusions
The development of open fractures in brittle rocks has a significant effect on the flow of any kind of fluid.The modelling of these structures has to date been a challenge, but with the use of fine grained cohesive powders in analogue models and dedicated numerical models the intricacies of these structures are beginning to be understood.
The direct comparison of numerical and analogue models have shown that these models show surprisingly robust structures, despite the differences in the model set-up, the details of the material properties and the inherent random response to small differences in packing, porosity and density.In the numerical models, these small differences stems from size differences in the particles, while in the analogue models, small inhomogeneities of the model fill will have he same effect.
The depth of nucleation of a fracture seems to be a function of the position of the fracture relative to the master fault.This results from bending of the material prior to failure.It is shown [16] that the model material is able to deform elastically before failing.The combination of analogue modelling and PIV and the visualisation in numerical models allow for excellent visualisation of detailed evolution of fractures in these models.

Fig. 1 a
Fig. 1 a) Normal fault zone in a competent carbonate (Ca), with approximately 3 m offset, showing strongly variable internal structure, width of the fault cavities and clasitc infill.Material from a mechanically weaker, slightly more clayey carbonate layer (Cl) is included in the fault zone both between the up-and downthrown parts of the clastic deposits (a), as well as in cavities further down dip (b).Also note the empty cavity in the bottom (c).a) Sediment-filled open fracture on Jebel Hafeet.Shown are: the wall rock (A), a calite precipitation zone (B) and stratified, clay-rich sediments in the fracture (C).b and c) Open fractures on the crest of the mountain, showing massive dilatancy.Note the rotated block in c).All photos are taken on Jebel Hafeet by M.Holland.

Fig. 2 a
Fig. 2 a) Result of an analogue model of 20 cm pure hemihydrate after roughly 3% of elongation.b) Result of a numerical model after roughly 3% elongation.In both models the three structural zones, and several features are indicated.See text for details.

Fig. 3 a
Fig.3a) The 90% porosity iso-surface at 3% elongation of the numerical model, with an oblique view from above.On the surface, linear tenile fracture have formed.b) The same as a) but now with view from below.Note the pockets of porosity along the fault plane, and the linear shape of the tensile fractures.c) Same view as a), but now at 4% elongation.The continued deformation has caused the block between the tensile fractures to collapse and close the tensile fracture at its base (see arrow).d) Same view as b), but now at 4% elongation.The number and size of porosity pockets on the fault plane has increased, but no percolating network has formed.

Fig 4 a
Fig 4 a) PIV calculations of strain from an analogue experiment showing a surface nucleated fracture moving down.b) PIV calculations of strain from an analogue experiment showing a depth nucleated fracture moving up.Although these calculations are from different .experiments,they show the general rule that the surface nucleated fracture forms first, and away from the centre, in comparison to the depth nucleated fracture.The position of these two fractures correspond to fracture A and B of Fig. 2 respectively, and the plots are mirrored in respect to c).For a) and b), the zone of high strain is outlined in blue.c) Plot of a numerical experiment, where the broken bonds are colored for their time of fracturing.Note how the two vertical tensile fractures move down, and the curved fractures nucleate at depth.