Rock fractures organize as networks, exhibiting natural variation in their spatial arrangements. Therefore, identifying, quantifying, and comparing variations in spatial arrangements within network geometries are of interest when explicit fracture representations or discrete fracture network models are chosen to capture the influence of fractures on bulk rock behaviour. Treating fracture networks as spatial graphs, we introduce a novel approach to quantify spatial variation. The method combines graph similarity measures with hierarchical clustering and is applied to investigate the spatial variation within large-scale 2-D fracture networks digitized from the well-known Lilstock limestone pavements, Bristol Channel, UK. We consider three large, fractured regions, comprising nearly 300 000 fractures spread over 14 200 m2 from the Lilstock pavements. Using a moving-window sampling approach, we first subsample the large networks into subgraphs. Four graph similarity measures – fingerprint distance, D-measure, Network Laplacian spectral descriptor (NetLSD), and portrait divergence – that encapsulate topological relationships and geometry of fracture networks are then used to compute pair-wise subgraph distances serving as input for the statistical hierarchical clustering technique. In the form of hierarchical dendrograms and derived spatial variation maps, the results indicate spatial autocorrelation with localized spatial clusters that gradually vary over distances of tens of metres with visually discernable and quantifiable boundaries. Fractures within the identified clusters exhibit differences in fracture orientations and topology. The comparison of graph similarity-derived clusters with fracture persistence measures indicates an intra-network spatial variation that is not immediately obvious from the ubiquitous fracture intensity and density maps. The proposed method provides a quantitative way to identify spatial variations in fracture networks, guiding stochastic and geostatistical approaches to fracture network modelling.

The Lilstock outcrop in the southern Bristol Channel provides exceptional exposures of several limestone beds displaying stratabound fracture networks, providing the opportunity to create a very large, complete, and ground-truthed fracture model. Here we present the result of automated fracture extraction of high-resolution photogrammetric images (0.9 cm/pixel) of the full outcrop, obtained using an unmanned aerial vehicle, to obtain a spatially extensive, full-resolution map of the complete fracture network with nearly 350,000 ground-truthed fractures. We developed graph-based functions to resolve some common issues that arise in automatic fracture tracing such as incomplete traces, incorrect topology, artificial fragmentation, and linking of fracture segments to generate geologically significant trace interpretations. The fracture networks corresponding to different regions within the outcrop are compared using several network metrics and the results indicate both inter- and intra-network (layer to layer) structural variabilities. The dataset is a valuable benchmark in the study of large-scale natural fracture networks and its extension to stochastic network generation in geomodelling. The dataset also highlights the intrinsic spatial variation in natural fracture networks that can occur even in weakly-deformed rocks over relatively short length scales of tens of metres.

Natural fractures in subsurface rocks are a source of heterogeneity that impacts flow and transport behaviour. The presence of fracture discontinuities needs to be modelled explicitly due to observed deviations from the continuum assumption of porous media. The departures are due to both individual properties (such as aperture, infill, and roughness) and global network properties (such as topological summary and length distribution). Understanding flow patterns due to effects of rock fractures networks is essential formany applications such as exploiting hydrocarbons, geothermal heat extraction, subsurface nuclear waste storage, and water aquifer development. Assessing the impact of fractures in modelling studies requires fracture network data which is difficult to sample from seismic data (due to image resolution issues) and borehole data (owing to sparse sampling). Outcrop analogue data provide a means to sample networkswhile honouring both spatial position and topological relationships.

The 100 000 m 2 wave-cut pavement in the Bristol Channel near Lilstock, UK, is a world-class outcrop, perfectly exposing a very large fracture network in several thin limestone layers. We present an analysis based on manual interpretation of fracture generations in selected domains and compare it with automated fracture tracing. Our dataset of high-resolution aerial photographs of the complete outcrop was acquired by an unmanned aerial vehicle, using a survey altitude optimized to resolve all fractures. We map fractures and identify fracture generations based on abutting and overprinting criteria, and we present the fracture networks of five selected representative domains. Each domain is also mapped automatically using ridge detection based on the complex shearlet transform method. The automatic fracture detection technique provides results close to the manually traced fracture networks in shorter time but with a bias towards closely spaced Y over X nodes. The assignment of fractures into generations cannot yet be done automatically, because the fracture traces extracted by the automatic method are segmented at the nodes, unlike the manual interpretation in which fractures are traced as a path from fracture tip to fracture tip and consist of several connected segments. This segmentation makes an interpretation of relative age impossible, because the identification of correct abutting relationships requires the investigation of the complete fracture trace by following a clearly defined set of rules. Generations 1 and 2 are long fractures that traverse all domains. Generation 3 is only present in the southwestern domains. Generation 4 follows an ENE-WSW striking trend, is suborthogonal to generations 1 and 2, and abuts on them and generation 3, if present. Generations 5 is the youngest fracture set with a range of orientations, creating polygonal patterns by abutting at all other fracture generations. Our mapping results show that the northeastern domains only contain four fracture generations; thus, the five generations of the outcrop identified in the southwestern domains are either not all present in each of the five domains or vary locally in their geometry, preventing the interpreter from linking the fractures to their respective generation over several spatially separate mapping domains. Fracture intensities differ between domains where the lowest is in the NE with 7.3 m -1 and the highest is in the SW with 10 m -1, coinciding with different fracture orientations and distributions of abutting relationships. Each domain has slightly different fracture network characteristics, and greater connectivity occurs where the development of later shorter fractures is not affected by the stress shadowing of pre-existing longer fractures.

Multiscale karst systems, typically below seismic resolution, are important contributors to permeability but also provide significant challenges in drilling operations and in flow simulations. A substantial portion of subsurface karsts is of hypogenic origin, that is, is associated with upward movement of chemically aggressive fluids. The Morro Vermelho Karst System and Cave developed within Proterozioc shallow water carbonates and is an inspiring analog for buried hypogenic karsts. MVKS formed in association with a regional strike-slip fault. In an initial stage, deformation in the “softer” carbonates overlying more rigid quartzites and basement was accommodated by a 100s of m wide semi-ductile shear zone. Wholesale dolomitization took place and numerous dolomite-filled veins formed. With progressing strike-slip, the fault propagated upward creating an anticline parallel to the fault zone, causing the development of m-long extensional fractures, the arrival of SiO2 rich fluids in the karts and, eventually, the formation of the cave. Vertical fracture corridors and/or faults capturing deep, generally layer parallel flow are present in different settings, from foredeeps to rifted margins and are primary candidates for widespread hypogenic karst development.

The Morro Vermelho Cave (MVC) (Brazil) developed within the Morro Vermelho karst system, which affected Neoproterozoic limestones (Salitre Formation). The MVC experienced little interactions with meteoric processes and is an example of a hypogenic cave formed during strike-slip deformation. The Salitre carbonates in the MVC experienced distributed deformation along an elongated domain overlying a buried strike-slip fault. Gently dipping, semiductile shear zones formed with decimeter-scale (3.9 in.) dolomitic veins. In our model, Mg-rich fluids flowing along the Salitre aquifer caused at the same time extensive dolomitization of the body of rock (100-m [328-ft] scale) experiencing distributed deformation. With progressive displacement, the deep strike-slip fault propagated upward causing the development of an anticline pop-up, steepening sedimentary layers, and steep 1-10-m-long (3.3-33.8-ft) fractures, which served as pathways for upward fluid flow. These steep extensional fractures made it possible for fluids flowing in lower, quartzitic aquifers to enter the carbonate aquifer causing silica deposition in rock cavities and in fractures and fault planes. Following the main stage of speleogenesis, silica deposition took over again depositing on the cave walls a continuous silica crust, rarely observed in other settings worldwide. The interplay between regional bedding-parallel flow and focused circulation of fluids along steep faults and dipping layers, and the associated rock-fluid interactions are not unique to the contractional settings presented but can also occur in association with similar faults in rifted continental margins.

Natural fractures conduct fluids in subsurface reservoirs. Quick and realistic predictions of the fracture network organization and its fluid flow efficiency from limited amount of data is critical to optimize resources productivity. We recently developed a method based on multiple point statistics (MPS) technique to produce geologically-constrained fracture network simulations. The method allows to account for the intrinsic non-stationarity of these networks by considering a multivariate input data instead of averaged distribution of fracture parameters. In addition, the method considers probability maps reflecting the influence of fracture drivers in the network variability. Consequently, the simulated fracture networks derived from the innovative MPS approach are geologically better constrained than in classical discrete fracture network modelling approaches. This paper proposes to apply this method in subsurface conditions where available data are sparsely distributed. We developed a workflow where data are gathered from wellbore and from additional sources (outcrops). These data are used to extrapolate a network around the borehole as training images and themselves are extrapolated at the reservoir scale following a geological probability map.This work also presents innovations on the way how training images and probability maps that may integrate more geology constrain than relying almost entirely on available data.

Natural fracture network characteristics can be establishes from high-resolution outcrop images acquired from drone and photogrammetry. Such images might also be good analogues of subsurface naturally fractured reservoirs and can be used to make predictions of the fracture geometry and efficiency at depth. However, even when supplementing fractured reservoir models with outcrop data, gaps will remain in the model and fracture network extrapolation methods are required. In this paper we used fracture networks interpreted from two outcrops from the Apodi area, Brazil, to present a revised and innovative method of fracture network geometry prediction using the multiple-point statistics (MPS) method.

The MPS method presented in this article uses a series of small synthetic training images (TIs) representing the geological variability of fracture parameters observed locally in the field. The TIs contain the statistical characteristics of the network (i.e. orientation, spacing, length/height and topology) and allow for the representation of a complex arrangement of fracture networks. These images are flexible, as they can be simply sketched by the user.

We proposed to simultaneously use a set of training images in specific elementary zones of the Apodi outcrops in order to best replicate the non-stationarity of the reference network. A sensitivity analysis was conducted to emphasise the influence of the conditioning data, the simulation parameters and the training images used. Fracture density computations were performed on selected realisations and compared to the reference outcrop fracture interpretation to qualitatively evaluate the accuracy of our simulations. The method proposed here is adaptable in terms of training images and probability maps to ensure that the geological complexity in the simulation process is accounted for. It can be used on any type of rock containing natural fractures in any kind of tectonic context. This workflow can also be applied to the subsurface to predict the fracture arrangement and fluid flow efficiency in water, geothermal or hydrocarbon fractured reservoirs.

Representing fractures explicitly using a discrete fracture network (DFN) approach is often necessary to model the complex physics that govern thermo-hydro-mechanical-chemical processes (THMC) in porous media. DFNs find applications in modelling geothermal heat recovery, hydrocarbon exploitation, and groundwater flow. It is advantageous to construct DFNs from the photogrammetry of fractured outcrop analogues as the DFNs would capture realistic, fracture network properties. Recent advances in drone photogrammetry have greatly simplified the process of acquiring outcrop images, and there is a remarkable increase in the volume of image data that can be routinely generated. However, manually digitizing fracture traces is time-consuming and inevitably subject to interpreter bias. Additionally, variations in interpretation style can result in different fracture network geometries, which, may then influence modelling results depending on the use case of the fracture study. In this paper, an automated fracture trace detection technique is introduced. The method consists of ridge detection using the complex shearlet transform coupled with post-processing algorithms that threshold, skeletonize, and vectorize fracture traces. The technique is applied to the task of automatic trace extraction at varying scales of rock discontinuities, ranging from 10° to 102m. We present automatic trace extraction results from three different fractured outcrop settings. The results indicate that the automated approach enables the extraction of fracture patterns at a volume beyond what is manually feasible. Comparative analysis of automatically extracted results with manual interpretations demonstrates that the method can eliminate the subjectivity that is typically associated with manual interpretation. The proposed method augments the process of characterizing rock fractures from outcrops.

In carbonate rocks, channelized fluid flow through fracture conduits can result in the development of large and connected karst networks. These cavity systems have been found in multiple hydrocarbon and geothermal reservoirs, and are often associated with high-permeability zones, but also pose significant challenges in drilling and reservoir management. Here, we expand on the observed interplay between fractures, fluid flow and large cave systems, using outcrop analysis, drone imagery and fluid-flow modelling. The studied carbonate rocks are heavily fractured and are part of the Salitre Formation (750–650 Ma), located in central Bahia (NE Brazil). Firstly, the fracture and cave network data show a similar geometry, and both systems depict three main orientations, namely; NNE–SSW, NW–SE and ESE–WNW. Moreover, the two datasets are dominated by the longer NNE–SSW features. These observed similarities suggest that the fractures and caves are related. The presented numerical results further acknowledge this observed correlation. These results show that open fractures act as the main fluid-flow conduits, with the aperture model defining the fracture-controlled flow contribution. Furthermore, the performed modelling highlights that geometrical features such as length, orientation and connectivity play an important role in the preferred flow orientations.

Improving the understanding of fluid flow through fracture networks is crucial to the optimum recovery of hydrocarbons in naturally fractured reservoirs. In explicit dual continuum representations of matrix and fractures, it is not just the intensity of fracturing and the topological relationships in the network but also connectivity of hydraulically open fractures that determines effective permeability of the system. Fractured folds often form prolific reservoirs owing to the structural closure they afford and the additional porosity and permeability due to the fold related fracturing. The fracture patterns are of specific interest owing to the complex geometries associated with folding. In this work we utilize a combined outcrop based and numerical approach to characterize fracture patterns, fracture apertures and fluid flow sensitivities using a folded ‘box-type’ anticlinal structure example from the Pag Island, Croatia. The Pag Island is part of the main Adriatic-Dinaridic Carbonate Platform in the External Dinarides region that persisted from the Triassic to the Eocene with Paleocene-Miocene Dinaridic deformation phases. The current day geology consists of up to 1 km thick Cenomanian - Senonian shallow water rudists bearing carbonates and about 650m thick Eocene-Oligocene Nummulitic carbonate successions. We present a 3D structural model of the Pag Island with a Discrete Fracture Network (DFN) digitized representation built from drone photogrammetry. The multiscale fracture patterns in this multilayer folded reservoir analog are classified on the basis of folding stages (pre-, syn-, post folding). FEM geomechanical modeling is used to quantify stress-sensitive heterogeneous fracture apertures with the Barton-Bandis empirical model. Permeability sensitivities of the fracture network are tested using coupled single phase fluid flow and poroelasticity simulations. The results show that hydraulic apertures are extremely sensitive to fracture orientation with respect to the shortening direction and reservoir pressure transients. Our novel approach for coupled flow & geomechanics simulations yields effective permeability tensors that are functions of both matrix and fracture conductivities and also fracture aperture variations in time. We test dynamic closure & dilation of fractures in depletion and injection scenarios and our results stress the importance of dynamic fracture apertures in structurally complex fractured reservoirs.

Accurate prediction of net hydraulic fracture propagation pressure is often impossible. Under some conditions the pressure drop along the fracture dominates while in other cases the pressure at the tip determines the net pressure. Properly scaled lab tests and modeling indicate that effective stress determines the tip propagation pressure, but this is hard to confirm with field data in specific cases.

We therefore gathered a large data base of fracture treatments from many areas to investigate the correlation between net pressure and effective stress. In order to avoid any spurious effect from fluid friction, tortuosity and height containment we limited the data to relatively small injections with water or linear gel in vertical wells. All treatments were in conventional clastic reservoirs, but over a large range of permeability, rock stiffness and geological age.

The data show a remarkably good correlation between net pressure and effective reservoir stress, with a slope of 0.46. Lower net pressure of 200-300 psi was found in over pressured reservoirs and higher net pressure of about 1500 psi was seen in depleted reservoirs. We checked that this is not due to another underlying parameter, such as modulus or depth, which could explain the correlation. It is concluded that the correlation is due to a true relation between net pressure (controlled by fracture propagation) and effective stress. Simulation of representative treatments with a new model that includes a cohesive zone at the fracture tip shows excellent agreement with the observed correlation, supporting a physical relation.

The relation between net pressure and effective stress in the reservoir can contribute to improved treatment design in green fields and also will aid in understanding fracture height growth, since effective stress will differ between formation layers. Calibrated models will still be important in view of lack of detailed formation knowledge, but a correct description of the physics of fracture propagation, based on effective stress at the tip, will facilitate more accurate model predictions.

Pore pressure and its relationship with fracture net pressure has been reported qualitatively from both field and experimental observations. From a modeling perspective, the ubiquitously used pseudo 3D (P3D) models that are based on linear elastic fracture mechanics (LEFM) do not include the effect of reservoir depletion (or overpressure). Models that utilize effective stress as propagation criteria with a cohesive zone description, introduce the pore pressure directly into the simulation and hence can potentially capture the effect of pore pressure on fracture propagation. This work investigates the effect of pore pressure on hydraulic fracturing net pressure and geometry using empirical and numerical simulation approaches. We carried out an analysis of more than 400 datafrac injections spanning a wide range of geological ages and depositional environments in order to investigate the relationship between observed net pressure and reservoir pore pressure. The net fracture propagation pressure from the fracture treatment analysis was seen to be correlated with the effective stress in the reservoir. Fracture propagation simulations were performed using a coupled finite element – finite difference fracture simulator. The code uses a cohesive zone model (CZM) to describe fracture propagation. Four different effective stress scenarios were used to study the effect of effective stress on net pressure. The simulation results closely match the empirical relation between net pressure and effective stress as obtained from the analysis of actual frac treatment data. It is observed from the simulations that the magnitude of the effective stress also has an effect on the fracture geometry with a high effective stress leading to wider, shorter and more radial fractures. The derived empirical correlation is hence useful as a fracture design parameter. The datafrac net pressure diagnostics workflow in the pseudo 3D models can incorporate local tip pore pressure as a new pressure matching parameter. The pore pressure effect can thus explain high net pressures routinely observed in frac operations and also as a containment mechanism.