Coupled hydro-mechanical evolution of fracture permeability in sand injectite intrusions

Qun Gn, Derek Elsworth, Yixin Zho, Antonio Gripp, Andrew Hurst

a School of Geosciences, University of Aberdeen, Aberdeen, UK

b State Key Laboratory of Coal Resources and Safe Mining, China University of Mining and Technology, Beijing, China

c Department of Energy and Mineral Engineering, EMS Energy Institute and G3 Center, Pennsylvania State University, University Park, PA, USA

ABSTRACT Sandstone “injectite” intrusions are generally developed by the fluidization of weakly cemented sandstones and their subsequent injection into fractured reservoirs.In this work,a continuum coupled hydromechanical model TOUGH-FLAC3D is applied to simulate the discrete fracture network in large-scale sand injectite complexes. A sand production constitutive model is incorporated to consider the formation of sand through plastic deformation and its influence on evolution of fracture permeability. Overpressures in the fluidized sand slurry drives the injection with sand dikes intruded upwards, typically into previously low permeability “tight” mudstone formations. The contrast in poroelastic properties of the underlying weak sandstone and overlying injectite receptor directly affects the evolution of fracture aperture both during and after intrusion. Fluid drainage into the unconsolidated matrix may reduce the extent of fracture aperture growth, through the formation of shear bands. The results of this work have broad implications related to the emplacement of sandstone intrusions and subsequent hydrocarbon accumulation, maturation and then production.

Keywords:Sandstone intrusions Fracture permeability Sand fluidization Hydro-mechanical coupling

1. Introduction

Sandstone intrusions within sedimentary successions have been systematically investigated for decades, including the sand injection process (Beaudoin and Friès,1982; Parize,1988; Huang,1988;Parize et al.,2007a,b)and the mechanism of sand fluidization and dyke formation(Newsom,1903;Jenkins,1930;Bureau et al.,2013;Scott et al., 2013). Normally three different processes promote the development of active sand injection and subsequent intrusion:(1)Sand body over-pressurization by extensive forceful sand injection during the burial stage, until the pore pressure in the sandstone formation overcomes the strength of rock matrix (Jolly and Lonergan, 2002). In this process, the scale of induced fractures range from millimeter to kilometer in length,filled by a sand-fluid mixture(Vigorito and Hurst,2010);(2)Rapid infilling of unlithified sediments, with this leading to meter to kilometer scale deep fractures/faults suddenly opened or re-opened (Mode I) (Monnier et al., 2015). Here, the increasing volumetric strain leads to depressuring, which promotes settling and deposition of the fluid mixture into the fractures/faults (Scholz et al., 2009); and (3)Induced seismicity may facilitate the fluidization/liquefaction of shallow unconsolidated sand bodies where they lay beneath consolidated sediments (Obermeier, 2009). This may occur as the small horizontal stress difference at shallow depth promotes fluidization/liquefaction and the upward propagation into dykes/fissures, potentially forming sand volcanoes if the surface is reached (Montenat et al., 2004). Depending on the stress regime and lithology,the resulting sandstone intrusion will present itself in different forms and structures, including dykes, sills, cones, and saucer-shaped sills with distinct scale networks (Cartwright et al.,2008; Bureau et al., 2013). Usually, dykes in sand injectites are formed from intrusions into host rock,while residing in subsurface or exposing as outcrops with meter-scale width and length(Grippa et al.,2019).Consequently,the dykes could be feasibly represented by discontinuities as fractures in the model, defined by various sizes of aperture and length.From observations of injectite dykes at outcrop, local stress field during injection, presence and form of pre-existing faults, modes of tectonic compaction and host rock lithology appear to be the main factors controlling the architecture and morphology of the sand intrusion complexes (Monnier et al.,2014). Experimental investigations have attempted to recreate sand intrusions using compressed air in dry silica powder(Rodrigues et al., 2009). However, the low viscosity of the compressed air and high permeability of the silica powder render the system as insufficient to transport a significant sand concentration in the slurry,and the dykes terminate prematurely.Substitution by a cohesive mixture of sand and gelatin has been more successful(Bureau et al., 2013); however, results again show that sand intrusion is still limited in extent relative to that observed in situ.

Fig.1. Development of upward-migrating dyke intrusions due to fluid pressurization. The increase of fluid mixture pressure in the basal sandstone formation results in vertical fracturing and growth in Mode I failure. Shear bands develop in the fracture branch in the upper inclined dyke zone in a mixed mode of shear-tension failure.

In past characterizations of sand production, two- (2D) and three-dimensional(3D)finite difference models have been used to simulate the production of sand around the wellbore regime, by applying a strain-softening Mohr-Coulomb constitutive model(Detournay et al., 2006; Kim et al., 2011; Alquwizani, 2013). The previous investigations of injectite dyke intrusions decoupled the influence of pore pressure evolution with the stress regime,without considering the response of sand production in modifying the properties of the fluidized sand mixture.Furthermore,properly accommodating the evolution of the permeability in dykes is essential in understanding the path of fluid flow. The improved predictions of the architecture and permeability of such dynamically evolving injectite complexes will provide important insights in improving the efficiency of hydrocarbon recovery.Consequently,we address the following questions using a fully coupled hydromechanical model:

(1) What is the influence of sand fluidization in altering the evolution of fracture aperture?

(2) How the poroelastic contrast between the consolidated and unconsolidated formations affects the path of fluid flow and the resulting permeability change in the injectite complex?

The following investigates the evolution of aperture of dyke due to the transport of fluidized sand by incorporating the deformation of sand followed by fluidization,and accommodating the influence of the poroelasticity in modulating the fracture pressure and aperture growth.

2. Mechanical response to sand fluidization

In sand injectite complexes, unconsolidated sedimentary formations (sandstones) usually underlie mudrock seals with low permeabilities of the order of ～10-2mD to 10-5mD(Osipov et al.,2004). The overlying mudrock seal promotes the development of overpressurization in the underlying sand, resulting from compaction within the sequence. If sufficiently high, relative to in situ stresses, these overpressures will drive the fluidized sand mixture as a hydraulic fracture (Fig. 1). Fractures typically grow vertically by opening in Mode I, until a critical intrusion length is reached. The development of included dyke branches may also result in the development of shear bands. A sill zone is characterized by intrusions which are dominated by tabular sills. The increase in permeability resulting from the dykes provides conduits for fluid transport. The ensemble process is triggered by pore pressures in excess of the minimum principal stress and evolves only if the underlying sand is sufficiently friable that it may fluidize.

2.1. Fluidization triggers

The overpressurization in the sand body will facilitate failure of the seal, mostly by tensile/shear fracturing, through reducing the effective stress. A typical failure criterion involves the combined Griffith-Coulomb failure envelope, depending on the relationship between the differential principal stresses (σ1and σ3) and tensile strength of host rock,T(Jaeger and Cook,1979; Cosgrove,2001):

(1) σ1- σ3>5.66T:Shear failure occurs under high differential stress;

(2) σ1- σ3＜ 4T: Tensile failure under low differential stress;and

(3) 4T＜ σ1- σ3＜ 5.66T:Mixed mode tension and shear failure.

The sand fluidization process is governed by the fluid pressure,stress state, and fluid velocity adjacent to the pore surface in the sand body.

2.2. Fluidization mechanisms

Rocks in sandstone layer may disaggregate into sands in a form of shear or tensile failure.In non-cohesive sands at relative shallow depth where the total stresses are low, when the vertical pore pressure gradient within the sand exceeds the lithostatic gradient,the fluid fluidization velocity(Vmf)is reached(Bureau et al.,2013).The minimum fluidization velocity is defined in Eq.(1),modulated by grain diameterd, acceleration of gravityg, density of fluid ρf,density of solid grain ρs, and viscosity of the fluid mixture μ(Richardson,1971):

Typical minimum fluidization velocity for fine-grained and well-sorted sand formations in the size ranges (≤0.05 mm) are 0.3—1 cm/s in the water saturated medium (Eichhubl and Boles,2000). In this state, the sand will fluidize to dramatically increase hydraulic conductivity. Fig. 2 summarizes the schematic steps integrating the production of sand and subsequent fluidization process, transporting with fluids in pore networks. The sand grains may be transported through a connected pore system(Lowe, 1975). The elliptic particles surrounded by brown dots represent the processes of particle mechanical deformation and fluidization. The sand particle will be disaggregated after mechanical deformation, followed by the fluidization process. In this work,the sand production is first modulated by strain-weakening model with Mohr-Coulomb failure criterion. Usually the failure is triggered by fluid overpressure. Besides the artificial injection leading to overpressure in sand body, the disequilibrium compaction tectonic process could also result in excess fluid pressure. During the fast burial of low permeable sediments, the rapid burial rates prevent fluid escape from the local pore spaces,which leads to overpressure in a short time either. After the occurrence of failure, as soon as the fluid velocity reaches the minimum fluidization velocity, which means the drag force overcomes the effective weight of grains, the sand particles could become entrained with the fluid. Production of fluidized sand may drastically change the porosity of the fluidized medium,characterized by the mass balance relation (Vardoulakis et al.,1996):

where φ is the porosity;and Δmis the change of sand mass in the rock, which could be calculated as the difference between the fluidized sand massmfluidizeand the original,as deposited,sand massmdepo:

2.3. Fluidization flux

In this work, a continuum constitutive model is developed to describe sand fluidization, including formation failure and sand erosion process, and solved using the coupled flow-deformation code FLAC3D-Tough (Gan and Elsworth, 2016a, b) (Fig. 2). It is assumed that once being fluidized, the sand migrates when the local fluid pressure gradient overcomes a local stress-dependent threshold, which is determined by the grain size and stress state(Detournay,2009).In this work,a uniform distribution of grain size in the bottom formation is defined,as the primary location of sand fluidization target at the sandstone layer.In Darcy flow,this critical stress-dependent value can be calculated by the critical specific discharge velocityqcr(Detournay, 2009):

whereais a dimensionless coefficient determined experimentally,kis the permeability, φcris the critical porosity to trigger collapse and sand removal,Rpis the sand grain radius,cris the residual cohesion, φris the residual friction angle, and σeffnormis the effective normal stress. As the sand particle is removed from the boundary layer, the porosity of the rock increases.

When the local discharge magnitude normal to the exterior boundary surface exceeds the critical magnitude, sand production initiates. The rate of produced sand has a direct relationship with the evolution of rock porosity. Based on the assumption that the sand begins to be physically removed from the boundary layer as soon as the porosity reaches the critical porosity φcr,applying mass conversation enables evolution of the rock porosity φ using a modified version of Eq. (2) as

where λ is the erosion coefficient, used to characterize the spatial distribution of erosion potential in the rocks,which is modulated by the damage and sorting of the grains;quis the local discharge rate;SLis the local zone surface area exposed to flow;andVLis the local element volume.The erosion constant is defined by the governing equation, assuming a linear relationship between the erosion coefficient λ and the induced local plastic strain γp(Wan et al.,2003):

Fig. 2. Schematic steps for sand particles production due to mechanical failure, and followed fluidization process mixing with fluid.

where γcritis the threshold plastic strain to indicate the timing for the onset of erosion. When the plastic shear strain γpis less than γcrit, the erosion coefficient λ remains as zero; when the shear strain γpis greater than the upper bound strain value γmaxwhich indicates complete collapse of the rock, the erosion coefficient remains constant at γmax. The erosion coefficient is used to describe the degree of deformation occurred in rock matrix,which links the plastic deformation with the actual erosion process.Consequently,the weight of the produced sand massMsin local element is evaluated as a function of porosity change and erosion coefficient as(Detournay, 2009):

Once being fluidized,the density of the mixtureis updated based on the mass concentration of the fluidized solidscs:

2.4. Stress-dependent dyke aperture

Based on field observations,dike apertures usually range from a few centimeters up to 7 m.Despite that most dykes are oriented at steep angles, most of the sandstone volume is accommodated in the permeable (thick aperture) low angle intrusions, which represents that most of the sandstone volume is dominated with saucershaped intrusions (Grippa et al., 2019). In order to examine the evolution of dyke permeability due to the perturbation of pore pressure,the permeability of dyke is associated with the thickness of aperture as

whereVratiois the representative volume ratio of dykes imbedded in sediment elements;bis the aperture of dyke,which is subjected to evolve by different processes, including normal closure, shear dilation, and tensile opening; andbiniis the initial aperture of the fracture. In this work, the constitutive model for predicting dyke aperture is based on the simplified Barton-Bandis model(Baghbanan and Jing, 2007):

Fig. 3. Flowchart of sand production and stress-dependent dyke aperture evolution constitutive models implemented in the Tough-FLAC3D simulator.

2.5. Poroelasticity of dual medium

Furthermore, the pressure in matrix and fracture will be adjusted due to the deformation/compaction and the change of fluid modulus.The constitutive model for governing poroelasticity response in dual porosity system is relevant to the matrix volume and fluid modulus,defined as (Detournay and Cheng,1993):

whereBis the Skempton coefficient,Kis the bulk modulus,Kfis the fluid modulus,vis the matrix volume ratio,Krrepresents the rock modulus for fracture (1) and matrix (2), α is the coefficient in calculating the average pore pressure for the dual medium,andbis the Biot coefficient. In the scenario of porous mudstone, the high porosity engenders a lower magnitude Skempton coefficient,resulting in a lower magnitude of pore pressure build-up from Skempton mechanism.

3. Geologic setting and interpretation

Sand injectite complexes,including the Panoche Giant Injection Complex (PGIC), have been identified in different stratigraphic formations in the northern San Francisco area, Sacramento Basin,and San Joaquin Basin, USA (Huuse and Mickelson, 2004; Vigorito et al., 2008). The candidate reservoir for this work is selected from unit 3 (PGIC) in an intrusion of fluidized sand into hydraulically fractured slope mudstones of the Great Valley sequence,where is one of two adjacent exceptionally well-exposed largest outcrops(Grippa et al.,2019).Unit 3 comprises a 330—440 m thick mudstone-dominated interval of exclusively high-angle sandstone dykes. The sandstone intrusions in PGIC generally are formed by fluidization and injection of sand predominantly from turbiditic sandstone units in the Dosados Member(Vigorito and Hurst,2010).

Fig. 4 shows the 2D vertical geometry of high-angle intruded sandstone dykes in a surrounding tight mudstone formation in the PGIC.A permeability map is determined by measuring permeability on samples collected along and across the intrusions.The porosity and permeability for the parent sandstone formation are 0.3 mD and 480 mD, respectively, relative to the mudstone porosity of 0.1.Dykes penetrated the mudstone through intrusion with permeability and porosity of 30 mD and 0.25, respectively.

Fig. 4. Permeability distribution of intruded sandstone dykes pattern and sedimentary formation, including mudstone (grey) and sandstone (brown) (Grippa et al., 2019).

To interpret and quantify the pattern of dykes,such as orientation and spatial variation,an open source tool-box FracPaQ is applied to the 2D digital image (Healy et al., 2017; Rizzo et al., 2017). Fig. 5 shows the statistical distribution of dyke orientations and lengths,identifying that the major dykes are oriented in two orthogonal sets,70°and -20°against the horizontal direction, respectively.Manzocchi (2002) introduced the ternary plot of fracture connectivity with the 3 vertices of a triangle denoting I,Yand X nodes in the fracture network. Nodes are classified as ‘I’ (for isolated ends of traces),‘Y’(for branch points,splays or abutments)or ‘X’(for crosscutting intersections). FracPaQ plots two ‘contour’ lines of connectivity, forCL= 2 and 3.57, whereCLis the number of connections(intersections) per line (or trace). Fig. 5b shows the connectivity of trace segments in the dyke network with a major portion of I nodes,suggesting a relatively low connectivity of the network.

Fig.5. (a)Rose diagram for sand dykes’orientation data in the mudstone formation of the Great Valley sequence,California,USA,and the total analyzed fracture segments n=287;(b) Connectivity dyke segments in the target zone (the connectivity of dykes is dominated by I nodes); (c) Frequency of dyke segment orientations; and (d) Density of trace segments in target formation.

Fig. 6. Distribution of dykes in mudstone formation. Dykes are in two sets, 70° and -20° relative to the positive X-direction (horizontal), respectively. The lithology of bottom formation is sandstone,and the top layer is mudstone.The injection well is located at the sandstone formation with a constant injection pressure of 20 MPa.Reservoir dimensions are in meter.

Fig.6 shows a simplified 2D pattern of dykes with dip angles in the mudstone formation that follows the general undulating distribution of Fig. 4 and that honors the fracture statistics noted previously. This dyke pattern, comprising the pre-defined fracture permeability,porosity,and modulus,is input into the coupled flowdeformation simulator FLAC3D-Tough. The scale of 2D vertical section of reservoir geometry is 2300 m × 700 m × 10 m(length × width × depth). The lithology for the top layer is tight mudstone formation with intruded dyke network, while the underlying formation is porous sandstone.The injector located at the lower bottom formation provides enough constant pressure source to generate over-pressure,to push the flow moving upward into the fracture networks. Additionally, there is a constant stress gradient equal to 1.5×104Pa/m applied uniformly from the top layer to the bottom, while the hydrostatic pressure gradient is 104Pa/m. The observation point is defined in the coordinate of (920, 364) m,which intends to monitor the evolution of fracture permeability in top fracture networks,due to pressurization from downward layer.

4. Reservoir model set-up

Fig. 7 shows the initial distribution of stress state, mean permeability, and pore pressure across the entire reservoir, based on the prescribed gradients of stress and pressure, and geometry listed in Fig.6.Sandstone dykes are the main conduits for fluid flow in the upper mudstone formation to the surface. In this study, we explore the roles of sand fluidization and poroelastic contrast induced from the differences in matrix and fluid properties in evolution of dyke permeability. Consequently, the simulation scenarios are designed as shown in Table 1, highlighting the effect of both sand fluidization and matrix mechanical properties. The comparison between Case 2 and the Base Case 1 illustrates the influence of poro-elasticity induced from the lithologic variation,while Case 3 highlights the influence of sand fluidization and transport on the evolution of permeability. Table 2 shows the applied rock and fluid properties for the Base Case 1 used in this work. A high friction angle of 40°and cohesion of 4 MPa in mudstone formation are defined, in order to constrain the major plastic deformation in the bottom sandstone formation, since the major concerns in this work focus on the formation and transport of sand flux from the bottom layer.

Fig.7. (a)Applied horizontal boundary stress,Sxx;(b)Applied vertical boundary stress,Szz;(c)Initial pore pressure distribution(5.28—12 MPa);and(d)Distribution of initial mean permeability.

Table 1Modeling design schedule.

Table 2Mechanical properties of matrix and fractures used in the simulations (Rutqvist et al., 2013; Grippa et al., 2019).

Fig. 8. Evolution of fracture permeability at the observation point coordinate (920,364) m for Cases 1 (Base Case, black), 2 (red), and 3 (blue).

5. Results and discussion

Fig. 9. Evolution of matrix porosity comparison at the observation point coordinate(920, 364) m for Cases 1 (Base Case, black) and 2 (red).

Based on the prescribed three simulation cases, the analysis primarily focuses on the magnitude of permeability evolution and the mass of the produced sand, to illustrate the influence of sand fluidization. The properties of the dual-medium of the mudstone where the fracture network is embedded exert a key influence on determining the response of the fracture aperture.To highlight the effect of the poroelastic contrast due to lithologic variation, the Case 2 scenario is employed to replace mudstone with the same unconsolidated lithology as the basal sandstone, using an initial porosity of 0.3 and low modulus.

Fig. 8 shows the evolution of fracture permeability at the coordinate (920, 364) m within the upper mudstone formation. Based on the gradient of permeability enhancement in temporal scale,the plastic shear failure occurrence after 106s generated significant potential in enhancing fracture permeability from 4.2×10-14m2to 1×10-14m2.The comparison indicates that the poroelastic effect is significant in controlling the growth of the fracture aperture.Even though the initial permeability of Case 3 is identical to the two other scenarios,the higher porosity matrix in Case 3 yields a larger potential in exchanging fluid between fracture and matrix. The greater dissipation of pore pressures into the matrix reduces the potential of pressure build-up in the fracture. Consequently, the aperture growth gradient in Case 3 (red curve) is slower than that in the other two cases.

Fig.9 compares the evolution of matrix porosity between Cases 1 and 2 at coordinate(920,364)m.Since the porosity increment is caused by plastic failure during sand production,there is a reduced potential in the porous matrix in inducing plastic failure. The porosity enhancement in Case 2 is flatter than that in Case 1.As the pore pressure increases sufficiently to develop failure,the porosity increases significantly after ～106s (～12 d).

Fig. 10 presents a comparison of the fracture network permeability distribution after 1 month of injection. It indicates that the porous matrix in Case 2 counters the trend in increasing fracture network permeability. The majority of the fracture network permeability has not been elevated to the same degree as the other two cases. The development of shear bands in the matrix of the upper dyke zone actually inhibits further increase in permeability.A comparison of the changes in permeability between Cases 3 and 1 indicates that the sand production and fluidization process mediate the permeability enhancement. This response may result from the process of sand fluidization,which triggers an increase in porosity.Furthermore, fluidized sands, mixed with water, will increase the average density of flow and the fluid moduli. Consequently, the applied stress will lead to a larger pressure build-up. In this work,we did not consider the effect of fluidization in modifying the viscosity of the flowing mixture, although it is anticipated that effective viscosity will increase with the concentration of sand particles (Gibilaro et al., 2007).

Fig.10. Spatial distribution of fracture network permeability (m2) in the upper mudstone formation after 1 month of injection for Cases 1, 2, and 3.

Fig.11. Comparison of matrix porosity in the fracture network zone after 1 month of injection for Cases 1 and 2.

To highlight the influence of sand production and fluidization in modifying matrix porosity, the evolutions of matrix porosity and fracture porosity are listed separately. The sand will be produced and fluidized once the porosity φ exceeds the critical porosity value φcr=0.5(Eq.(7)).Fig.11 shows a comparison of the matrix porosity for Cases 1 and 2 in the fracture network zone. This indicates that the matrix porosity in Case 1 evolves significantly from the initial state, while the matrix porosity in Case 2 only increases slightly from 0.3 to 0.33. Most of the porosity-enhanced areas are within the zones of the fractured block, since the fluid mixture with produced sands is transported along the major fractures.Fig.12 further indicates the mass distribution of the fluidized sand after 1 month of injection activity. Since the injector is located in the basal sandstone formation, the fluidized sand particles primarily originate from this zone, driven by pressurization. The relative porous and low-modulus material in Case 2 increases the potential for fluidization,in terms of the magnitude of the induced plastic strain to the critical threshold magnitude.

Fig.12. Distribution of the mass of fluidized sand (g) in the reservoir after 1 month of injection for Cases 1 and 2. Reservoir dimensions are in meter.

6. Conclusions

This work explores the mechanisms of sand emplacement by fluidization and injection contributing to the formation of sand injectite complexes. This is accomplished by inserting an appropriate constitutive model into the framework of the coupled flowdeformation continuum model Tough-FLAC3D. This study focuses on the evolution of permeability for intruded dyke networks,particularly accommodating the transport of sand particles. The fluidization process is considered as a subsequent step following the disaggregation of sand after its failure.When the fluid velocity reaches a critical fluidization velocity, the sand will deform uncontrollably as the grains become entrained within the fluid. The magnitude of erosion process actually determines the timing and mass of the fluidized sands. The higher the fluid velocity and the more weakly cemented the matrix, the stronger the fluidization.

Reservoir-scale sand injectite simulations reveal that poroelastic effects exert a strong control on the growth of fracture aperture.Higher permeabilities of the matrix tend to dissipate overpressures from the fractures and reduce the propagation of the injectite dykes.The increment of porosity increase is not as significant for a porous matrix as for the tight matrix condition. Moreover, the development of shear bands in the dyke zone also inhibits the growth in fracture permeability. In this work, sand fluidization initiates once the porosity increases to a critical magnitude. Since the erosion coefficient is linearly related to the magnitude of the induced plastic strain, the propensity for sand fluidization is principally within the unconsolidated sandstone formation after pressurization. If the matrix is unconsolidated with low deformation modulus, the potential for fluidization becomes more significant.Moreover, the fluidization and subsequent sand transport process prompt the development of fracture permeability by increasing either fluid density or modulus. A higher fluid modulus tends to build up greater pore pressure,to yield higher growth of aperture.This study provides an understanding of the controls on the evolution of permeability and porosity in sand intrusion reservoirs —contributing to our understanding of permeability structure that may aid in the recovery of hydrocarbons.

Declaration of Competing Interest

The authors wish to confirm that there are no known conflicts of interests associated with this publication and there has been no significant financial support for this work that could have influenced its outcome.

Acknowledgments

The authors would like to thank the support in using FracPaQ from Roberto Rizzo in the University of Aberdeen. We also appreciate the financial support from the Laboratory of Coal Resources and Safe Mining (China University of Mining and Technology, Beijing) (Grant No.SKLCRSM16KFC01).