## Publications

253 results found

Bird R, Paluszny A, Thomas RN,
et al., 2023, Modelling of fracture intensity increase due to interacting blast waves in three-dimensional granitic rocks, *International Journal of Rock Mechanics and Mining Sciences*, Vol: 162, Pages: 1-15, ISSN: 0020-7624

The complexity of the physics of rock blasting is a longstanding modelling challenge. This work presents in detail a three-dimensional, material non-linear finite element based model for wave propagation, combined with a postprocessing procedure to determine the fracture intensity caused by blasting. The rock is described with the Johnson–Holmquist-2 constitutive model, an elastoplastic-damage model designed for brittle materials undergoing high strain rates and high pressures and fracturing; it is also combined with an instantaneous tensile failure model. Additionally, material heterogeneity is introduced into the model through variation of the material properties at the element level, ensuring jumps in strain. A detailed algorithm for the combined Johnson–Holmquist-2 and tensile failure model is presented and is demonstrated to be energy-conserving, and is complemented with an open-source MATLABTM implementation of the model. A range of sub-scale numerical experiments are performed to validate the modelling and postprocessing procedures, and a range of materials, explosive waves and geometries are considered to demonstrate the model’s predictive capability quantitatively and qualitatively for fracture intensity. Fracture intensities on 2D planes and 3D volumes are presented. The mesh dependence of the method is explored, demonstrating that mesh density changes maintain similar results and improve with increasing mesh quality. Damage patterns in simulations are self-organising, and form thin, planar, fracture-like structures that closely match the observed fractures in the experiments. The presented model is an advancement in realism for continuum modelling of blasts as it enables fully three-dimensional wave interaction, handles damage due to both compression and tension, and relies only on measurable material properties.

Walding JC, Paluszny A, Zimmerman RW, 2023, Numerical modelling of the influence of tidal stresses on qualitative fracture patterns on the surface of Europa

The ice crust of Jupiter's Galilean satellite Europa exhibits a number of large-scale lineae features, as well as numerous regions with smaller scale lineae patterns. In this work, a three-dimensional finite-element simulator is used to model lineae as fractures in the crust that nucleate, grow, and interact. The medium is assumed to be isotropic and linearly elastic. Fractures are assumed to grow primarily in tension due to tidal stresses in the ice crust, and a damage criterion is used to model the weakening of the ice matrix that occurs concurrently with fracturing. The growth of multiple fractures is modelled geometrically as a function of multi-modal stress intensity factors computed at the fracture tips. The tidal forces that drive this fracturing process are computed according to the model of Wahr et al. (2009), and fracturing is evaluated over multi-scale periods from days to millions of years. Fracture nucleation and growth are modelled within the span of the satellite, with emphasis on characterizing the fracture patterns in the equatorial region. Multiple three-dimensional non-planar fractures are seen to grow and interact within each region. The simulated patterns are qualitatively compared against images obtained by NASA's Galileo mission.

Chaiwan P, Burtonshaw JEJ, Thomas RN, et al., 2023, A three-dimensional, finite element-based study of the effect of heterogeneities on thermo-hydro-mechanical deformation during cold fluid injection

This paper investigates the effects of the injection of (relatively) cold CO2 into a saline aquifer, focusing on the effects that local heterogeneities may have on the geomechanical changes observed during injection. Representative properties of the Smeaheia site, a full-scale carbon storage site in the North Sea, are used in the simulation. The geometry of the Smeaheia site is represented by a simplified three-dimensional geometric model containing four rock layers, and a fault that spans multiple layers. Thermo-poro-elastic deformation is modeled using a finite element approach, with the coupled solution of fluid flow, heat transfer, and mechanical deformation of the Smeaheia site. The boundary conditions are based on in situ conditions at the Smeaheia site, along with fixed injection temperature and flow rate. The effects of injection rates, injection periods, injection temperatures, heterogeneities, fault, and fracture properties are investigated. The results show that the changes in temperature, pressure, and stresses are concentrated in the vicinity of the injection point, with a limited effect on the caprock away from the injection well. Thermal stresses dominate the overall stress reduction, shifting the stress state towards the failure condition. The flow rate causes a localized effect on the pressure field, due to the high permeability of the reservoir. Heterogeneity in the Young's modulus yields variations in deformation which are of comparable magnitude to the effects of cold injection.

Navidtehrani Y, Betegón C, Zimmerman RW,
et al., 2022, Griffith-based analysis of crack initiation location in a Brazilian test, *International Journal of Rock Mechanics and Mining Sciences*, Vol: 159, Pages: 1-16, ISSN: 0020-7624

The Brazilian test has been extremely popular while prompting significant debate. The main source of controversy is rooted in its indirect nature; the material tensile strength is inferred upon assuming that cracking initiates at the centre of the sample. Here, we use the Griffith criterion and finite element analysis to map the conditions (jaws geometry and material properties) that result in the nucleation of a centre crack. Unlike previous studies, we do not restrict ourselves to evaluating the stress state at the disk centre; the failure envelope of the generalised Griffith criterion is used to establish the crack nucleation location. We find that the range of conditions where the Brazilian test is valid is much narrower than previously assumed, with current practices and standards being inappropriate for a wide range of rock-like materials. The results obtained are used to develop a protocol that experimentalists can follow to obtain a valid estimate of the material tensile strength. This is showcased with specific case studies and examples of valid and invalid tests from the literature. Furthermore, the uptake of this protocol is facilitated by providing a MATLAB App that determines the validity of the experiment for arbitrary test conditions.

Setiawan NB, Zimmerman RW, 2022, Semi-analytical method for modeling wellbore breakout development, *Rock Mechanics and Rock Engineering*, Vol: 55, Pages: 2987-3000, ISSN: 0723-2632

Borehole breakout initiation, progression, and stabilization are modeled using a semi-analytical method based on Melentiev’s graphical conformal mapping procedure, and Kolosov–Muskhelishvili complex stress potentials. The only input data required are the elastic moduli of the rock, the Coulomb strength parameters (cohesion and angle of internal friction), and the far-field stresses. The stresses around the borehole wall are computed, the region in which the rock has failed is then “removed”, creating a new borehole shape. This process is iterated until a shape is obtained for which the breakout will progress no further, and a stable state has been reached. This modeling shows that stresses around the flank of the breakout evolve so as to reduce the propensity for shear failure, which helps to explain why the breakout width remains relatively constant throughout the process, even as the breakout region deepens radially. The failed area around the borehole becomes smaller and more localized, as the breakout tip sharpens and deepens. Using the Mogi–Coulomb failure criterion, a good match is obtained between the modeled breakout geometry and the geometry observed by Herrick and Haimson in laboratory experiments on an Alabama limestone. The new method leads to a correlation between breakout geometry, rock strength properties, and in situ stress. The paper ends with a critical discussion of the possibility of inferring the in situ stress state from observed breakout geometries.

Saceanu MC, Paluszny A, Zimmerman RW,
et al., 2022, Fracture growth leading to mechanical spalling around deposition boreholes of an underground nuclear waste repository, *International Journal of Rock Mechanics and Mining Sciences*, Vol: 152, ISSN: 0020-7624

This study presents a three-dimensional numerical analysis of multiple fracture growth leading to spalling around nuclear waste deposition boreholes. Mechanical spalling due to stress amplification after drilling is simulated using a finite element-based fracture growth simulator. Fractures initiate in tension based on a damage criterion and grow by evaluating stress intensity factors at each fracture tip. Tip propagation is multi-modal, resulting in final fracture patterns that are representative of both tensile and shear failure. Their geometries are represented by smooth parametric surfaces, which evolve during growth using lofting. The corresponding surface and volumetric meshes are updated at every growth step to accommodate the evolving fracture geometries. The numerical model is validated by comparing simulated fracture patterns against those observed in the AECL Underground Rock Laboratory Mine-By Experiment. It is subsequently calibrated to simulate fracture initiation and growth around boreholes drilled in the Forsmark granodiorite, subjected to a far-field anisotropic triaxial stress that corresponds to the in situ stress model from the Swedish Forsmark site. The deposition tunnel is implicitly simulated by attaching the deposition borehole to a free domain boundary.Several geomechanical cases are investigated, in which fracture growth is numerically evaluated as a function of in situ stress state, tunnel orientation, borehole geometry, total number of boreholes and borehole spacing. Numerical results show that spalling occurs in all cases, given the underground conditions at Forsmark, with borehole geometry, spacing and stresses affecting the extent of fracture nucleation and growth patterns.The uncertainty in underground stress conditions is evaluated through varying stress magnitudes and orientations relative to the tunnel floor. Whereas tunnel orientation influences the relative locations where fractures initiate with respect to the tunnel floor, frac

Burtonshaw JEJ, Paluszny A, Thomas RN, et al., 2022, The influence of hydraulic fluid properties on induced seismicity during underground hydrogen storage

Geological hydrogen storage is a promising technique to seasonally store surplus renewable energy at the large scale. Depleted gas reservoirs are currently being considered as potential storage sites for hydrogen. Induced seismicity is a key aspect of the sustainability and acceptance of both onshore and offshore gas storage and production. This study conducts a preliminary numerical assessesment of the deformation of a single three-dimensional fault that extends from the reservoir and into the caprock during the injection and storage of hydrogen over a period of weeks. Simulations are conducted with a three-dimensional finite element, fully coupled poroelastic code. Induced seismicity is quantified both in terms of tangential displacements on the fault surfaces, cumulative moment release, and maximum seismic event magnitudes. The reservoir is an anticlinal structure bound from above by a low permeability caprock, and below from a low permeability underburden. Mechanical properties are defined based on North Sea depleted gas reservoirs. The effect of injection fluid parameters such as density and viscosity are investigated. For hydrogen storage in depleted oil and gas reservoirs with proximate faults, some minor seismicity can be expected. According to these preliminary simulations, geological hydrogen storage appears to be operationally sustainable in terms of the magnitude of induced seismicity, at least for the injection volumes, rates and storage periods considered in this work. For the studied case, the maximal potential seismicity observed is of moment magnitude 2.1, with very few seismic events.

Thomas RN, Bird RE, Paluszny A, et al., 2022, Numerical Study of Three-dimensional Blast-induced Damage Patterns resulting from Simultaneous Borehole Blasting of Hard Rocks

Simultaneous detonation of charges in closely spaced boreholes is a commonly used blasting technique for fragmentation and construction. Blasting is a challenging phenomenon to model, due to the complexity of the mechanical deformation, fracturing, and fragmentation, and the spatial scales involved. Blast wave models must consider the possibility of constructive interference between separate waves in three-dimensions. A three-dimensional finite element method is used to study the induced damage in a general hard rock tunnel blast setting. The Johnson Holmquist-2 elastoplastic-damage model is used to quantify shear and tensile failure. In simulations with two blastholes separated by up to one meter, damage patterns emanating from boreholes interact to form self-organizing 'fracture-like' structures. Mechanical interaction between the two blasts is a function of the input charge wave properties, blasthole separation, and distance along the charge, with high concentrations of damage at the free boundary representing the tunnel wall. Constructive interference between the two blast waves is not shown to directly induce additional damage zones, and instead, interaction results from the overlap and slight extension of each blasthole's damage zone towards the other.

Zimmerman R, Bhullar A, Stewart G, 2021, One-dimensional flow to a constant-pressure boundary in a pressure-sensitive porous medium, AGU 2021 Fall Meeting, Publisher: Springer

In some “pressure-sensitive” rock formations, the variation of permeability with pore pressure is sufficiently large that it needs to be accounted for. Accounting for this variation of permeability renders the pressure diffusion equation nonlinear, and not amenable to exact analytical solutions. A constitutive model that fits many data sets is that of a permeability that varies exponentially with pressure. One-dimensional flow to a constant-pressure boundary in a medium whose permeability varies exponentially with pressure is investigated in this work, using numerical and approximate approaches. As part of the analysis, a dimensionless “pressure-sensitivity” parameter β is defined as the logarithm of the ratio of the far-field permeability to the permeability at the outflow boundary.The problem can be approached by applying a Boltzmann transformation to the governing PDE, and converting the resultant ODE into an integral equation, which can be solved by an iterative procedure. This numerical solution can serve as a benchmark against which to judge the accuracy of approximate solutions. One simple approximate model is based on the “sorptivity approximation” developed by Philip for unsaturated flow, which can be adapted to any nonlinear diffusion problem. This model linearizes the problem by defining an effective permeability as a drawdown-weighted average of the permeability as a function of pore pressure. The Philip approximation predicts that the flowrate, normalized against the case in which the permeability is constant, is 1–0.333β, to first-order in β.This problem can also be approached as perturbation problem, with β as the perturbation parameter. The zeroth-order solution is the classical error-function solution for a 1D diffusion problem with a constant-pressure boundary. The first-order correction is obtained in closed form, from which the flowrate at the outflow boundary is found. To first-order

Lutz MP, Zimmerman RW, 2021, The effect of pore shape on the Poisson ratio of porous materials, *Mathematics and Mechanics of Solids*, Vol: 26, Pages: 1191-1203, ISSN: 1081-2865

A brief review is given of the effect of porosity on the Poisson ratio of a porous material. In contrast to elastic moduli such as K, G, or E, which always decrease with the addition of pores into a matrix, the Poisson ratio ν may increase, decrease, or remain the same, depending on the shape of the pores, and on the Poisson ratio of the matrix phase, νo. In general, for a given pore shape, there is a unique critical Poisson ratio, νc, such that the addition of pores into the matrix will cause the Poisson ratio to increase if νo<νc, decrease if νo>νc, and remain unchanged if νo=νc. The critical Poisson ratio for spherical pores is 0.2, for prolate spheroidal pores is close to 0.2, and tends toward zero for thin cracks. For two-dimensional materials, νc=1/3 for circular pores, 0.306 for squares, 0.227 for equilateral triangles, and again approaches 0 for thin cracks. The presence of a “trapped” fluid in the pore space tends to cause νc to increase, and for the range of parameters that may occur in rocks or concrete, this increase is more pronounced for thin crack-like pores than for equi-dimensional pores. Measurements of the Poisson ratio therefore may allow insight into pore geometry and pore fluid. If the matrix phase is strongly auxetic, small amounts of porosity will generally not cause the Poisson ratio to become positive.

Burtonshaw J, Thomas R, Paluszny A, et al., 2021, Effect of Viscosity and Injection Rate on the Tensile and Shear Deformation Experienced By Interacting Hydraulic and Natural Fractures, 55th U.S. Rock Mechanics/Geomechanics Symposium - American Rock Mechanics Association (ARMA)

Alsuwaidi ES, Xi G, Zimmerman RW, 2021, Mechanical characterization of Laffan and Nahr Umr anisotropic shales, *Journal of Petroleum Science and Engineering*, Vol: 200, Pages: 1-16, ISSN: 0920-4105

Geomechanics-related wellbore instability has become a major source of non-productive time in highly deviated wells drilled in oil and gas fields in offshore Abu Dhabi. These wells are drilled through two highly anisotropic shale formations, namely the Laffan shale formation and the Nahr Umr shale formation. Most of the models used in the oil and gas industry do not account for the strength and elastic anisotropy of shale. Therefore, a laboratory study was conducted to examine the strength and elastic anisotropy of shales using uniaxial compression tests and triaxial compression tests.Jaeger's Plane of Weakness (JPW) model was used to understand the anisotropic failure behavior of highly laminated shales of the Laffan and Nahr Umr formations. This model assumes that for an anisotropic rock, there exists a plane of weakness that has strength properties (cohesion and angle of internal friction) different than the strength properties of the intact rock. By minimizing the Root Mean Square Error (RMSE), the experimental strength values of the samples, as measured at different orientations, were fitted to the JPW model.Elastic moduli were also measured on these shales, as a function of orientation angle. The results showed that the moduli vary with angle according to the expected tensor transformation law. Therefore, the transverse isotropy assumption is a reasonable model to be used when dealing with these laminated sedimentary rocks.

Bhullar AS, Stewart GE, Zimmerman RW, 2021, Perturbation solution for one-dimensional flow to a constant-pressure boundary in a stress-sensitive reservoir, *Transport in Porous Media*, Vol: 137, Pages: 471-487, ISSN: 0169-3913

Most analyses of fluid flow in porous media are conducted under the assumption that the permeability is constant. In some “stress-sensitive” rock formations, however, the variation of permeability with pore fluid pressure is sufficiently large that it needs to be accounted for in the analysis. Accounting for the variation of permeability with pore pressure renders the pressure diffusion equation nonlinear and not amenable to exact analytical solutions. In this paper, the regular perturbation approach is used to develop an approximate solution to the problem of flow to a linear constant-pressure boundary, in a formation whose permeability varies exponentially with pore pressure. The perturbation parameter αD is defined to be the natural logarithm of the ratio of the initial permeability to the permeability at the outflow boundary. The zeroth-order and first-order perturbation solutions are computed, from which the flux at the outflow boundary is found. An effective permeability is then determined such that, when inserted into the analytical solution for the mathematically linear problem, it yields a flux that is exact to at least first order in αD. When compared to numerical solutions of the problem, the result has 5% accuracy out to values of αD of about 2—a much larger range of accuracy than is usually achieved in similar problems. Finally, an explanation is given of why the change of variables proposed by Kikani and Pedrosa, which leads to highly accurate zeroth-order perturbation solutions in radial flow problems, does not yield an accurate result for one-dimensional flow.

Bird RE, Paluszny A, Zimmerman RW, 2021, Application of the Displacement Discontinuity Model to a finite pulse wave impinging on a fracture

For a harmonic wave of single frequency impinging on a fracture, the Pyrak-Nolte displacement discontinuity model (DDM) can be used to determine the reflection and transmission coefficients of the resultant harmonic waves. However, it is often the case in the field that a wave impinging on a fracture will be a pulse of finite duration which contains a spectrum of frequencies. Here, for a P-wave impinging a fracture at normal incidence, the DDM is extended to consider arbitrarily shaped pulse waves using the discrete Fourier Transform. The method presented calculates two resultant spectra for the transmitted and reflected wave as well as their corresponding energy coefficients. The impinging pulse wave is made dimensionless with respect to a defined "characteristic"frequency; this allows for a dimensionless comparison to the energy response of a single frequency harmonic wave using the Pyrak-Nolte curves. The results show that the total reflected energy of a finite-pulse incoming wave will be lower than that of a harmonic wave having the same characteristic frequency.

Salimzadeh S, Zimmerman RW, Khalili N, 2020, Gravity hydraulic fracturing: a method to create self-driven fractures, *Geophysical Research Letters*, Vol: 47, ISSN: 0094-8276

In this study, we investigate the possibility of using a high‐density fluid to induce downward fracture growth in a hydraulic fracturing process. We propose a mathematical model to calculate the minimum amount of a dense fluid required to trigger downward fracture propagation under gravity forces, and we verify the calculated minimum volume of the fluid through numerical simulations. Results show that when the injected fluid exceeds the minimum amount, a steady downward growth of the hydraulic fracture is obtained. The fracture propagation consists of two distinct responses: The first response can occur under either toughness‐dominated, viscosity‐dominated, or an intermediate hydraulic fracturing regime, depending on fluid rheology, rock properties, and injection scenario. The second response occurs mainly under the toughness‐dominated regime, meaning the predominant energy dissipation mechanism is to overcome the fracture toughness and break the rock. In the latter, the speed of the downward fracture growth depends on the viscosity and fluid weight.

Setiawan NB, Zimmerman RW, 2020, A unified methodology for computing the stresses around an arbitrarily-shaped hole in isotropic or anisotropic materials, *International Journal of Solids and Structures*, Vol: 199, Pages: 131-143, ISSN: 0020-7683

A unified semi-analytical solution based on graphical conformal mapping and complex variable methods is proposed to calculate the in-plane stress around an arbitrarily-shaped hole in isotropic or anisotropic materials. The method requires only the outline coordinates of the hole, the elastic moduli of the material, and the magnitude and direction of the far-field stresses. Comparison with many published results for a wide range of shapes, such as triangles, squares, ovaloids, and ellipses, has been carried out to validate the method. The method has also been applied to a highly irregular geometry that has been observed in a breakout of a subsurface borehole. The solution is essentially closed-form, in the sense that it can be explicitly expressed in terms of the mapping coefficients, and parameters that depend only on the elastic moduli of the materials. With such a degree of flexibility, the method will be useful to study the effect of hole geometry on the stress distribution around holes in isotropic or anisotropic materials.

Paluszny A, Thomas RN, Saceanu MC,
et al., 2020, Hydro-mechanical interaction effects and channelling in three-dimensional fracture networks undergoing growth and nucleation, *JOURNAL OF ROCK MECHANICS AND GEOTECHNICAL ENGINEERING*, Vol: 12, Pages: 707-719, ISSN: 1674-7755

The flow properties of geomechanically generated discrete fracture networks are examined in the context of channelling. Fracture networks are generated by growing fractures in tension, modelling the low permeability rock as a linear elastic material. Fractures are modelled as discrete surfaces which grow quasi-statically within a three-dimensional (3D) volume. Fractures may have their locations specified as a simulation input, or be generated as a function of damage, quantified using the local variation in equivalent strain. The properties of the grown networks are shown to be a product of in situ stress, relative orientation of initial flaws, and competitive process of fracture interaction and growth. Fractures grow preferentially in the direction perpendicular to the direction of maximum tension and may deviate from this path due to mechanical fracture interaction. Flow is significantly channelled through a subset of the fractures in the full domain, consistent with observations of other real and simulated fractures. As the fracture networks grow, small changes in the geometry of the fractures lead to large changes in the locations and scale of primary flow channels. The flow variability and formation of channels are examined for two growing networks, one with a fixed amount of fractures, and another with nucleating fractures. The interaction between fractures is shown to modify the local stress field, and in turn the aperture of the fractures. Pathways for single-phase flow are the results of hydro-mechanical effects in fracture networks during growth. These are the results of changes to the topology of the network as well as the result of mechanical self-organisation which occurs during interaction leading to growth and intersection.

Paluszny Rodriguez A, Graham CC, Daniels K,
et al., 2020, Caprock integrity and public perception studies of carbon storage in depleted hydrocarbon reservoirs, *International Journal of Greenhouse Gas Control*, Vol: 98, ISSN: 1750-5836

Capture and subsurface storage of CO2 is widely viewed as being a necessary component of any strategy to minimise and control the continued increase in average global temperatures. Existing oil and gas reservoirs can be re-used for carbon storage, providing a substantial fraction of the vast amounts of subsurface storage space that will be required for the implementation of carbon storage at an industrial scale. Carbon capture and storage (CCS) in depleted reservoirs aims to ensure subsurface containment, both to satisfy safety considerations, and to provide confidence that the containment will continue over the necessary timescales. Other technical issues that need to be addressed include the risk of unintended subsurface events, such as induced seismicity. Minimisation of these risks is key to building confidence in CCS technology, both in relation to financing/liability, and the development and maintenance of public acceptance. These factors may be of particular importance with regard to CCS projects involving depleted hydrocarbon reservoirs, where the mechanical effects of production activities must also be considered. Given the importance of caprock behaviour in this context, several previously published geomechanical caprock studies of depleted hydrocarbon reservoirs are identified and reviewed, comprising experimental and numerical studies of fourteen CCS pilot sites in depleted hydrocarbon reservoirs, in seven countries (Algeria, Australia, Finland, France, Germany, Netherlands, Norway, UK). Particular emphasis is placed on the amount and types of data collected, the mathematical methods and codes used to conduct geomechanical analysis, and the relationship between geomechanical aspects and public perception. Sound geomechanical assessment, acting to help minimise operational and financial/liability risks, and the careful recognition of the impact of public perception are two key factors that can contribute to the development of a successful CCS project in a

Saceanu MC, Paluszny A, Zimmerman RW, et al., 2020, Numerical modelling of spalling around a nuclear waste storage deposition borehole using a fracture mechanics approach, 54th U.S. Rock Mechanics/Geomechanics Symposium

Fracture growth leading to mechanical spalling around a deposition borehole for disused nuclear fuel waste is modelled numerically. Simulations are conducted using a finite-element-based discrete fracture growth simulator, which computes deformation in the system based on the mechanical properties of the rock. Fractures are grown by computing stress intensity factors at each fracture tip, and the mesh is re-generated to accommodate the changing fracture geometries at every growth step. Several numerical models are created to explore the effect of boundary conditions on the initiation and development of spalling fractures at Forsmark, where the Swedish repository for nuclear waste is planned to be constructed. It is shown that the reported uncertainty in the principal stress magnitudes and orientations will affect the predicted fracture nucleation and growth patterns, and implicitly the final repository design. The potential effect of spalling on the structural integrity of the deposition borehole is illustrated for each stress scenario.

Thomas RN, Paluszny A, Zimmerman RW, 2020, Permeability of three‐dimensional numerically grown geomechanical discrete fracture networks with evolving geometry and mechanical apertures, *Journal of Geophysical Research: Solid Earth*, Vol: 125, ISSN: 2169-9313

Fracture networks significantly alter the mechanical and hydraulic properties of subsurface rocks. The mechanics of fracture propagation and interaction control network development. However, mechanical processes are not routinely incorporated into discrete fracture network (DFN) models. A finite element, linear elastic fracture mechanics‐based method is applied to the generation of three‐dimensional geomechanical discrete fracture networks (GDFNs). These networks grow quasi‐statically from a set of initial flaws in response to a remote uniaxial tensile stress. Fracture growth is handled using a stress intensity factor‐based approach, where extension is determined by the local variations in the three stress intensity factor modes along fracture tips. Mechanical interaction between fractures modifies growth patterns, resulting in nonuniform and nonplanar growth in dense networks. When fractures are close, stress concentration results in the reactivation of fractures that were initially inactive. Therefore, GDFNs provide realistic representations of subsurface networks that honor the physical process of concurrent fracture growth. Hydraulic properties of the grown networks are quantified by computing their equivalent permeability tensors at each growth step. Compared to two sets of stochastic DFNs, GDFNs with uniform fracture apertures are strongly anisotropic and have relatively higher permeabilities at high fracture intensities. In GDFN models, where fracture apertures are based on mechanical principles, fluid flow becomes strongly channeled along distinct flow paths. Fracture orientations and interactions significantly modify apertures, and in turn, the hydraulic properties of the network. GDFNs provide a new way of understanding subsurface networks, where fracture mechanics is the primary influence on their geometric and hydraulic properties.

Thomas R, Paluszny A, Zimmerman RW, 2020, Growth of three-dimensional fractures, arrays, and networks in brittle rocks under tension and compression, *Computers and Geotechnics*, Vol: 121, ISSN: 0266-352X

Concurrent growth of multiple fractures in brittle rock is a complex process due to mechanical interaction effects. Fractures can amplify or shield stress on other fracture tips, and stress field perturbations change continuously during fracture growth. A three-dimensional, finite-element based, quasi-static growth algorithm is validated for mixed mode fracture growth in linear elastic media, and is used to investigate concurrent fracture growth in arrays and networks. Growth is governed by fracture tip stress intensity factors, which quantify the energy contributing to fracture extension, and are validated against analytical solutions for fractures under compression and tension, demonstrating that growth is accurate even in coarsely meshed domains. Isolated fracture geometries are compared to wing cracks grown in experiments on brittle media. A novel formulation of a Paris-type extension criterion is introduced to handle concurrent fracture growth. Fracture and volume-based growth rate exponents are shown to modify fracture interaction patterns. A geomechanical discrete fracture network is generated and examined during its growth, whose properties are the direct result of the imposed anisotropic stress field and mutual fracture interaction. Two-dimensional cut-plane views of the network demonstrate how fractures would appear in outcrops, and show the variability in fracture traces arising during interaction and growth.

Setiawan NB, Zimmerman RW, 2020, Semi-analytical method for tracking the evolution of borehole breakouts

In response to the high-stress concentration caused by excavation, rocks at the wellbore wall will collapse if the shear failure limit is exceeded. When this occurs, the shape of the wellbore enlarges from circular to roughly elliptical (breakout). The standard stress calculation method assumes a wellbore with a circular cross-section, and hence can no longer be used. Yet, knowledge of the near-wellbore stress state after the wellbore has been damage may still be required for various purposes, such as stability assessment of the open hole. In this paper, a semi-analytical solution is presented for the calculation of the stresses around a wellbore having an essentially arbitrary cross-sectional shape, in an isotropic formation. The zone of initial shear failure around a circular wellbore is calculated for a given far-field stress, using the standard solution for a circular wellbore. The rock in the failed region is then assumed to be removed, and the stress state around this new wellbore contour is then calculated using the method of conformal mapping and Kolosov-Muskhelishvili complex stress potentials. This process is iterated until a shape is obtained for which the breakout will not progress any further, at which point a stable stress state, and stable borehole shape, is reached. The method is validated against several sets of experimental data from the literature. This semi-analytical solution allows a stability assessment around the wellbore to be carried out after the breakout develops. Similarly, in an inverse manner, the method can be used to predict the in situ stress state around a damaged wellbore, if the cross-sectional hole shape can be reconstructed by methods such as electrical logging tools.

Almulhim OA, Paluszny A, Thomas RN, et al., 2019, Fully-coupled three-dimensional finite element simulations of the interaction between a hydraulic fracture and a pre-existing natural fracture

The effects of approach angle and matrix poroelasticity on three-dimensional fracture interactions are investigated by examining changes in stress intensity factors around hydraulic fracture tips. Additionally, the effect of the compressive stress induced by the hydraulic fracture (the “stress shadow effect”) on the natural fracture is explored. Fracture interaction is captured using three-dimensional interaction maps based on two interaction measures that quantify the magnitude and type of interaction. The results show that the stress conditions at the hydraulic fracture tip are more favorable for growth when interacting fractures have shallow approach angles, and the interaction reduces as the approach angle increases. For “high” permeability reservoirs, fluid leak-off causes the rock matrix to dilate, and generates a tensile stress that amplifies the local stress field. The increase in stress in the region ahead of the hydraulic fracture (the stress amplification zone) can cause natural fractures to open in tension before the hydraulic fracture intersects the natural fracture. The increase in matrix permeability causes earlier activation of natural fractures due to the additional stress induced by the dilated rock matrix.

Almulhim OA, Paluszny A, Thomas RN, et al., 2019, Fully-coupled three-dimensional finite element simulations of the interaction between a hydraulic fracture and a pre-existing natural fracture

Copyright 2019 ARMA, American Rock Mechanics Association. The effects of approach angle and matrix poroelasticity on three-dimensional fracture interactions are investigated by examining changes in stress intensity factors around hydraulic fracture tips. Additionally, the effect of the compressive stress induced by the hydraulic fracture (the “stress shadow effect”) on the natural fracture is explored. Fracture interaction is captured using three-dimensional interaction maps based on two interaction measures that quantify the magnitude and type of interaction. The results show that the stress conditions at the hydraulic fracture tip are more favorable for growth when interacting fractures have shallow approach angles, and the interaction reduces as the approach angle increases. For “high” permeability reservoirs, fluid leak-off causes the rock matrix to dilate, and generates a tensile stress that amplifies the local stress field. The increase in stress in the region ahead of the hydraulic fracture (the stress amplification zone) can cause natural fractures to open in tension before the hydraulic fracture intersects the natural fracture. The increase in matrix permeability causes earlier activation of natural fractures due to the additional stress induced by the dilated rock matrix.

Setiawan NB, Zimmerman RW, 2019, The implications of using anisotropic elasticity and fully-triaxial failure criteria for borehole stability analysis in shales

Setiawan and Zimmerman (IJRMMS, 2018) recently developed a computational toolkit for analyzing the stability of boreholes drilled at an arbitrary angle through an anisotropic formation. The borehole stresses can be computed using either the Hiramatsu-Oka solution based on isotropic elasticity, or the Lekhnitskii-Amadei anisotropic elasticity solution. Shear failure along the borehole wall can be modeled using two variants of the Jaeger plane of weakness model, using either the Mohr-Coulomb or the Mogi-Coulomb failure criterion for failure along planes other than a bedding plane. In this paper, the implications of using these various models are investigated. All four variants of the above-mentioned models have been evaluated over a wide range of the relevant parameter space, including the elastic anisotropy ratio, the well deviation angle, the in situ stress anisotropy, etc. The breakout pressure will be influenced by the elastic anisotropy when the ratio Ev/Eh exceeds about 2.3, even in a nearly-vertical wellbore. This effect is further amplified as the degree of elastic anisotropy increases, or the wellbore becomes more deviated. Due to its incorporation of the strengthening effect of the intermediate principal stress, the Mogi-Coulomb criterion predicts lower values of the minimum required mud weight than does the Mohr-Coulomb criterion.

Setiawan NB, Zimmerman RW, 2019, The implications of using anisotropic elasticity and fully-triaxial failure criteria for borehole stability analysis in shales

Copyright 2019 ARMA, American Rock Mechanics Association. Setiawan and Zimmerman (IJRMMS, 2018) recently developed a computational toolkit for analyzing the stability of boreholes drilled at an arbitrary angle through an anisotropic formation. The borehole stresses can be computed using either the Hiramatsu-Oka solution based on isotropic elasticity, or the Lekhnitskii-Amadei anisotropic elasticity solution. Shear failure along the borehole wall can be modeled using two variants of the Jaeger plane of weakness model, using either the Mohr-Coulomb or the Mogi-Coulomb failure criterion for failure along planes other than a bedding plane. In this paper, the implications of using these various models are investigated. All four variants of the above-mentioned models have been evaluated over a wide range of the relevant parameter space, including the elastic anisotropy ratio, the well deviation angle, the in situ stress anisotropy, etc. The breakout pressure will be influenced by the elastic anisotropy when the ratio Ev/Eh exceeds about 2.3, even in a nearly-vertical wellbore. This effect is further amplified as the degree of elastic anisotropy increases, or the wellbore becomes more deviated. Due to its incorporation of the strengthening effect of the intermediate principal stress, the Mogi-Coulomb criterion predicts lower values of the minimum required mud weight than does the Mohr-Coulomb criterion.

Setiawan NB, Zimmerman RW, 2018, Wellbore breakout prediction in transversely isotropic rocks using true-triaxial failure criteria, *International Journal of Rock Mechanics and Mining Sciences*, Vol: 112, Pages: 313-322, ISSN: 0020-7624

This paper presents a unified approach through which the influence of the elastic and strength anisotropy on wellbore instability can be thoroughly examined. The stresses at the wellbore wall are first calculated using the Lekhnitskii-Amadei solution, which accounts for elastic anisotropy. Then, shear failure is treated by combining the Mogi-Coulomb criterion for intact rock, with the Jaeger plane of weakness concept. The developed model accounts for all three principal stresses in predicting the onset of shear failure.The results of the specific case investigated show that rock elastic anisotropy induce higher stress concentrations. The difference, compared with the stresses found using the isotropic elastic model, could reach as high as 25% for the highest degree of anisotropy that might be expected for rocks of practical interest. The strengthening effect of the intermediate stress, as reflected in the Mogi-Coulomb criterion, reduces the required mud weight density by approximately 1.0 pounds-per-gallon (ppg). Furthermore, it is demonstrated that the risk posed by bedding slippage, for a wellbore with an inclination between 15° and 50° from the vertical, is masked when an isotropic elastic stress model is used. In contrast, the fully anisotropic model shows that an extra mud weight of approximately 4.5 ppg would be required, in order to avoid bedding plane slippage for the case under investigation. Although these results apply for a particular choice of strength properties and elastic properties, they give an indication of the implications of fully accounting for anisotropy and the effect of the intermediate stress when doing borehole stability analysis.

Thomas RN, Paluszny A, Hambley D,
et al., 2018, Permeability of observed three dimensional fracture networks in spent fuel pins, *Journal of Nuclear Materials*, Vol: 510, Pages: 613-622, ISSN: 0022-3115

The three-dimensional fracture network within a spent fuel pin is characterised using sequential grinding, and its permeability is numerically estimated. Advanced Gas-cooled Reactors (AGRs) produce spent fuel pins consisting of an outer steel cladding enclosing ceramic uranium dioxide (UO2) pellets. During irradiation, fuel pellets may undergo fracturing due to thermal, densification and swelling effects. Fracture patterns are usually observed on the surface of the pellet or through a cross section or longitudinal plane along the pellet. In this work, the three-dimensional fracture pattern within the pellet is characterised using an optical microscope. The pellet is progressively ground and polished, providing sequential cross sections, which together yield a three-dimensional discrete fracture pattern. Multiple large fractures grow to connect the cladding to the internal region of the pellet. Multiple surface fractures are observed that do not penetrate into the matrix of the pellet. The porosity of the UO2 and apertures of the fractures are estimated by sampling microscopic images. Darcy flow is numerically solved using the finite element method, computing flow through the matrix and fractures simultaneously. The equivalent tensorial permeability of the system is estimated for various approximate fracture apertures. The fracture network raises the permeability of the pellet by an order of magnitude.

Lang P, Paluszny Rodriguez A, Morteza N,
et al., 2018, Relationship between the orientation of maximum permeability and intermediate principal stress in fractured rocks, *Water Resources Research*, Vol: 54, Pages: 8734-8755, ISSN: 0043-1397

Flow and transport properties of fractured rock masses are a function of geometrical structures across many scales. These structures result from physical processes and states and are highly anisotropic in nature. Fracture surfaces often tend to be shifted with respect to each other, which is generally a result of stress‐induced displacements. This shift controls the fracture's transmissivity through the pore space that forms from the created mismatch between the surfaces. This transmissivity is anisotropic and greater in the direction perpendicular to the displacement. A contact mechanics‐based, first‐principle numerical approach is developed to investigate the effects that this shear‐induced transmissivity anisotropy has on the overall permeability of a fractured rock mass. Deformation of the rock and contact between fracture surfaces is computed in three dimensions at two scales. At the rock mass scale, fractures are treated as planar discontinuities along which displacements and tractions are resolved. Contact between the individual rough fracture surfaces is solved for each fracture at the small scale to find the stiffness and transmissivity that result from shear‐induced dilation and elastic compression. Results show that, given isotropic fracture networks, the direction of maximum permeability of a fractured rock mass tends to be aligned with the direction of the intermediate principal stress. This reflects the fact that fractures have the most pronounced slip in the plane of the maximum and minimum principal stresses, and for individual fractures transmissivity is most pronounced in the direction perpendicular to this slip.

Salimzadeh S, Paluszny Rodriguez A, Zimmerman RW, 2018, Effect of cold CO2 injection on fracture apertures and growth, *International Journal of Greenhouse Gas Control*, Vol: 74, Pages: 130-141, ISSN: 1750-5836

The injection of cold CO2 is modelled in three dimensions using a two-stage coupled thermoporoelastic model, with the aim of evaluating changes in apertures and potential growth of fractures. Non-isothermal flow is considered within the fractures and the rock matrix, and the two flow domains are coupled through a mass transfer term. The numerical model has been developed using standard finite elements, with spatial discretisation achieved using the Galerkin method, and temporal discretisation using finite differences. A full-scale field case geometric model, based on the Goldeneye depleted hydrocarbon reservoir in the North Sea, is developed and used for simulations. The in situ faults are modelled discretely as discontinuous surfaces in a three-dimensional matrix, including basement, reservoir, caprock and overburden layers. The faults are assumed initially to be low-permeable faults, with the same permeability as the caprock. However, the simulations show that their apertures (and as a result, their permeabilities) vary due to the thermoporoelastic effects caused by the injection of the relatively cold CO2. The change in the fracture apertures is mainly due to thermal effects; the reservoir layer undergoes contractions due to the cooling, significantly increasing fault aperture in the region of the fault within the reservoir, whereas the fault aperture is reduced in regions within the caprock. Propagation of fractures under thermoporoelastic loading is investigated. Results show that the distance to the injection well, as well as spatial orientation of fractures with respect to the injection well, affect aperture evolution and potential growth of fractures. A sensitivity analysis is performed on the parameters affecting the fracture growth: minimum normal stress acting on the fracture plane, dip angle of the fracture, and the contact friction coefficient. It is found that low friction, low normal contact stress, or high in situ shear stress on the fracture surface

This data is extracted from the Web of Science and reproduced under a licence from Thomson Reuters. You may not copy or re-distribute this data in whole or in part without the written consent of the Science business of Thomson Reuters.