University of Southern California Dissertations and Theses

Synergistic coupling between geomechanics, flow, and transport in fractured porous media: applications in hydraulic fracturing and fluid mixing

(USC Thesis Other)

Synergistic coupling between geomechanics, flow, and transport in fractured porous media: applications in hydraulic fracturing and fluid mixing

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content SYNERGISTIC COUPLING BETWEEN GEOMECHANICS, FLOW, AND TRANSPORT IN FRACTURED POROUS MEDIA: APPLICATIONS IN HYDRAULIC FRACTURING AND FLUID MIXING by Minh Tuan Tran A Dissertation Presented to the FACULTY OF THE USC GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulfillment of the Requirements for the Degree DOCTOR OF PHILOSOPHY (PETROLEUM ENGINEERING) December 2023 Copyright 2023 Minh Tuan Tran , "Energy is the golden thread that connects economic growth, social equity, and environmental sus- tainability." - Ban Ki-moon This dissertation is dedicated to my beloved family, whose steadfast love, support, and encouragement have been the cornerstone of my journey. Throughout the years, you have been my pillars of strength, providing me with the motivation and inspiration to pursue my academic aspirations. Your boundless belief in my abilities has been a constant source of reassurance, pushing me to surpass my own expectations. To my parents, who have been my guiding lights and the epitome of selflessness, I am eternally grateful. Your sacrifices, both big and small, have paved the way for my success. You have instilled in me the values of hard work, determination, and perseverance, teaching me that no dream is too far-fetched and no obstacle is insurmountable. To my wife, Van Bui, words cannot adequately express the depth of my gratitude. Your belief in me has propelled me forward and shaped me into the person I am today. Without your resolute love and support, this dissertation would not have been possible. Your unwavering support and unfaltering faith in me during the toughest periods of the Ph.D. has given me the courage to embark on this intellectual journey and strive for excellence. To my two beautiful angels, Ryan and Evelyn, may it serve as a testament to the power of educa- tion, hard work, and perseverance. May it inspire you to pursue your dreams relentlessly, with the knowledge that you, too, have a strong foundation of love and support behind you. ii ACKNOWLEDGEMENTS During my academic pursuit to reach this significant milestone, numerous individuals have made valuable intellectual contributions to this endeavor, and I am truly grateful for their support. Above all, I extend my heartfelt appreciation to my mentor, Dr. Birendra Jha, for providing me with an incredibly enriching academic atmosphere over the past five years. Professor Jha’s dedication to fostering intellectual growth and motivating his students to achieve academic excellence, in my view, is unmatched. I express my gratitude for his generous allocation of time and resources, which have played a crucial role throughout the entire duration of this project. I have been fortunate to have had numerous opportunities to engage in discussions about my research with the esteemed faculty members at the University of Southern California. Professors Iraj Ershaghi and Fred Aminzadeh, in particular, offered valuable critiques and encouraged close collaboration with their research groups. My time as a member of the Geosystems Engineering and Multiphysics (GEM) lab has been highly enjoyable, and I will cherish the memories of our weekly meetings and conference travels with Dr. Magdalene Ante, Dr. Saro Meguerdijian, Renuhaa Asaithambi, and Dr. Xiaoxi Zhao. Being affiliated with such a remarkable academic environment, alongside exceptional students and faculty, has been a privilege. I am deeply thankful to the members of my dissertation committee, namely Professor Kristian Jessen, Professor Yehuda Ben-Zion, and Professor Derek Elsworth, for their critical reviews of various aspects of my work. Their insights have been invaluable. I would also like to express my gratitude for the indispensable support provided by the staff of the Viterbi School. I am appreciative ofthefinancialbackingfromtheindustrialaffiliatesoftheGEMlab, andIacknowledgetheirsupport with gratitude. My financial assistance came through the Provost Fellowship from the USC Viterbi School of Engineering. Lastly, I conducted numerical simulations using the computing resources at the High Performance Computing and Simulations, which played a crucial role in my research. I want to express my heartfelt gratitude to my extended family, relatives, and friends who have consistently provided me with encouragement and motivation. Your love and unwavering belief iii in my capabilities have been a constant source of strength. Even during moments of doubt and exhaustion, your steadfast support, understanding, and uplifting words have propelled me forward. Your presence in my life has been a reminder of the profound significance of human connections, and I am eternally thankful for your enduring faith in me. iv TABLE OF CONTENTS DEDICATION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii ACKNOWLEDGEMENTS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii LIST OF TABLES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii ABSTRACT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xvii CHAPTER 1 : Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Motivation and Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Overview of Chapters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 CHAPTER 2 : CouplingbetweenTransportandGeomechanicsAffectsSoluteSpreadingand Mixing during Viscous Fingering in Deformable Aquifers . . . . . . . . . . . 4 2.1 Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 2.2 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 2.3 Physical Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.4 Mathematical Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.5 Numerical Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 2.6 Numerical Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 2.7 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 2.8 Sensitivity Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 2.9 Summary and Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 CHAPTER 3 : Effect of Poroelastic Coupling and Fracture Dynamics on Solute Transport and Geomechanical Stability . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 3.1 Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 v 3.2 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 3.3 Physical Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 3.4 Mathematical Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 3.5 Numerical Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 3.6 Model Verification and Benchmarking . . . . . . . . . . . . . . . . . . . . . . . . . . 58 3.7 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 3.8 Sensitivity Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 3.9 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 CHAPTER 4 : An Embedded Discrete Fracture–Phase Field Model to Quantify Fracturing- induced Pore Pressure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 4.1 Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 4.2 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 4.3 Mathematical Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 4.4 Numerical Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 4.5 Model Verification and Benchmarking . . . . . . . . . . . . . . . . . . . . . . . . . . 107 4.6 Stress Shadow and Fracture Interaction . . . . . . . . . . . . . . . . . . . . . . . . . 113 4.7 Field Case Study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123 4.8 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128 CHAPTER 5 : Key Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130 BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132 vi LIST OF TABLES TABLE 2.1 Simulation parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 TABLE 2.2 Base Case parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 TABLE 2.3 Sensitivity cases for the permeability heterogeneity. . . . . . . . . . . . . . . . . 32 TABLE 3.1 Groundwater chemicals with viscosity µ 1 different from the water viscosity µ 0 . A water viscosity of 0.89 cp at 25 degC is assumed. . . . . . . . . . . . . . . . . 41 TABLE 3.2 Simulation parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 TABLE 3.3 Base Case parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 TABLE 3.4 Sensitivity parameter values. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 TABLE 4.1 Tension Test Parameters. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 TABLE 4.2 2D Fluid-filled fracture propagation. . . . . . . . . . . . . . . . . . . . . . . . . 111 TABLE 4.3 Model parameters for investigating stress shadow and fracture interaction . . . 114 TABLE 4.4 Sensitivity parameter values. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 TABLE 4.5 Field case study input parameters . . . . . . . . . . . . . . . . . . . . . . . . . . 124 vii LIST OF FIGURES FIGURE 2.1 A physical model of alternating slug injection into a reservoir under a tectonic state of stress created by compressions σ h and σ H on its two sides. Light- colored (viscosity µ 1 ) and dark-colored (viscosity µ 0 >µ 1 ) slugs are injected alternatingly into the reservoir from an injection well, and the fluids are pro- duced from the producer. Displacement of the more viscous fluid by the less viscous fluid leads to viscous fingering. Fingering-induced perturbation in the pressure field generates perturbation in the effective stress field, which affects deformation and hydraulic properties of the reservoir, which in turn can affect spreading of the less viscous fluid and production of the more viscous fluid. . 8 FIGURE 2.2 A Gaussian random permeability field. Values of the Gaussian model parame- ters are as follows: geometric mean k g =50 md, log-k variance σ 2 logk =0.01, and correlation lengths ξ x =ξ y =1 m. The color scale is in millidarcy unit. . . . . . 10 FIGURE 2.3 Base Case: (a) Concentration evolution is shown through snapshots in time. I see the development of fingers and classical viscous fingering mechanisms— tip-splitting, shielding, and merging of fingers. (b) Transversely-averaged con- centration c and change in pressure ∆ p (from t 0 ) are plotted along the flow direction at three time steps. Fingers spread the tracer downstream of the main slug as shown by c curves at t 1 and t br . ∆ p goes up and down with changing mobility of the fluid in the system; the yellow fluid is more mobile than the blue fluid. ∆ p at the producer is always zero because of its fixed pressure boundary condition. . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 FIGURE 2.4 Base Case results at the time of tracer breakthrough at the producer. The con- centration field c shows viscous fingers. Magnitudes of the pressure gradient |∇p|, the effective volumetric stress gradient |∇σ ′ v |, and the porosity gradient |∇ϕ | fields are shown. The white curves are the transversely averaged profiles of the 2D fields, scaled to fit within the domain width of [0,L x ]. The similar- ity between the average concentration profile and the average gradient profiles confirms that viscous fingering affects, and is affected by, the gradient of pres- sure, stress and porosity fields. The time rate of change in porosity, dϕ/dt , also shows a signature of the fingered region. . . . . . . . . . . . . . . . . . . 26 FIGURE 2.5 Base Case results at the observation point at the center of the domain. Change in concentration ∆ c causes, and is coupled to, the changes in pressure ∆ p, effective volumetric stress ∆ σ ′ v , von Mises stress ∆ σ vm , porosity ∆ ϕ , and de- formations ∆ U x and ∆ U y . t is the dimensionless time. Top: features of the rate of change of porosity, dϕ/dt , are related to the concentration evolution, maximum magnitudes of U x and U y , and tracer breakthrough which occurs at t br . The three vertical lines are time markers that relate features in the dϕ/dt curve to breakthrough and maximum deformation events. . . . . . . . . . . . 26 viii FIGURE 2.6 Spreading and mixing in the Base Case. (a) Degree of mixing χ evolves non- monotonically because of the balance between source and sink terms of the mixing process. The decline in χ begins shortly after breakthrough when the sink begins to dominate source. (b) The mean scalar dissipation rate ϵ c curve jumps when the trailing front enters the domain (source of concentration gradi- ents), and it declines rapidly aftert mc because viscous fingering has ceased and slug drains out (sink). The blue dashed curve is the analytical ϵ c (Eq. 2.23) for Fickian mixing in a periodic (no entry or breakthrough of slugs), homogeneous, viscously stable and geomechanically-uncoupled flow. The difference between analytical and numerical ϵ c are due to fingering-enhanced interface stretching during 0<t<0.7, entry of the slug’s trailing front at t≈ 0.7, and breakthrough at t br =1.16. (c) Wiggles in the breakthrough curve, c prd , indicate arrivals of fingers at the producer. t mc is the time at which c prd is maximum. . . . . . . 28 FIGURE 2.7 Coupling between transport and deformation. (a) The degree of mixing χ and the root mean square of porosity change, ∆ ϕ rms , increase and decrease synchronously as the slug is transported through the deforming reservoir. (b) The breakthrough concentration c prd and the root mean square of effective volumetric stress change,∆ σ ′ v,rms , are also related. Both relations are sensitive to the values of the Biot coefficient α and the drained modulus K. . . . . . . 29 FIGURE 2.8 Sensitivity of the degree of mixing at breakthrough (χ br ) and the breakthrough time (t br ) to deformation parameters (upper row) and transport parameters (lower row). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 FIGURE 2.9 Effect of Biot’s coefficient α on spreading. The breakthrough concentration c prd vs. the dimensionless time for different values of α shows that a less consolidated (higher α ) formation has an earlier arrival of the maximum con- centration (shorter t mc ) compared to a more unconsolidated one. c prd (t mc ) have similar values for different values of α . . . . . . . . . . . . . . . . . . . . . 31 FIGURE 2.10 Effect of Biot’s coefficient α , a measure of poroelastic coupling strength, on the tracer concentration field. The concentration field corresponding to α = 0.25 (relatively low poroelastic coupling strength) is subtracted from the concentra- tion field corresponding to α =0.5 (left), 0.7 (middle), and 0.9 (right) at their respective breakthrough times. As the poroelastic coupling strength increases, the difference in the concentration field also increases. See Table 2.2 for all other parameter values. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 FIGURE 2.11 (a)Heterogeneouspermeabilityfieldsfordifferentvaluesofthelogpermeability variance σ 2 logk . Color scale is log(k) with k in md. (b) Concentration fields at breakthrough for the three cases. Finger tips are more pointy (less rounded), split and diffused, and fingers are longer at highest σ 2 logk . Color scale is for dimensionless concentration. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 FIGURE 2.12 (a) Heterogeneous permeability fields for different values of the correlation length ratio ξ y /ξ x . (b) Concentration at breakthrough for the three cases. Fingers are longer at larger ξ y /ξ x . . . . . . . . . . . . . . . . . . . . . . . . . 33 ix FIGURE 2.13 Concentration fields at breakthrough for different values of the Peclet number. Pe=700 shows more tip-splitting and less diffusion (sharper fingers) compared to Pe=200 and 400. Compared to Pe=200, the trailing front of the slug in Pe=400 is farther from the leading front, which means more spreading and faster breakthrough. Also, Pe=400 shows more fingers compared to Pe=200. . 34 FIGURE 2.14 Concentration fields for different cases of the log viscosity ratio R at the respec- tive breakthrough times. The breakthrough, which is indicated by the position of the trailing front, occurs earlier at higher R because of faster traveling fingers. 34 FIGURE 3.1 A physical model of tracer injection, tracer soaking or chase fluid injection, and tracer-waterwithdrawalinareservoirpopulatedwithdiscretenaturalfractures of different lengths and orientations. The horizontal well is intersected by some of the fractures. The initial state of stress in the reservoir is determined by boundary compressions S h and S H and zero normal displacement roller boundary conditions. The fractures are denoted by the green lines. The cyclic well is denoted in the thick blue line. The yellow-shaded region around the well indicates the mixture of injected tracer (viscosity µ 1 ) and resident water (viscosity µ 0 ), which forms a short time after tracer injection begins at the well. The remaining black domain is saturated with the resident fluid. . . . . 44 FIGURE 3.2 The modified Bandis model to model nonlinear fracture mechanics in response to dynamically evolving normal and shear stresses Li et al. (2017). The top left quadrant shows the behavior of a fracture in tension withK ni denoting the initial normal fracture stiffness. The bottom right quadrant shows a fracture in compression, displaying the initial mechanical aperturea 0 under zero excess ef- fective normal stress, the maximum aperture of the fracture under compression (V m ), and the change in mechanical aperture ∆ a n due to normal stress pertur- bation. Hysteresis within a loading-unloading cycle on the fracture is shown with the triangle marker (loading or withdrawal period) and the circle marker (unloading or injection period) lines. The upper right quadrant shows the de- pendence of shear stress (τ ) on shear displacement (δ h ). The mobilization of JRC is divided into the pre-peak and post-peak periods, separated by point B corresponding to peak shear displacement δ peak and peak shear strength τ peak . As the shear displacement increases to the residual shear displacement δ res , the residual shear strength τ res is slowly reached at points C and D (upper right quadrant). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 FIGURE 3.3 Flow and geomechanics are coupled two-way through strain-dependent poros- ity ϕ (ϵ v ) in the fluid mass balance equation (Equation 3.1) and the fluid pres- sure dependent effective stress σ ′ (p) in the linear momentum balance equation (Equation 3.13). The flow and transport are coupled two-way through the concentration-dependent fluid viscosity µ f (c) and pressure-dependent veloc- ity field v(p). The geomechanics is coupled one-way with fracture mechanics by providing the nodal displacement, stress, and strain tensors to the Bandis model for slip activation, slip magnitude calculation, and fracture aperture up- date. The fracture mechanics delivers updated stress-dependent permeability k f (δ h ,σ ′ n ) to the flow solver based on shear slip δ h and normal dilation, which is dependent on the effective normal stress σ ′ n . . . . . . . . . . . . . . . . . . 56 x FIGURE 3.4 (a) The Mandel problem setup with the origin of the x-y coordinate system positioned at the center of a plane strain 2D slab of dimension 2L x × 2L y . Symmetry allows modeling of a quarter of the slab, producing a computational domain of size L x × L y . The system is under a compressive vertical force de- noted as 2F applied instantaneously at t=0. (b) The dimensionless pressure p/p 0 vs. the dimensionless x-coordinate x d = x/L x at different dimensionless times t d . The lines show the simulation results and the dots are the analytical solutions. The simulation agrees with the analytical result and captures the Mandel-Cryer effect at the early time step of t d =0.05. . . . . . . . . . . . . . 60 FIGURE 3.5 Top row: Tracer concentrations in fractures and in the matrix at three suc- cessive time steps of the base case model: (a) at tracer breakthrough (t BT ) which is also the end of injection period, (b) at the end of the chase fluid in- jection period (soaking), and (c) 30 days into the withdrawal period. The well is denoted by the blue horizontal line. Matrix concentration c m is shown via contours drawn at an increment of 0.1. Fracture concentration c f is shown via colored straight lines drawn along the fracture length. Both follow the same color scale shown on right. The shape and size of the contours are controlled by the connectivity of the well to the fractures and the transport processes. Bottom row: corresponding excess pore pressure distribution at the three se- lected time steps. Color scales of the bottom panels are different, implying the sign reversal of excess pore pressure during push and pull periods. . . . . . . 63 FIGURE 3.6 Evolution of the tracer concentration at selected points along the well. The curves at the right and left ends of the well are similar during the injection periodwhilethemiddlepointshowsahigherconcentrationduetoahigherfrac- ture density near the center. The soaking period does not show any distinction among the three curves because the tracer is neither injected nor withdrawn. The withdrawal period curves are more separated from each other compared to the injection period because both advection and diffusion processes feel the fracture morphology deeper in the aquifer during later periods. . . . . . . . . 65 FIGURE 3.7 The temporal evolution of excess pore pressure (∆ p) in the matrix, excess effective volumetric stress ( ∆ σ ′ v ), and excess Mohr-Coulomb stress Invariant (∆ I ′ MC ) at different locations in the domain, which are marked in the inset figures in the bottom left. The locations are (a) near the midpoint of the well, (b) center of the upper half of the domain, and (c) near the top boundary. Pressure (star marker), Mohr-Coulomb stress invariant (black marker), Von Mises stress (blue marker), and effective volumetric stress (diamond marker) evolve with time as the well goes through tracer injection, soaking (or chase fluid injection), and withdrawal episodes. . . . . . . . . . . . . . . . . . . . . . 68 FIGURE 3.8 Change in direction, β S H (degrees), of the maximum horizontal stress S H , measured from its initial direction along the y-axis.A positive angle means the principal stress rotates counter-clockwise while a negative angle indicates clockwise rotation. The dashed lines indicate the new S H direction. Green thick lines are the embedded fracture network. The thick blue line is the well trajectory. Two time steps are shown: (a) at the end of injection, and (b) at the end of withdrawal. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 xi FIGURE 3.9 Empirical cumulative distribution function (ECDF) of the fracture permeabil- ityk f after initialization (square marker), at the end of injection (star marker), and at the end of withdrawal (diamond marker). Fracture permeabilities in- crease during injection and decrease during withdrawal. The two inset figures show the spatial distributions ofk f at the respective timesteps: injection in the right inset and withdrawal in the left inset. The insets use different color scales to highlight the different fractures activated during different episodes–injection and withdrawal. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 FIGURE 3.10 Thetemporalevolutionoffracturepermeabilityanditsdependenceonfracture orientation θ f . The fracture orientation θ f (degrees) is defined as the angle between the fracture plane and the initial S H direction (i.e. true North). It ranges from− 90 ◦ (clockwise from true North) to 90 ◦ (counter-clockwise from true North). The true North is the direction of S H in Figure 3.1. Panel (a) plots the permeability of three representative fracture segments that slipped andexperiencedshear-induceddilationduringtheinjectionperiod. Thelegend showsdifferentorientationsofeachfracturesegments. Thehorizontaltimeaxis is plotted on log scale to illustrate the early injection period better. The mini- panel shows the locations of the fracture segments that are plotted. Panel (b) plots fracture permeability and orientation for all fracture segments where circle markers represent unslipped segments and square markers denote slipped segments. The points are color-coded by the absolute magnitude of excess normal effective stress acting on the fracture segments ∆ σ ′,f n . . . . . . . . . . 71 FIGURE 3.11 Tornado plot of the percentage change of breakthrough time t BT (panel (a)) and degree of mixing ζ BT (panel (b)) compared to the base case of different sensitivity parameters listed on the y-axis. From bottom to top on the vertical axis: Biot Coefficient α , Fluid Compressibilityc f , Young ModulusE, Log Vis- cosity Ratio R, Joint Roughness Coefficient JRC, Joint Compressive Strength JCS, Diffusion Coefficient D d . Red columns to the left show reduction while blue columns to the right show increase of either breakthrough time or de- gree of mixing compared to the base case. The end values in each sensitivity bar indicate the corresponding values of each sensitivity parameter, with units shown in Table 3.4. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 FIGURE 3.12 Effect of Biot coefficient α on the tracer concentration distribution. The thick horizontal line is the cyclic well and the thin black lines are the fracture planes. Difference between two concentration fields for (a) α = 0.5 and α = 0.2, and (b) α = 0.9 and α = 0.2 are shown at the breakthrough time. The loss of tracer concentration is intensified as the flow-geomechanical coupling strength increases with increasing α . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 xii FIGURE 3.13 (a) Tracer retrieval curve data from a case study Li et al. (2016). The inset figure shows the location of the well (thick blue line), fracture (thinner green lines), and tracer monitoring stations (midpoint of well and midpoint of the right well segment). The horizontal axis is the time elapsed from the start of the withdrawal period, ∆ t wd , in days. (b) and (c) show tracer retrieval curves during tracer withdrawal periods at different well locations for different values of the fluid compressibility C f (b) and rock Modulus E (c). Lower panels demonstrate the effect of flow-geomechanical coupling (different colors but same line styles) and fracture morphology (same colors but different line styles). Allcurvesareshiftedtowardacommonstartingpointonthehorizontal axis for better visual comparison. . . . . . . . . . . . . . . . . . . . . . . . . . 77 FIGURE 3.14 Empirical cumulative distribution function (ECDF) of fracture permeability (k f ) plotted for different values of the Biot coefficient α (a), fluid compress- ibility C f (b), and rock’s modulus E (c). Injection and withdrawal responses are separated by plotting them as solid and dash lines, respectively. . . . . . . 79 FIGURE 3.15 Effect of the log viscosity ratio R, a transport parameter, on the stress state. (a) the root mean square of the square root of the second invariant of devi- atoric stress tensor p J 2,rms is plotted against the degree of mixing ζ . (b) The root mean square of Mohr-Coulomb stress invariant I ′ MC,rms is plotted against the domain-averaged concentration c m . Evolution during the injection and withdrawal periods are shown via solid and dotted lines, respectively. The direction of time for each period is shown via arrows. The shaded regions are magnified in (c) and (d) for better visualization. . . . . . . . . . . . . . . . . 81 FIGURE 3.16 Evolution in the stress states (shear stress τ versus effective normal stress σ ′ n ) of three representative fracture segments is shown using a Mohr circle diagram. Three timesteps are shown–initial (black square), end of injection (solid circles), and end of withdrawal (asterisks)–for three different JCS values. The static failure envelope is plotted in red with the static frictional coefficient of0.35. Theshadedregioninthetopleftdiagramisshownintherightdiagram with a zoom-in. The arrows indicate the direction of change in the stress, from the initial state, during injection and withdrawal. . . . . . . . . . . . . . . . . 83 FIGURE 3.17 Stress states of three representative fracture segments are plotted in a Mohr diagramatdifferenttimes: initial(blacksquare), endofinjection(solidcircles), and end of withdrawal (asterisks) for different JRC values. The static failure envelope is plotted in red with the static frictional coefficient of 0.35. The shaded region in the top left diagram is shown in the right diagram with a zoom-in. The arrows indicate the directional change of stress state from the initial state during injection and withdrawal. . . . . . . . . . . . . . . . . . . 83 xiii FIGURE 4.1 Flow and geomechanics are two-way coupled through strain-dependent poros- ity ϕ (ϵ v ) and the stress-dependent mass transfer function ψ in the fluid mass balance equation and the pressure-dependent effective stress σ ′ (p) in the linear momentum balance equation. The flow and fracture mechanics modules are coupled through the phase-field-dependent geomechanical parameters M(φ) and α (φ) and the pressure-dependent strain energy functional η + (p). In ad- dition, the fracture mechanics module delivers stress-dependent permeability k f (δ h ,σ ′ n ) to the flow module based on shear slip δ h and normal dilation, which is dependent on the effective normal stress σ ′ n . The geomechanics module is coupled to fracture mechanics by providing the displacement, stress, and strain to the Bandis model for slip activation, slip magnitude calculation, and frac- ture aperture update, as well as evaluating the strain energy functional η + (u). ThefracturemechanicsmoduleupdatestheBiotcoefficient α (φ)andthestress degradation function g(φ) in the geomechanics module. . . . . . . . . . . . . 108 FIGURE 4.2 Process diagram of the coupled flow-geomechanics-fracture propagation frame- work. In the initialization step, the initial local strain history functional is specified as well as the initial pressure/displacement and permeability field. Inside the time loop, there is an inner loop to iterate through the solutions of the two solvers, flow-geomechanics solver and the phase field solver, to a specified error tolerance tol. Poroelastic parameters are updated for use in the next time step based on the converged solution of the phase field in the inner loop. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 FIGURE 4.3 Left: A square plate with a straight edge crack subjected to vertical pull on the top boundary to simulate a pure mode-I fracture. The y-displacement on the top boundary is prescribed and increases monotonically with time. Right: Resultant force at the bottom boundary increases with the prescribed displace- ment until the onset of failure, which is marked with a sharp drop in the force. The proposed model agree with the literature. . . . . . . . . . . . . . . . . . . 110 FIGURE 4.4 Model setup to test the proposed framework on the fluid-filled fracture prop- agation problem. A hydraulic fracture propagating in an impermeable elastic medium along the direction orthogonal to the orientation of minimum com- pressive far-field stress σ . Fracture propagation is driven by injection of an incompressible, Newtonian fluid at a constant volumetric rate q w Detournay (2004). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111 FIGURE 4.5 From top to bottom: Spatial distributions of the phase field variable, fluid pressure, andhorizontalandverticaldisplacementsatthreeselectedtimesteps. Weobservealargepressuredifferencewithlowpressures(inbluecolor)infront of the fracture tips. Displacements are zero at the tips of the fracture, and the largest displacement is observed in the middle, a typical behavior of crack opening. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112 FIGURE 4.6 Evolution of the fracture pressure and fracture length with time. The re- sults are compared with the analytical solution Detournay (2004). The slight mismatch stems from the difference in boundary condition which affects the solution when the fracture approaches the edge of the domain. . . . . . . . . 113 xiv FIGURE 4.7 Propagation of a hydraulically stimulated fracture (Stim. Frac.) monitored via a static observation fracture (Obs. Frac.). Spatial distribution of the phase fieldvariableandthefluidpressureatthreeselectedtimesteps(alongcolumns) are shown. The thin zig-zag line inside the Stim. Frac. in the bottom row panels show the fracture propagation trajectory determined by my fracture tracking algorithm in section 4.3.5. . . . . . . . . . . . . . . . . . . . . . . . . 115 FIGURE 4.8 Snapshots of the spatial distribution of the phase field at five selected time steps for different approach angles β . The rows are for different approach angles which change the fracture orientation from horizontal to vertical in the domain. The columns from left to right indicate the evolution in time. . . . . 116 FIGURE 4.9 Effect of the approach angle β of the stimulated fracture to the observation fracture. (a) As the stimulated fracture’s orientation with respect to the ob- servation fracture changes from parallel (β = 0) to perpendicular (β = 90), the stimulated fracture pressure gets smaller on average. For any given value of the angle, the pressure first drops and then rises as the fracture grows to sense boundaries in the domain. (b) Observation fracture pressure (normal- ized) versus the normalized time for different approach angles. . . . . . . . . . 117 FIGURE 4.10 Effect of the fracture approach angle β on stress trajectories in the invariant space of von Mises stress versus effective volumetric stress. Stress trajectory at a point also depends on the point’s position relative to the stimulated fracture. Left: Point A is at the center of the domain and falls on the path of the β = 90 ◦ stim. frac. (stimulated fracture becoming perpendicular to the obs. frac.). Since stimulated fractures are predominantly tensile-mode failures, the effective volumetric stress is tensile (positive) for β = 90 ◦ and remains tensile even after the tip moves past point A. For other values of β , the effective volumetric stress is compressive (negative), which is the origin of the stress shadow effect. Right: Point B falls closest to the path of the β = 30 ◦ stim. frac. As the fracture tip approaches point B, the effective volumetric stress becomes tensile. Further growth of the stimulated fracture causes compression at point B (turning of the trajectory) due to the stress shadow effect. . . . . . 118 FIGURE 4.11 Spatial distribution of the total volumetric stress at five selected time steps for different values of the fracture approach angle β . . . . . . . . . . . . . . . . . 119 FIGURE 4.12 Spatial distribution of the von Mises stress at five selected time steps for dif- ferent values of the approach angle β . Deviatoric stresses are affected by the approach angle. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120 FIGURE 4.13 Evolution of pore pressure at five selected time steps (along columns) for three values of the fracture spacing parameter (along rows). . . . . . . . . . . . . . 122 FIGURE 4.14 Observation fracture pressure (normalized) versus the normalized time for three different fracture spacing values. Fracture spacing is the distance be- tween the stimulated fracture and the observation fracture. At closer spacing, the stimulated fracture has lower pressure and shorter length because of the pressurization of the observation fracture. . . . . . . . . . . . . . . . . . . . . 122 xv FIGURE 4.15 Effect of change in mechanical properties (Young’s modulus E and Biot’s co- efficient α ) of the rock surrounding the stimulated fracture region (colored box on the left). (a) Observation fracture pressure, and (b) length of the stimulated fracture. A 50% reduction inα causes almost 54% reduction in the observation fracture pressure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123 FIGURE 4.16 Schematic of the field case study. W1, W2, and W3 are three horizontal wells in a horizontal section of a reservoir rock. Fractures are shown by oval shapes placed on the wells. Well W2 is hydraulically stimulated to create fracture stages 1 and 2 while pressure is monitored at wells W1 and W3. W1 and W3 contain one hydraulic fracture each from a prior stimulation job. . . . . . . . 125 FIGURE 4.17 Field case study data: Pressure and rate data during hydraulic fracturing of well W2 and its monitoring at wells W1 and W3. Two stages of fracturing are shown with similar injection rates for both stages (approximately 11.13 cubic meter per minute) and similar pressures (29.65 MPa). Stage 1 injection lasts about 30 minutes while stage 2 injection lasts about 20 minutes. . . . . . . . 126 FIGURE 4.18 Field case study: Spatial distribution of fracture pressure and the phase field variable at four selected time steps during stage-1 and stage-2 fracturing at well W2. Left side panels show that injection-induced overpressure in stage- 1 (red color) causes pressure diffusion and stress shadow that affects stage-2’s trajectoryandlength. Rightsidepanelsshowthatthestage-1fracturegrowsin a direction that is predominantly parallel to the observation fractures. Stage-2 fracture senses the presence of relatively lower compression inside stage-1 and higher compression in the region between the two stages, which dictates stage- 2’s growth in the orthogonal direction and the intersection of the two fractures. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126 FIGURE 4.19 Fieldcasestudy: Spatialdistributionofeffectivevolumetricstress(leftcolumn) and Von Mises stress (right column) at four selected time steps during stage-1 and stage-2 fracturing at well W2. The time steps are the same as the figure above. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127 FIGURE 4.20 Field case study simulation: Observed poroelastic signal at wells W1 and W3 when well W2 undergoes stage-1 and stage-2 hydraulic fracturing stimulation. W2 stage-1 and W2 stage-2 curves are injection pressures monitored at two separate locations corresponding to the stage-1 and stage-2 locations; see Fig- ure 4.16. The peak in an injection pressure curve corresponds to tensile failure oftherock. Post-peakevolutioncorrespondstothefracturepropagationperiod for the respective fracture. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129 xvi ABSTRACT SYNERGISTIC COUPLING BETWEEN GEOMECHANICS, FLOW, AND TRANSPORT IN FRACTURED POROUS MEDIA: APPLICATIONS IN HYDRAULIC FRACTURING AND FLUID MIXING by Minh Tuan Tran The main focus of my research is to develop a holistic computational modeling framework that couples reservoir geomechanics, fluid flow, solvent transport, and fracture dynamics. The first objective is to quantitatively characterize the effect of geomechanical coupling on key macroscopic transport phenomena such as spreading, mixing, and viscous fingering. The second objective is to investigate how flow-transport coupling and fracture dynamics modulate the stress state and geomechanical stability of stress-sensitive subsurface formations. The third objective is to probe the impact of poroelasticity on the temporal and spatial evolution of embedded fracture networks in terms of induced permeability, dynamic propagation, and fracture interactions. The key conclusions are summarized here. I examined how viscous fingering-induced fluctuations in the flow field can increase the strength of coupling between poroelastic stresses and the mix- ing/spreading characteristics of the tracer slug. Failure to account for the stress-transport coupling will result in an inaccurate prediction of the tracer breakthrough time, tracer recovery, and the degree of mixing. I characterized the effect of geomechanical coupling on key transport metrics and quantified the impact each coupling parameter had on tracer breakthrough and the degree of mixing. I explored how the characteristics of the tracer retrieval curve were modulated by both fracture morphology and the strength of geomechanical coupling. I discovered that the transport- geomechanics synergy exists and can be modeled via numerical coupling to the pressure field. I also found that the geomechanical coupling strength governed the evolution of fracture permeability with different quantifiable impacts. I investigated how flow-transport coupling and fracture dynam- ics influence the stress state and fracture stability in the domain. Last but not least, I had shown xvii that the poroelastic pressure signal resulting from deformation could serve as a diagnostic tool for identifying fracture morphology and reservoir heterogeneity. xviii CHAPTER 1 Introduction 1.1. Motivation and Objectives Despite recent advances in modeling flow and geomechanics coupling, a holistic multiphysics ap- proach to capturing the coupling synergy between fluid flow, injectant transport, induced stresses, and fracture mechanics is lacking. There is a growing interest in developing an integrated multi- physics simulator to investigate physical processes during fluid injection into a fractured reservoir, i.e., injection-induced deformation and stresses, dissolution and spreading of the injected fluid, and convective instabilities originating from the contrast in injected and resident fluid properties, e.g. viscous fingering. The novelty of this dissertation lies in uncovering the rich interplay among these interactive dynamic physical processes through a novel computational scheme that I developed to solvethecoupledproblemofflow, transport, geomechanics, andfracturemechanicsatthefieldscale. Motivated by the difference in the mathematical nature of the governing equations, I propose a nu- merical formulation of a monolithic flow-geomechanics solver sequentially coupled with a transport solver via the velocity field and the concentration-dependent mobility model. Embedded Discrete Fracture Modeling (EDFM) is used to model the flow and transport processes in a fractured porous media while an improved Bandis model is employed to simulate the fracture-mechanical response to the changing stress state. Poroelastic coupling and EDFM implementation are validated against classic analytical and numerical problems. A case study, inspired by a huff-n-puff tracer flowback study, will be conducted to investigate the applicability of the proposed simulator in the field (Fig- ure 1). A sensitivity analysis will be performed to evaluate the dependence of global transport characteristics, permeability evolution, and fracture stability on critical parameters dictating the strength of coupling between geomechanics, flow, and transport. Finally, I employ phase field mod- eling and EDFM to simulate the growth and interaction of multiple randomly oriented fractures during fluid injection. I characterize the impact of fracture properties (i.e. spatial distribution, orientation) and formation mechanical heterogeneity on critical variables during fracture propaga- tion such as fracture length and induced poroelastic signal. A case study shows that the model is 1 capable of modeling complex underlying physical processes and displaying field-tested observations during fracture propagation with success. The findings of this study are relevant in many practical applications such as waste disposal, carbon dioxide sequestration, contaminant transport, enhanced oil recovery, and tracer surveillance. For example, integration of poroelastic coupling and fracture dynamics in transport simulations can better guide us on where to inject expensive injectant to attain ideal mixing zone, especially in reservoir embedded with stress-sensitive fractures in which the loss of circulation is a possibility. Furthermore, understandingtherelationshipbetweenfracturedynamicsandgeomechanicalcoupling and improving the mechanistic modeling of fracture permeability evolution help determine where to optimally drill infill wells in tight shales by locating regions that are most conducive to fracture opening. Likewise, the optimal conditioning of hydraulic fracture (i.e. maintaining propped opening from closure stress) is much dependent on the direction and magnitude of the principal stress state, which constantly evolves during dynamic processes such as injection and production of fluids in the development area. Besides, the relationship between fracture dynamics and poro-elasticity is critical to improving the operational design of injection/production operations by maintaining subsurface formation integrity. My research also contributes to refining the state of knowledge in tracer behavior in a fractured reservoir, assisting in the qualitative and quantitative deduction of fracture and reservoir properties from tracer survey studies. 1.2. Overview of Chapters The following chapters take turns dissecting various aspects of the main research questions stated in the abstract. The flow structure of each chapter starts with outlining the overarching challenges and introducing background literature on the topics. Then, constructions and assumptions of the phys- ical model, mathematical model, and numerical model are described in detail. Finally, numerical simulation results are presented with necessary benchmarks as well as possible applications to field data where available. Chapter 2 addresses how poroelasticity impacts solute transport processes during viscously unstable flows in a deformable aquifer. Chapter 3 investigates the role of transport- geomechanics coupling in modulating the spreading and miscibility of a solute slug during viscously 2 unstable flows. Chapter 4 quantifies the impact of poroelastic coupling and stress shadowing on dynamic fracture behaviors (i.e. nucleation, propagation, and interaction). 3 CHAPTER 2 Coupling between Transport and Geomechanics Affects Solute Spreading and Mixing during Viscous Fingering in Deformable Aquifers 2.1. Abstract We investigate the influence of poroelasticity on solute transport processes during viscously unstable flows in a deformable aquifer. I hypothesize that viscous fingering-induced fluctuations in the flow field can increase the strength of coupling between poroelastic stresses and mixing and spreading characteristics of the solute slug. To test the hypothesis, I develop a coupled flow, transport and geomechanics framework to simulate spreading, mixing and deformation mechanisms during vis- cous fingering in a stress-sensitive aquifer. I show quantitatively that concentration-dependent fluid mobility affects the pressure and effective stress fields, and deformation-driven changes in porosity affects advective and diffusive fluxes. I perform a sensitivity analysis to evaluate the dependence of global transport characteristics—the solute breakthrough time, breakthrough curve and the degree of mixing—on the strength of flow-deformation coupling defined in terms of poroelastic moduli. Results are obtained for different magnitudes of the Peclet number, solute viscosity contrast, Biot coefficient, rock bulk modulus, fluid compressibility, and the permeability field variance. I find that, in presence of less consolidated formations, softer rocks, or more compressible fluids, the injected solute breaks through earlier and mixes with the resident fluid better. Our results provide impor- tant insights for modeling hydrodynamically unstable contaminant remediation and fluid injection processes in stress-sensitive reservoirs. 2.2. Introduction Spreading and mixing of solutes in subsurface reservoirs has garnered significant interest, espe- cially in real-time surveillance applications such as tracer analysis and waste disposal monitor- ingWelty and Gelhar(1991);Kapoor and Kitanidis(1998);Le Borgne et al.(2010);Jha et al.(2013); Nicolaides et al.(2015). Whilespreadingdescribesthespatialstretchinganddeformationofasolute slugdrivenbyadvectivemechanisms, mixingisthemechanismthatincreasestheactualfluidvolume 4 occupied by the solute behind the spreading front Kapoor and Gelhar (1994); Kapoor and Kitanidis (1998);Le Borgne et al.(2010). Spreadingusuallyintensifiesthesoluteconcentrationcontrastinthe domainwhilemixingtendstosmoothouttheconcentrationgradientwithinthesolute-contaminated region Le Borgne et al. (2013); Jha et al. (2011b). These processes are also critical for enhanced oil recovery (EOR) technologies such as chemical flooding and injection of CO 2 , nitrogen, methane, and ethane gases Waggoner et al. (1992); Chen and Meiburg (1998); Lake et al. (2014). The injected fluid can be less viscous than the in-situ fluid, e.g. crude oil or brine, to improve wellbore injectivity at lower injection pressures. It can also be more viscous than the resident fluid to increase the sweep efficiency of the fluid-fluid displacement process. Also, the injected and in-situ fluids can mix and react within the reservoir, which can change the physical properties of the mixture, e.g. viscosity or density. This can impact the fluid velocity field and lead to two-way coupling between flow and transport Welty and Gelhar (1991); Waggoner et al. (1992); Tan and Homsy (1992); Jha et al. (2011a). The displacement of a more viscous fluid by a less viscous fluid leads to the Saffman- Taylor instability also known as viscous fingering Koval (1963); Homsy (1987). It is well-known that fingering reduces the recovery efficiency of fluid-fluid displacement processes such as water- flooding Lake et al. (2014). In case of groundwater contaminant treatment and transport, viscous fingering leads to a rapid increase in the macrodispersivity Welty and Gelhar (1991) and affects contaminant breakthrough and removal Nicolaides et al. (2015). The length of the mixing region is shown to increase linearly with time and monotonically with the permeability variance in one study Tan and Homsy (1992), and in another study it is shown to increase linearly with the square root of time and behave non-monotonically with the variance De Wit and Homsy (1997). In terms of the interaction between the permeability correlation length scale and the scale of fingering, a resonance mechanism was proposed to explain the maximum growth of the mixing length at inter- mediate correlation values Tan and Homsy (1992). In summary, the effect of fingering on solute mixing in heterogeneous rocks is still unclear. Often, the reservoir of interest is geomechanically sensitive and deformable due to the presence of pre-existing fractures, unconsolidated formations, excessive flow-induced stresses, or low compress- ibility fluids. The last scenario can result from an increase in the water saturation associated with 5 water injection or hydrocarbon production from an oil reservoir because water is usually less com- pressible than oil or gas Zhao and Jha (2019). The role of Biot coefficient in groundwater pumping- induced fissure formation has been hypothesized recently Budhu (2011). The effect of consolidation on solute transport has been studied in the context of loosely packed landfill columns with incom- pressible grains, where it is reported that compaction-induced decrease in porosity and increase in the particle velocity lead to a faster breakthrough Peters and Smith (2002); Alshawabkeh et al. (2005); Fox et al. (2011a); Zhang et al. (2012) compared to the case without compaction. How- ever, I do not know how the breakthrough time varies for rocks at depth, where the grains are tightly packed and cemented together, and their compressibility cannot be ignored. I hypothesize that the compaction-induced breakthrough acceleration mechanism is weaker in rocks because the magnitude of compaction is smaller in rocks than in soils. As a result of these observations, there is a growing interest in integrating rock deformation with fluid flow processes in existing reservoir models to obtain a more realistic representation of the subsurface and to improve the agreement between model predictions and field observations. Several computational frameworks of coupled flow and geomechanics have been proposed for ap- plications in hydrocarbon recovery, groundwater extraction, underground waste disposal, and CO 2 storage Lewis and Sukirman (1993a); Gawin et al. (1995); Rutqvist et al. (2013); White and Borja (2007); Ferronato et al. (2010); Jha and Juanes (2014); Castelletto et al. (2015a). Effect of stress on permeability under the assumption of uncoupled pressure and stress fields has been studied recently using a pore scale model Roded et al. (2018). However, most of these studies neglect the two-way coupling mechanisms between transport and deformation i.e. deformation-induced changes in the solute transport behavior and concentration-dependent fluid mobility effects on deformation. The effect of geomechanical coupling on various transport mechanisms such as spreading, mixing, and viscous fingering in stress-sensitive reservoirs, is unexplored Tran et al. (2018). A study on the effect of stress on contaminant transport in fractured crystalline rock showed that the residence time of contaminant is a strong function of the boundary stress ratio Wang et al. (2017). Wefocusonunravelingthetransport-deformationcouplingmechanismsatreservoirconditionswhen 6 the two-way coupling between flow, transport and deformation cannot be neglected. In particular, I focus on the mechanisms that affect mixing and spreading because these processes control the ultimate recovery during contaminant removal or enhanced oil recovery. We present a numerical framework to model and analyze the effect of poroelasticity on tracer transport under viscous finger- inginstability. Iusetheframeworktoevaluatethesensitivityoftransportmetricstoporomechanical properties including the Biot coefficient and the bulk modulus. 2.3. Physical Model We consider a 2D section of an aquifer bounded by a horizontal injection well and a horizontal production well (Figure 2.1). The reservoir is fully saturated and under hydromechanical equilib- rium with external stresses σ h and σ H , which are the minimum and maximum principal horizontal compressions from regional tectonics. The injection well injects two fluids of different viscosities alternatingly in the form of slugs. The fluids are perfectly miscible with each other and with the in-situ fluid. One of the fluids has the same viscosity as the in-situ fluid, and the other fluid is less viscous due to the presence of a chemical tracer. The production well is open to flow at a fixed bottom hole pressure lower than the reservoir pressure. As the less viscous slug moves through the reservoir and displaces the more viscous in-situ fluid, viscous fingering ensues at the leading front of the less viscous slug. This affects the spreading and mixing of the slug as well as the recovery of the in-situ fluid from the producer. Since the reservoir pressure is coupled to both slug transport and the state of stress, I can use this model to probe the coupling between transport and deformation processes and extract possible physical mechanisms arising from such coupling. We make the follow- ing assumptions, which are appropriate for the physical problem at hand: (1) slightly compressible fluid, (2) small strain, (3) quasi-static equilibrium, (4) linear elasticity, (5) single-phase flow, (6) isotropic material behavior, (7) permeability not dependent on deformation, (8) isothermal system (no temperature change), and (9) first contact miscibility between the fluids. 2.4. Mathematical Model The starting point is to state the governing equations of flow, deformation and transport in the strong form, which will be converted into the respective weak forms and discretized later. I adopt 7 Injector Producer Figure 2.1: A physical model of alternating slug injection into a reservoir under a tectonic state of stress created by compressions σ h and σ H on its two sides. Light-colored (viscosity µ 1 ) and dark-colored (viscosity µ 0 >µ 1 ) slugs are injected alternatingly into the reservoir from an injection well, and the fluids are produced from the producer. Displacement of the more viscous fluid by the less viscous fluid leads to viscous fingering. Fingering-induced perturbation in the pressure field generatesperturbationintheeffectivestressfield,whichaffectsdeformationandhydraulicproperties of the reservoir, which in turn can affect spreading of the less viscous fluid and production of the more viscous fluid. 8 a continuum representation of a macroscopic model of porous media, where fluid and solid are viewed as overlapping continua Bai and Elsworth (2000b); Coussy (2004). The model domain has an area Ω = ( L x × L y ), where L x and L y are the dimensions in the x and y directions. I denote the closed boundary of the domain by Γ . The initial porosity is ϕ 0 , and the permeability k is assumed to be isotropic and constant in time. The domain is saturated with the fluid of density ρ f,0 . Subscript 0 indicates initial condition. The fluid mass conservation equation for single-phase slightly compressible flow and infinitesimal deformation with no sources or sinks in the domain is 1 M ∂p ∂t +α ∂ϵ v ∂t +∇· v =0, (2.1) where p is pore pressure, t is time, 1/M = ϕ 0 c f + (α − ϕ 0 )/K s is the inverse Biot modulus, α = 1− (K/K s ) is the Biot coefficient indicating the degree of cementation between grains, K is the drained bulk modulus of rock, K s is the bulk modulus of solid grains, c f =(1/ρ f )dρ f /dp is the fluid compressibility, ρ f is the fluid density, v is the Darcy fluid velocity vector relative to the solid skeleton, and ϵ v is the volumetric strain. The volumetric strain rate term in the equation above couples fluid flow to mechanical deformation and acts as a local source term. Denoting the solid skeleton displacement vector with u and the skeleton velocity with v s = ∂u/∂t, and noting that ϵ v =∇· u, I obtain, 1 M ∂p ∂t +∇· (v+α v s )=0. (2.2) The fluid momentum balance equation or Darcy 's law is: v = k µ (c) (ρ f g−∇ p), (2.3) where µ (c) is the concentration-dependent fluid viscosity, and g is the gravitational vector. I use a Gaussian random field model with the following parameters to model permeability heterogeneity: geometric mean permeability k g , dimensionless log-permeability variance σ 2 logk , and permeability correlation lengths ξ x and ξ y . An example is shown in Fig. 2.2. 9 Figure 2.2: A Gaussian random permeability field. Values of the Gaussian model parameters are as follows: geometric mean k g =50 md, log-k variance σ 2 logk =0.01, and correlation lengths ξ x =ξ y =1 m. The color scale is in millidarcy unit. The quasi-static linear momentum balance equations for a fluid-saturated rock is ∇· σ +ρ g =0, (2.4) where σ = σ ′ − αp I d is the Cauchy total stress tensor, σ ′ is the effective stress tensor, I d is the Identity tensor, ρ = ϕρ f + (1− ϕ )ρ s is the bulk density, and ρ s is the solid grain density. Here, ϕ is the Eulerian or ‘true’ porosity in the current configuration. A sign convention where the normal stresses are positive in tension is used. The total volumetric stress σ v = trace(σ )/n d and the effective volumetric stress σ ′ v are related by σ ′ v =σ v +αp where n d is the number of space dimensions (n d =2 in 2D).σ ′ v andp are also related by the compatibility relation∇ 2 σ ′ v =α ∇ 2 p. The isotropic part of the deformation ϵ v (volumetric expansion or contraction) is controlled by σ ′ v and the deviatoric part (rock failure and plasticity) is controlled by the von Mises stress σ vm = √ 3J 2 , where J 2 is the second invariant of the deviatoric stress (σ − σ v I d ). Substituting the effective stress into Equation (2.4), I obtain ∇· σ ′ − α ∇p+ρ g =0. (2.5) Assuming a linear elastic behavior of the rock,σ ′ =Dϵ , whereϵ =(∇u+∇u T )/2 is the infinites- imal strain tensor defined as the symmetric gradient of the displacement vector. The 2D domain represents a thin layer of rock, thus a state of plane stress is assumed. In the compact engineering 10 notation, σ = [σ xx ,σ yy ,σ xy ] T , ϵ = [ϵ xx ,ϵ yy ,2ϵ xy ] T , and I d = [1,1,0] T are the stress, strain and identity vectors. In the plane stress configuration, the isotropic drained elasticity tensor D can be written as D = E 1− ν 2       1 ν 0 ν 1 0 0 0 (1− ν )/2       (2.6) where drained Young’s Modulus E and Poisson’s Ratio ν are the two parameters of the isotropic elastic material model. In 2D, σ ′ v = Kϵ v , with K = E/(3(1− 2ν )). The fluid mass conservation and the mechanical equilibrium equations are coupled through Biot’s poroelasticity equation, which relates the fluid mass increment per unit initial bulk volume, δm = δ (ρ f ϕ (1+ϵ v )), to strain and pressure as δm/ρ f,0 = αδϵ v +δp/M . We can obtain the equation for porosity evolution from the solid mass conservation equation Lewis and Schrefler (1998) and partitioning of the total stress on porousskeletonbetweenthestressactingonthesolidgrainsandthepressureoftheporefluidCoussy (2004); Jha and Juanes (2014): dϕ dt = ∂ϕ ∂t +∇ϕ · v s =(α − ϕ ) ∂ϵ v ∂t + α − ϕ K s ∂p ∂t (2.7) which can be integrated in time and approximated under infinitesimal strain as δϕ =(α − ϕ )(δϵ v +(1− α )δp/K ). (2.8) There are several interesting limiting cases of Equation (2.1). The first case is that of almost incompressible solid and fluid, i.e. c f → 0, K s → ∞, M → ∞. To maintain finite values of the pressure, this results into an instantaneous and local equilibration of pressure such that δm = ρ f,0 αδϵ v or δϕ =αδϵ v . This implies ∇· (v+α v s )=0 ∇· v =− α (∂ϵ v /∂t)=− ∂ϕ/∂t. (2.9) 11 This is relevant to model flow in shallow aquifers or underneath dams, away from wells, where flow is driven by the volumetric expansion or contraction of the pores. The second case is that of flow through low porosity consolidated rocks i.e. ϕ ≪ 1, K → K s , α → 0, which results into the traditional flow equation uncoupled to mechanical deformation: (ϕ 0 c f +c r )∂p/∂t+∇· v =0, (2.10) where c r is the rock compressibility Aziz and Settari (1979). This approximation is also valid in highly stiff rock ( K → K s ) saturated with highly compressible fluid (large c f ) such as gas. This case is relevant for modeling flow in conventional, strongly-consolidated petroleum reservoirs. The third case is undrained deformation, low permeability i.e. δm →0, k→0, i.e. v→0 such that (1/M)∂p/∂t+α (∂ϵ v /∂t)=0 δp u =− αMδϵ v =− Bδσ v (2.11) where B = αM/K u is the Skempton’s coefficient, K u = K +α 2 M is the undrained bulk modulus, andδp u is the undrained pressure change. This case is relevant for monitoring stress changes due to production, andpressurechangesduringhydraulicfracturingDaneshy(2015);Kampfer and Dawson (2016a); Roussel and Agrawal (2017a) in unconventional rocks such as shales. The fourth case is drained deformation i.e. δp → 0, which is the same as the first case in terms of the resulting equations. The Advection Dispersion Equation (ADE) for tracer transport is ∂(ϕc ) ∂t +∇· (vc− ϕD d ∇c)=0 (2.12) Here,c is the normalized local concentration of tracer, andD d is the hydrodynamic dispersion coeffi- cient assumed to be a constant. I assume that velocity-dependent dispersion is less important for the rectilinear injection case considered here compared to the radial injection case Chen and Meiburg (1998); Nicolaides et al. (2015) where the average velocities decrease rapidly with the distance from 12 injection source. Re-arranging the terms to extract the effect of geomechanical coupling on trans- port, ϕ ∂c ∂t +∇· (vc)− ϕD d ∇ 2 c=D d ∇ϕ ·∇c− c ∂ϕ ∂t (2.13) The left-hand side is the residual of the traditional ADE uncoupled to geomechanics, and the right- hand side is the mechanics-to-transport coupling term composed of source/sink due to variations in porosity, which can be calculated using Equation (2.8). We use an exponential model for the concentration-dependent viscosity of the mixture, µ (c)=µ 1 e R(1− c) , (2.14) where R=log e µ 0 /µ 1 is the log viscosity ratio of the in-situ fluid viscosity µ 0 and the injected fluid viscosity µ 1 . When 0<c<1, depending on the value of R, the viscosity of the mixture is higher than the viscosity of injected tracer. This promotes viscous fingering at the front between the tracer and in-situ fluids, when the former displaces the latter. The mathematical model requires boundary and initial conditions for well-posedness. For the rect- angular domain considered in Figure 2.1, the normal displacement is zero on the bottom and right boundaries, compressions are applied on the left and top boundaries, and normal fluid flux is zero on the left and right boundaries. The producer is controlled by a prescribed bottom hole pressure p out , and the injector is controlled by an injection rate q in (per unit length in the third dimension) prescribed at the bottom boundary. The flow rate at the producing boundary q out can be calculated from p out , and the flow area (length) L x . Natural outflux boundary condition is imposed on the top, left, and right boundary. Initial pressure, concentration, stress, displacements and porosity are prescribed to achieve equilibrium with the boundary and body forces at t = 0. Governing param- eters of the problem are R, Pe, α , K, c f , σ 2 logk , and ξ y /ξ x . Boundary compressions, pressure, flow rate, k g and slug length prescribed as part of the boundary conditions can also be included in the list of control parameters. 13 2.4.1. Dimensionless forms It is instructive to look at the dimensionless forms of the governing equations. Choosing the do- main length L x as the characteristic length, p c as the characteristic pressure (or stress), µ 1 as the characteristic viscosity, k g as the characteristic permeability, v c = q in /(h x h z ) as the characteristic Darcy velocity (assuming h z = 1 meter), ϕ 0 as the characteristic porosity, t c = ϕ 0 L x /v c as the characteristic time, and neglecting gravity, I obtain the dimensionless flow equation: p c M ∂p D ∂t D +α ∂ϵ v ∂t D − ϕ 0 ∇ D · k D e R(c− 1) ∇ D p D =0, (2.15) where the subscript D indicates dimensionless quantities. k D (x D ,y D ) is a multi-Gaussian function of σ 2 logk and ξ y /ξ x . The dimensionless form of the mechanics equations becomes ∇ D · σ ′ D − α ∇ D p D =0 (2.16) The dimensionless form of the porosity evolution equation becomes δϕ D =((α/ϕ 0 )− ϕ D )(δϵ v +(1− α )(p c /K)δp D ). (2.17) The dimensionless form of the ADE is ∂(ϕ D c) ∂t D +∇ D · (v D c)− 1 Pe ∇ D · (ϕ D ∇ D c)=0, (2.18) where v D = − k D e R(c− 1) ∇ D p D and Pe = v c L x /(ϕ 0 D d ) is the Peclet number. I can choose p c as equal toK or define it from the prescribed injection rate as p c =q in µ 1 /k g . Here, I choose the latter. Below, I drop the subscript D when referring to the dimensionless quantities. 2.4.2. Spreading and mixing We will characterize spreading of the tracer slug by the tracer breakthrough time t br and the break- through curve c prd (t)=c(L y ,t) at the producer, where the transversely-averaged concentration c is 14 defined as: c=(1/L x ) Z Lx 0 c(x,y)dx (2.19) Since the general flow is along the y-direction or the longitudinal direction, evolution of transversely- averaged quantities from the flow, transport and geomechanics sub-problems can help us under- stand the macroscopic effect of viscous fingering and poroelasticity on transport. In my previous work Jha et al. (2011a); Jha (2013), I have analyzed the transversely averaged concentration in detail. We can see the relation between transversely averaged pressure p, which controls stress changes, the Darcy velocity v y , which controls the breakthrough time, and the Biot parameters α and M, by applying the transverse averaging operator to Equation (2.2), 1 M ∂p ∂t + ∂ ∂y (v y +α v s,y )=0 (2.20) q in and p out , which are driving longitudinal flow through the domain, are constants in time. We know that as the permeability variance σ 2 logk or the viscosity ratio R increases, the tracer breaks through at the producer earlier i.e. t br becomes smaller Waggoner et al. (1992); Chen and Meiburg (1998); Berkowitz et al. (2006); Le Borgne and Gouze (2008); Nicolaides et al. (2015). However, the effects of σ 2 logk and R on the breakthrough curve are hypothesized to be different: higher σ 2 logk should lead to enhanced mixing at the tip of tracer channels and, therefore, lower breakthrough concentrations, whereas higher positive values of R should lead to reduced mixing at the finger tip and higher breakthrough concentrations. This reasoning is aligned with the observation that local- scale mixing influences global spreading measures as reported previously Kapoor and Gelhar (1994); Le Borgne et al.(2010). Themixingprocessinporousmediaflowscanbedescribedbytheevolution ofthedegreeofmixingandtherateofmixingKapoor and Kitanidis(1998);Le Borgne et al.(2010); Jha et al. (2011a); Jha (2013); Nicolaides et al. (2015), which is based on the idea of mixing in turbulent flows Jones and Launder (1972); Eswaran and Pope (1988); Pope (2000). The degree of mixing χ , which quantifies the amount of mixing in the domain, can be defined in terms of the concentration variance σ 2 c or the scalar energy, and the rate of mixing can be quantified in terms of the mean scalar dissipation rate ϵ c Jones and Launder (1972); Eswaran and Pope (1988); 15 Le Borgne et al. (2010); Jha et al. (2011a,0,0): χ (t)=1− σ 2 c (t) σ 2 c,max , σ 2 c =⟨(c−⟨ c⟩) 2 ⟩, ϵ c =⟨|∇c| 2 ⟩/Pe, (2.21) where⟨·⟩ denotes spatial averaging over the entire domain. I fix ⟨c⟩=0.5 and σ 2 c,max =0.25 in the definition above to understand lack of mixing as deviation from the perfectly mixed state of 50% tracer concentration, which will correspond to σ 2 c =0 and χ =1. In a completely segregated state, e.g. c=0 or c=1, σ 2 c =σ 2 c,max and χ =0. I denote the degree of mixing at the breakthrough time as χ br = χ (t br ). ϵ c quantifies the amount of concentration gradient or, equivalently, the fluid-fluid interfacial area in the flow. A higher value of ϵ c means a longer and sharper interface, which leads to a larger value of the total diffusive flux across the interface, which leads to a larger value of dχ/dt . Therefore, ϵ c can be understood as a measure of the mixing rate. ϵ c can be measured during a lab or field experiment by measuring the average viscosity of the fluids in the domain at discrete time steps and calculating the time rate of change of the average viscosity Jha (2013). The degree of mixing evolves as follows: dχ/dt = 2ϵ c σ 2 c,max +source due to injection− sink due to production +source/sink due to porosity coupling. (2.22) Similarly, the evolution equation for ϵ c should have two source terms due to fingering-induced interface stretching and entrance of new interfaces from the injection boundary, and it should have three sink terms due to merging of fingers, diffusion, and the production boundary. Stretching and merging are nonlinear interactions that arise due to coupling between the velocity field and the concentration gradient field. For analysis of the simulation results below, I also define σ 2 c =c 2 − c 2 , χ = 1− (σ 2 c /σ 2 max ), and ϵ c =|∇c| 2 /Pe similar to the definition of the breakthrough concentration in Equation (2.19). 16 Fickian mixing For Fickian or purely diffusive mixing of linear slugs or blobs transported through a homogeneous porous medium in absence of geomechanical coupling, i.e. R = 0,σ 2 logk = 0, α = ϕ ≪ 1, the evolution of ϵ c has been studied. The 1D analytical solution for ϵ c in such diffusive mixing suggests: ϵ 1D c (t) = C 2 0 L x t − 3/2 /(8 √ 2π Pe), i.e. a power-law slope of − 3/2 in the logϵ c vs. logt plot with C 0 as a constant Le Borgne et al. (2010). However, the evolution of ϵ c in 2D periodic flow transitions from early time power-law to intermediate time exponential to late time error-function behavior Jha (2013). This corresponds to different regimes of slow and fast decay of the scalar energy σ 2 c , and it is captured by the following analytical model Jha (2013): ϵ c (t)=e − C 2 t /(C 1 + √ tPe), (2.23) where C 1 and C 2 are positive constants. This is a solution of the following evolution equation: dϵ c dt +C 2 ϵ c + 1 2 √ Pe ϵ 2 c √ te − C 2 t =0, (2.24) 2.5. Numerical Model Numerical simulation of hydrodynamically unstable processes such as viscous fingering is sensitive to discretization or truncation errors, which can grow in time due to the nonlinearity of the equa- tions. Because of different smoothness and continuity requirements of the pressure, displacement, and concentration functions, the three PDEs must be treated differently during discretization. For example, the miscible flow pressure field is expected to be smooth at all times whereas the concen- tration field is expected to show sharp gradients due to fingering that evolves in time. Higher-order discretization schemes are required in ADE to resolve the sharp concentration gradients Jha et al. (2011a,0). Also, the fluid mass conservation must be honored locally for each element, which can become a challenge in heterogeneous flow fields. Given these and other considerations around simu- lation runtime, the finite volume method has proven to be a robust method for solving the flow and transport equations. The mechanics equation is suited for a nodal finite element method because it 17 can satisfy force balance at each node in the mesh. In absence of faults/fractures, the displacement field is expected to be free of any discontinuities. However, a smooth stress field requires the dis- placement to be discretized with at least piecewise linear functions. In mixed field problems such as the one I have, the choice of discrete function spaces for pressure and displacement must satisfy the stability and solvability conditions, e.g. the Ladyzhenskaya-Babuska-Brezzi (LBB) condition Brezzi (1974); White and Borja (2007), to simultaneously solve the equations when both the fluid and solid become nearly incompressible, e.g. immediately after wells start to produce or inject. To satisfy or avoid the LBB condition at all times, I consider compressible systems with finite positive values of K s and c f as well as reasonably large time step sizes, which are representative of many petroleum reservoirs and deep aquifers. Motivated by these considerations, I use a finite volume–finite element formulation to solve my cou- pled problem Jha and Juanes (2014). I derive the weak forms of the flow and mechanics equations as follows: multiplication of the strong form by a suitably smooth test function, integration over the domain, integration by parts to reduce the order of differentiation on the unknown function, appli- cation of the divergence theorem, and imposition of the boundary conditions. The flow problem is discretized in space with cell-centered piecewise constant pressures and two-point flux approxima- tion (TPFA). The mechanics problem is discretized with bilinear elements and nodal displacement vectors Jha and Juanes (2014). For compressible systems considered here, this choice is not subject to the LBB or inf-sup condition. The transport problem is discretized using a sixth-order accurate compact finite difference method Lele (1992); Jha et al. (2011a). The flow problem is integrated in time using the first-order accurate implicit Backward Euler scheme. The transport equation is integrated in time using the third-order accurate explicit Runge-Kutta scheme. We solve the di- mensionless forms of the equations over a dimensionless time period 0<t<3, which allows complete drainage of the injected slug. 2.5.1. Coupling strategy and solution scheme We utilize a fully-coupled simultaneous solution approach to solve for pressure and displacement fieldsateachtimestepLewis and Sukirman(1993b);Pao et al.(2001);Lewis et al.(2003);Jha and Juanes 18 (2007); Ferronato et al. (2010). The system of linear equations becomes:    K − Q Q T S+∆ tT       U p    n+1 =    f u ∆ tf p +Q T U n +Sp n    (2.25) whereU is the vector of nodal degrees of freedom from the finite element interpolation u=N T U, N is the matrix of bilinear shape functions,p is the vector of cell-centered pressure unknowns, n is the previous time step, and n+1 is the current time step at which displacements and pressures are sought such that ∆ t=t n+1 − t n . The global stiffness matrix is calculated by assembling individual element stiffness matrices, K = Z Ω i B T i D i B i dΩ (2.26) whereB is the strain-displacement matrix computed by applying the symmetric gradient operator toN, and i is the element or cell index. For 2D quadrilateral elements with eight nodal degrees of freedom per element,N is a 2× 8 matrix,B is a 3× 8 matrix, andK i is a 8× 8 matrix. The global transmissibility matrix T is calculated by assembling the individual interfaces transmissibilities. For the interface ij between cell i and cell j, the interface transmissibility is T ij = Γ ij /[(l i µ i /k i )+ (l j µ j /k j )], which is the harmonic mean of fluid mobilities k i /µ i and k j /µ j weighted by the cell centroid-to-face centroid distances l i and l j within the two cells. Γ ij is the interface area between two cells. The global storativity or compressibility matrix is a diagonal matrix with components, S ii = R Ω i 1/M i dΩ . The flow-geomechanics coupling matrix at the element level is, Q i = Z Ω i α B T i I d dΩ . (2.27) The right hand side of the weak form of the mechanics problem consists of traction and body forces, f u = Z Γ t N T σ dΓ+ Z Ω N T ρ g dΩ (2.28) where Γ t denotes Neumann boundaries of the domain with the prescribed tractionsσ =[0,σ h ] T on the top boundary and σ = [σ H ,0] T on the left boundary. The right hand side of the weak form of 19 the flow problem consists of prescribed fluxes and well terms. For cell i, f p i =q in − q out . (2.29) We solve Equation (2.25) using a direct solver. Next, I discuss the solution of the ADE Eq. (2.18). Coupling between the poroelastic solver and the transport solver is implemented using a sequential one-way coupled approach. For relatively small values of R as considered here, the transport-to-deformation coupling strength is expected to be weak because there is no explicit appearance of a concentration functional in Eq. (2.16), and concentration-dependence of pressure is limited to fluid mobility. This suggests that a one-way coupled approach to solve the two-way coupled transport-deformation problem can be implemented to reduce computational cost without sacrificing accuracy. I will test this hypothesis in Section 5 by analyzing the magnitude of fingering-induced stress and porosity changes resulting from transport- to-deformation coupling. We discretize Eq. (2.18) in space using a sixth-order accurate compact finite difference method Lele (1992); Jha (2013). I compute the discrete interface velocities using a two point flux approximation of Eq. (2.3), i.e. v ij =T ij (p i − p j )/Γ ij , where v ij is the velocity across the ij interface between cells i and j, and p i and p j are the pressures of the two neighboring cells. We compute the N x × N y matrices of cell-centered face-normal velocitiesv x andv y using a linear interpolation of the interface velocities. We integrate the semi-discrete ODE system over the time step ∆ t of Eq. (2.25) using a third-order explicit Runge-Kutta (RK) method. The residual function of the RK integrator, which is the discretized form of the flux divergence term in Equation (2.18), is: f c =− v x .∗ (D x c)− v y .∗ (D y c)− c.∗ (D x v x +D y v y )+ 1 Pe ϕ .∗ (D xx +D yy )c + 1 Pe [(D x ϕ ).∗ (D x c)+(D y ϕ ).∗ (D y c)]− 1 ∆ t c.∗ δϕ where D x , D y , D xx , and D yy are the discrete N x × N y differentiation matrices for the ∂/∂x, ∂/∂y, ∂ 2 /∂x 2 , and ∂ 2 /∂y 2 operators, respectively. The differentiation matrices are pre-computed 20 Parameter Symbol Value Number of cells in x direction N x 100 Number of cells in y direction N y 200 Domain length in x direction L x 50 meter Domain length in y direction L y 100 meter Boundary compressions σ h , σ H 10 MPa Initial pressure p 0.1 MPa Producer bottom hole pressure p out 0.5 MPa Injector bottom hole flow rate q in 1.8× 10 − 7 m 3 /sec/m Table 2.1: Simulation parameters outside the time loop using standard methods Lele (1992). c, ϕ , and δϕ are N x × N y matrices of cell-centered concentration, porosity and porosity increment fields evaluated at t n . v x and v y at t n+1 are available from the poroelastic solver. δϕ n+1 is calculated using Eq. 2.17, and I update porosity as ϕ n+1 = ϕ n +δϕ n+1 . The .∗ operator above denotes the element-wise multiplication of two arrays. An adaptive time stepping method based on the Courant–Friedrichs–Lewy (CFL) condition is used to ensure stability of the explicit transport solver. The dimensionless time step is calculated as ∆ t = γL x /max(|v|) where max(|v|) is the maximum magnitude of the dimensionless velocity over the entire domain. The value of γ is chosen to ensure stability and smoothness of the coupled solution. 2.6. Numerical Simulation Here, I present simulation parameter values and design decisions that I made to define my Base Case model. I discretize the rectangular domain in Figure 1 using a mesh with square elements. The mesh has a sufficient number of elements to ensure that the fingers are resolved, and the simulation finishes within a reasonable time frame (2-6 hours) on my cluster. See Table 2.1 for values of the mesh, initial condition, and boundary condition parameters. Values of other parameters in the Base Case model are given in Table 2.2. 2.6.1. Initialization Initialization of coupled flow-geomechanical models is a non-trivial task Settari and Walters (2001). Ideally, the initial state (initial displacements and pressure) should be in mechanical and hydrostatic equilibria such that the initial strains and fluid fluxes inside the domain are zero. Theoretically, 21 Parameter Symbol Value Young’s modulus E 5 GPa Poisson’s ratio ν 0.25 Initial grain density ρ s,0 2650 kg/m 3 Drained bulk modulus K 3.33 GPa Biot modulus M 6.25 GPa Diffusion coefficient D d 4.6× 10 − 7 m 2 /sec Peclet number Pe 400 In-situ fluid viscosity µ 0 1.2× 10 − 3 Pa-sec Initial fluid density ρ f,0 1000 kg/m 3 Fluid compressibility c f 1× 10 − 9 Pa − 1 Permeability k g 50 millidarcy Initial porosity ϕ 0 0.1 Biot coefficient α 0.5 Log viscosity ratio R 2 Log permeability variance σ 2 logk 0.01 Table 2.2: Base Case parameters this is achievable by prescribing a total initial stress tensor for each grid cell, which balances the prescribed boundary stresses, initial pressure, and body forces in each direction, and by prescribing an initial pressure distribution that honors the hydrostatic potential gradient (including the effects of depth-dependent compressibility and fluid-contacts in case of multiple fluids). In complex cases, e.g. heterogeneous distribution of ρ f,0 , α , plastic deformation, or time-dependent boundary condi- tions, achieving mechanical equilibrium at t=0 can become a challenge and may require a separate initialization step during the simulation. This is the approach I follow here to create a more robust formulation that is suitable for real-world applications. The stresses are initialized from the results of the first timestep, during which the domain deforms under the action of boundary loads, initial pressure, and body forces. The displacement field at the end of the initialization timestep is des- ignated as the reference displacement field. Deformation at subsequent time steps are calculated from that reference configuration. Ourinitialconcentrationfieldconsistsofapartiallyinjectedtracersluginthedomainwithadiffused interface between the two fluids. Initiation of fingering at the unstable interface requires a small perturbation in the velocity field, which is achieved by perturbing either the permeability field or the concentration-dependent viscosity field. I choose the latter alternative by superposing small ampli- 22 tude, random wave number perturbations on the c=0.5 contour at the interface Jha et al. (2011b). BasedontheparametersoftheproblemandthephysicsofSaffman-TaylorinstabilityHomsy(1987), a few of these wave numbers or modes are selected as the unstable modes, which grow to become fingers, while other modes are damped and suppressed. This is a result of the nonlinear dispersion relation between the growth rate and wavenumber Tan and Homsy (1986), which predicts short- wavelength stability (cutoff wavenumber) and long-wavelength instability Yortsos and Hickernell (1989); Loggia et al. (1995) and the presence of a selection mechanisms for the finger size in the transverse direction. I utilize this mechanism in initiating the instability from a broadband random perturbation of the initial front, which does not affect the fingering dynamics in my simulations at intermediate and late times. 2.6.2. Alternating injection We define an alternating injection cycle as the injection of one slug of the tracer fluid for a dimen- sionless duration of t inj followed by one slug of the more viscous fluid for the dimensionless duration t gap . The total duration of one cycle is t cyc = t inj +t gap . I implement this boundary condition at the injector as follows: c| y=0 =        1, if(i cyc − 1)t cyc <t≤ (i cyc − 1)t cyc +t inj , 0, ifi cyc t cyc − t gap <t≤ i cyc t cyc , (2.30) where 1≤ i cyc ≤ n cyc is the alternating injection cycle number, which is an integer, and n cyc is the total number of cycles. I choose t inj =0.7 and t gap =1.4. In this study, I fix n cyc =1. The simulation stops when the entire tracer slug has exited the domain. The breakthrough time t br is defined when the maximum tracer concentration in the producer exceeds 0.1, which occurs when the fingers of the first tracer slug arrive at the producer. The drainage time t dr is defined when no cell in the domain has concentration exceeding 0.01, meaning most of the tracer has exited the domain. 2.6.3. Sensitivity simulations To better understand the coupling mechanisms between transport and deformation processes in real reservoirs, whose properties can vary over a wide range, I perform a sensitivity analysis to 23 investigate the influence of rock-fluid properties on the breakthrough and mixing results of the Base Case defined above. I choose the following parameters of the problem as the sensitivity parameters: the Biot coefficient α , Peclet number Pe, log viscosity ratio R, permeability field parameters ( σ 2 logk andξ y /ξ x ), fluidcompressibility c f , andtherockdrainedbulkmodulusK. Eachsimulationprovides spatial maps of nodal displacements U x and U y , cell-centered pressure p, effective volumetric stress σ ′ v , Von Mises stress σ vm , concentration c, porosity ϕ , and face-centered velocities v x and v y at user-selected timesteps. 2.7. Results Below I discuss the Base Case results, extract the transport-deformation coupling mechanisms, and identify mechanisms’ sensitivity to important poroelastic parameters. 2.7.1. Base Case Figure 2.3 shows the development of viscous fingers during the simulation through nonlinear interac- tionsoftip-splitting,shielding,andcoalescenceoffingersTan and Homsy(1988);Zimmerman and Homsy (1992). Transversely-averaged profiles of concentration and change in pressure are also shown at threetimesteps, includingatthebreakthroughtimet br =1.16. Withtime, boththeleadingandthe trailing fronts of c stretch along the flow direction. The leading front stretches because of fingering, and the trailing front stretches due to diffusion. The pressure profiles can be explained as follows. As the volume occupied by the less viscous slug increases due to spreading, the pressure decreases in the domain. Similarly, when t>t br , the domain is increasingly occupied by the more viscous fluid and pressure increases because I maintain the same injection rate. This is the first-order effect of tracer transport, which does not show the detailed effect of fingering on deformation. For the cho- sen values of the governing parameters, the pressure profiles are smooth because of the dominantly diffusive nature of the pressure equation as opposed to the dominantly hyperbolic nature of the ADE. This also holds true for profiles or 2D maps of the effective volumetric stress σ ′ v (not shown) because of the stress compatibility relation between laplacians of σ ′ v and p mentioned earlier. In general, the pressure changes are negative, which lead to negative or compressive change in σ ′ v and ϕ leading to reservoir compaction as also shown in Figure 2.5. 24 -3.5 -2 0 0 1 0.5 (MPa) 0 0.5 1 time a b Figure 2.3: Base Case: (a) Concentration evolution is shown through snapshots in time. I see the development of fingers and classical viscous fingering mechanisms—tip-splitting, shielding, and merging of fingers. (b) Transversely-averaged concentration c and change in pressure ∆ p (from t 0 ) are plotted along the flow direction at three time steps. Fingers spread the tracer downstream of the main slug as shown by c curves at t 1 and t br . ∆ p goes up and down with changing mobility of the fluid in the system; the yellow fluid is more mobile than the blue fluid. ∆ p at the producer is always zero because of its fixed pressure boundary condition. 2.7.2. Transport–deformation coupling mechanisms We are interested in observing the effects of poroelasticity on the tracer transport behavior and thereby extracting the transport-deformation coupling mechanisms. Figure 2.4 shows important quantities from the three processes–flow, deformation and transport–that are expected to define such coupling mechanisms. To that effect, I analyze the gradient fields of pressure, stress and porosity. We choose ∇p because it drives advection in ADE and appears as a body force in the mechanical equilibrium Equation 2.5. I choose∇σ ′ v because it drives volumetric deformation as per Hooke’s law. σ vm can also be analyzed to evaluate the shear deformation tendency. I did not include a separate plot for σ vm in Figure 2.4 because, for the Base Case parameter values, σ vm behaves almost identical to σ ′ v . This happens because, initially, the reservoir is in a principal state of stress with zero shear, and hydrostatic loading from tracer transport does not induce significant shear. Finally, I choose∇ϕ and dϕ/dt because they appear in the deformation-to-transport coupling term on the RHS of ADE Equation 2.13. Next, I investigate the results at the observation point, which is at the center of the domain, (L x /2,L y /2). As the tracer slug passes through the observation point, the concentration first 25 0 0.2 0.4 0.6 0.8 1 -5 -3 -2 -1 -4 5 4 3 2 1 -2 -1 -3 -0.1 -0.2 -0.3 -0.4 -0.5 c 0 x10 -3 d φ/dt Figure 2.4: Base Case results at the time of tracer breakthrough at the producer. The concentration field c shows viscous fingers. Magnitudes of the pressure gradient |∇p|, the effective volumetric stress gradient |∇σ ′ v |, and the porosity gradient |∇ϕ | fields are shown. The white curves are the transversely averaged profiles of the 2D fields, scaled to fit within the domain width of [0,L x ]. The similarity between the average concentration profile and the average gradient profiles confirms that viscous fingering affects, and is affected by, the gradient of pressure, stress and porosity fields. The time rate of change in porosity, dϕ/dt , also shows a signature of the fingered region. -1 0 1 [-] -0.4 -0.2 0 -2 -1 0 1 MPa meter 0.5 1 1.5 2 2.5 -0.01 0 0.01 % Figure 2.5: Base Case results at the observation point at the center of the domain. Change in concentration ∆ c causes, and is coupled to, the changes in pressure ∆ p, effective volumetric stress ∆ σ ′ v , von Mises stress ∆ σ vm , porosity ∆ ϕ , and deformations ∆ U x and ∆ U y . t is the dimensionless time. Top: features of the rate of change of porosity, dϕ/dt , are related to the concentration evolution, maximum magnitudes of U x and U y , and tracer breakthrough which occurs at t br . The three vertical lines are time markers that relate features in the dϕ/dt curve to breakthrough and maximum deformation events. 26 increases and then decreases in time as shown by ∆ c in Figure 2.5. A drop in pressure at the ob- servation point leads to compressive stress change ∆ σ ′ v <0, compaction ∆ U y <0, and loss of porosity ∆ ϕ< 0. ∆ U x becomes positive due to the Poisson effect. Due to equal values of the boundary stresses and the hydrostatic loading from slug injection, σ ′ xx ≈ σ ′ yy , and the shear stress change is negligible ∆ σ xy ≈ 0. This leads to almost identical changes in ∆ σ ′ v magnitude and ∆ σ vm . Figure 2.6 shows the evolution of mixing and spreading metrics: χ , ϵ c , and c prd . As the injection progresses, χ increases due to fingering-induced spreading of the tracer, merging of the fingers, and diffusion across the finger interface. These sources of mixing affect the length and sharpness of the fluid-fluid interface, which are captured quantitatively by ϵ c . As mentioned in Section 3.2, ϵ c is a dynamic quantity itself; spreading of the slug and entry of the trailing front act as sources of ϵ c , and finger coalescence and diffusion act as sinks of ϵ c . After breakthrough, the mixed fluid exits the domain, and the producer acts as a sink for both χ and ϵ c . The balance between the sources and sinks of mixing determines when χ starts to decrease in Figure 2.6(a). Entrance of the trailing front causes a spike in ϵ c because the concentration gradient is large across the front, which causes a sudden jump in ϵ c followed by a drop due to smearing of the front from diffusion across the front Figure 2.6(b). The next event is breakthrough after which ϵ c starts to drop monotonically. The breakthrough curve c prd (t) shows arrivals of the viscous fingers at the producer as wiggles. All the fingers exit the domain by t ≈ 2.4 and mixing in the domain is purely diffusive after that, which results in another sharp drop in the value of scalar dissipation rate (Figure 2.6b). The breakthrough timet br , thetimet mc atwhichc prd ismaximum, andtheconcentrationc prd (t mc )arethreeimportant quantities during a contaminant remediation operation because their values indicate when, and how much, remediation effort is necessary at a given location. In enhanced oil recovery applications, e.g. waterflooding or CO 2 injection, the same quantities determine the estimated ultimate recovery (EUR), and the life of the well before it becomes uneconomical due to excessive production of the injectant. Compared to a purely diffusive ( R = 0) mixing and spreading process, viscous fingering shortens t br and t mc , and increases c prd (t mc ) Nicolaides et al. (2015). We identify the coupling between transport and deformation processes by plotting mixing and 27 Figure 2.6: Spreading and mixing in the Base Case. (a) Degree of mixing χ evolves non- monotonically because of the balance between source and sink terms of the mixing process. The decline in χ begins shortly after breakthrough when the sink begins to dominate source. (b) The mean scalar dissipation rate ϵ c curve jumps when the trailing front enters the domain (source of concentration gradients), and it declines rapidly after t mc because viscous fingering has ceased and slug drains out (sink). The blue dashed curve is the analytical ϵ c (Eq. 2.23) for Fickian mixing in a periodic (no entry or breakthrough of slugs), homogeneous, viscously stable and geomechanically- uncoupled flow. The difference between analytical and numerical ϵ c are due to fingering-enhanced interface stretching during 0<t<0.7, entry of the slug’s trailing front at t≈ 0.7, and breakthrough att br =1.16. (c) Wiggles in the breakthrough curve, c prd , indicate arrivals of fingers at the producer. t mc is the time at which c prd is maximum. spreading metrics against deformation metrics. I plot the degree of mixing as a function of the root mean square (RMS) of the change in porosity, ∆ ϕ rms , and the breakthrough concentration as a function of the root mean square of the change in effective volumetric stress, ∆ σ ′ v,rms (Fig- ure 2.7). As the transport-deformation coupling becomes stronger, e.g. α increases or K decreases (from left-most curve to right-most curve in Figure 2.7), the slope of χ vs. ∆ ϕ rms changes. This signifies the dependence of the degree of mixing on the change in porosity due to different degrees of geomechanical coupling. Porosity is related to the volumetric strain of the reservoir, and the effective volumetric stress is related to the total volumetric stress and pressure. Together, they can characterize the state of volumetric deformation in the reservoir. Also, the RMS of change from initial porosity or stress is numerically comparable to χ because σ 2 c is defined as a mean of squared deviation (Equation 2.21). Based on Figure 2.4, I can also choose the RMS of change in the von Mises stress as a metric for shear deformation to investigate coupling between transport and induced shear failure events such as hydroshearing. 28 0 0.05 0.15 0.25 0.35 0 0.5 1 1.5 2 2.5 0 0.2 0.4 0.6 0.8 (0.5, 3.33) Base Case = (0.5, 6.66 GPa) (0.25, 3.33) (0.95, 3.33) (MPa) ×10 -3 0 2 4 6 a b time time Figure 2.7: Coupling between transport and deformation. (a) The degree of mixing χ and the root mean square of porosity change, ∆ ϕ rms , increase and decrease synchronously as the slug is transported through the deforming reservoir. (b) The breakthrough concentration c prd and the root mean square of effective volumetric stress change, ∆ σ ′ v,rms , are also related. Both relations are sensitive to the values of the Biot coefficient α and the drained modulus K. 2.8. Sensitivity Analysis Figure 2.7 already shows how the degree or strength of coupling varies as α andK are varied. Here, I analyze the sensitivity of the coupling mechanisms to various poroelastic parameters in detail. 2.8.1. Biot coefficient The Biot coefficient α captures the strength of cementation between the grains relative to the grain compressibility. I saw in the mathematical model section (Equation (2.9)) that for an incompressible system, or in the drained deformation case, α =δϕ/δϵ v , and the Biot coefficient can be interpreted as the change in porosity for a unit change in the volumetric strain. α also dictates the degree of coupling between flow and mechanics problems through the coupling terms in Equations (2.1) and (2.5). According to the compatibility relation of α =∇ 2 σ ′ v /∇ 2 p, α determines the deviation in σ ′ v from its local average value for a given value of deviation-from-local-average in pressure. In other words, for smaller values of α , the effective volumetric stress is smoother and more independent of pressure fluctuations, such as the fluctuations arising from fingering or permeability field. α is less than 1.0 for rocks with high degrees of cementation and close to 1.0 for soft soil. I vary it from 0.95 (loose soil) to 0.1 (rigid, tight rocks). Based on Equation 2.13, I hypothesize that the effect of α on spreading and mixing is limited by the maximum change in porosity, which is less than 1% in all the simulations considered here. 29 200 400 600 0.32 0.33 0.34 1.1 1.2 1.3 10 -1 10 0 10 1 0.28 0.3 0.32 0.34 0.36 1 1.05 1.1 1.15 1.2 10 -3 10 -2 10 -1 0.32 0.33 0.34 0.35 1 1.1 1.2 1.3 Pe 0 1 2 0.1 0.2 0.3 0.4 0.5 1 1.5 2 0.2 0.4 0.6 0.8 1 0.332 0.334 0.336 0.338 1.157 1.158 1.159 1.16 (1/GPa) 0 2 4 6 (GPa) 0.2 0.4 0.6 0.8 1 0.336 0.3364 0.3368 1.14 1.15 1.16 1.17 1.18 0.33 0.34 0.35 0.36 1.05 1.1 1.15 1.2 Figure 2.8: Sensitivity of the degree of mixing at breakthrough (χ br ) and the breakthrough time (t br ) to deformation parameters (upper row) and transport parameters (lower row). Figures 2.8 and 2.9 analyze the effect of α on slug spreading and mixing behavior. The main observation is that the effect is weak for the range of parameter values considered here. t br decreases with α and the degree of mixing at breakthrough, χ br , increases slightly with α . α ’s effect on breakthrough is more significant than its effect on mixing. The effect on spreading appears at late times after breakthrough (Figure 2.9). The reason is that the primary effect of α is on the finger velocity. As α increases and the fingers move faster, breakthrough occurs earlier. The effect of α on χ is not significant, because fingers do not spend enough time in my relatively short domain to cause a difference in mixing. To better observe the effect of poroelastic or geomechanical coupling on the spatial distribution of tracer, I plot the difference of concentration fields corresponding to different values of the Biot coefficient, which is a measure of the poroelastic coupling strength, in Figure 2.10. The fields are at respective breakthrough times. As α increases, I see that the concentration at the trailing front of the injected slug, i.e. in the lower part of the figure, increases. This occurs because the finger velocities are different for different α . More interestingly, I observe concentration differences in the upper part of the domain where fingering occurs. This difference becomes more pronounced as poroelastic coupling becomes stronger from left to right in Figure 2.10, which suggests that neglecting poroelasticity can lead to an erroneous concentration field. 30 1 1.5 2 2.5 3 0 0.2 0.4 0.6 0.8 1 = 0.25 0.5 0.95 t Figure 2.9: Effect of Biot’s coefficient α on spreading. The breakthrough concentration c prd vs. the dimensionless time for different values of α shows that a less consolidated (higher α ) formation has an earlier arrival of the maximum concentration (shorter t mc ) compared to a more unconsolidated one. c prd (t mc ) have similar values for different values of α . -0.1 -0.05 0 0.05 0.1 Figure 2.10: Effect of Biot’s coefficient α , a measure of poroelastic coupling strength, on the tracer concentration field. The concentration field corresponding to α = 0.25 (relatively low poroelastic coupling strength) is subtracted from the concentration field corresponding to α = 0.5 (left), 0.7 (middle),and0.9(right)attheirrespectivebreakthroughtimes. Astheporoelasticcouplingstrength increases,thedifferenceintheconcentrationfieldalsoincreases. SeeTable2.2forallotherparameter values. 31 Case σ 2 logk ξ x (m) ξ y (m) Geometric mean k g (md) Base Case 0.01 1 1 50 Case 1 0.001 1 1 50 Case 2 0.1 1 1 50 Case 3 0.01 10 1 50 Case 4 0.01 1 10 50 Table 2.3: Sensitivity cases for the permeability heterogeneity. (a) (b) Figure 2.11: (a) Heterogeneous permeability fields for different values of the log permeability vari- ance σ 2 logk . Color scale is log(k) with k in md. (b) Concentration fields at breakthrough for the three cases. Finger tips are more pointy (less rounded), split and diffused, and fingers are longer at highest σ 2 logk . Color scale is for dimensionless concentration. 2.8.2. Permeability heterogeneity We vary the log-permeability variance σ 2 logk and the ratio of correlation lengths ξ y /ξ x as shown in Table 2.3. σ 2 logk controls the range of permeability values in the domain. At higher variances e.g. Case 2, the more permeable regions of the reservoir will be connected by streamlines leading to channelingofthetracerandashorterbreakthroughtime, whereasthelesspermeableregionswillact as stagnation or hold-up zones leading to a long residence or slug removal time Chen and Meiburg (1998); Nicolaides et al. (2015). Permeability heterogeneity controls the spatial location of both the tracer and in-situ fluids equivalently (Equation (2.3)), and the viscosity contrast controls the growth rate of their interface Tan and Homsy (1992). This leads to rapid stretching of the interface and subsequent mixing. Figure 2.11 shows c(x,y,t br ) for the three values of σ 2 logk . I observe more spreading of the slug and, hence, expect an earlier breakthrough and higher degrees of mixing at higher heterogeneity levels. This is confirmed in Figure 2.8. Presence of stagnation zones at higher 32 ξ y /ξ x ξ y /ξ x Figure 2.12: (a) Heterogeneous permeability fields for different values of the correlation length ratio ξ y /ξ x . (b) Concentration at breakthrough for the three cases. Fingers are longer at larger ξ y /ξ x . σ 2 logk also contributes to larger spreading of the tracer slug and higher drainage time. Next, I vary the ratio of permeability correlation lengths ξ y /ξ x to evaluate the effect of channelized features on the transport characteristics. I expect that Case 4 with a longer correlation length along the flow direction will lead to tracer channeling and earlier breakthrough. Figure 2.12 shows that as ξ y /ξ x increases, the fingers are elongated in the y-direction and thinner in the x-direction. Figure 2.8 verifies the hypothesis of earlier breakthrough at ξ y /ξ x > 1. Mixing of the slug is also enhanced across the longer interface achieved with larger ξ y /ξ x . 2.8.3. Peclet number Peclet number is the relative strength of advection compared to diffusion. Here, I consider the range of 200≤ Pe≤ 700. A lower value of Pe, due to lower injection rate q in or higher diffusion-dispersion coefficient D d , leads to diffused fingers that are wider and fewer in number compared to another case with higher Pe. Higher Pe also means faster travelling fingers (higher growth rate of the viscous instability Homsy (1987)), i.e. earlier breakthrough as seen in Figure 2.13 and confirmed in Figure2.8. Extensivetip-splittingandmergingoffingersleadtoenhancedspreadingandlower σ 2 c or higher χ . Merging of fingers increases the advective energy of the dominant finger to further spread the solute in the domain. This balances against the lower diffusive flux at higher Pe, which causes loss of diffusive mixing and lower χ . The balance between the two mechanisms of splitting-merging and diffusion results in a non-monotonic profile of χ br vs. Pe as shown in Figure 2.8. 33 Figure 2.13: Concentration fields at breakthrough for different values of the Peclet number. Pe=700 showsmoretip-splittingandlessdiffusion(sharperfingers)comparedtoPe=200and400. Compared to Pe=200, the trailing front of the slug in Pe=400 is farther from the leading front, which means more spreading and faster breakthrough. Also, Pe=400 shows more fingers compared to Pe=200. Figure 2.14: Concentration fields for different cases of the log viscosity ratio R at the respective breakthrough times. The breakthrough, which is indicated by the position of the trailing front, occurs earlier at higher R because of faster traveling fingers. 2.8.4. Viscosity ratio The log viscosity ratioR between the injected and in-situ fluids determines the number, size (width) and speed of the fingers Homsy (1987). Since fingering plays a critical role in the development of miscibilitybetweenthetwofluidsandinrecoveringthein-situfluid(incaseofmiscibleCO 2 injection in EOR) and the tracer fluid (in case of contaminant remediation), R is an important parameter for sensitivity analysis. When I inject a less viscous fluid into the aquifer ( R > 0), the degree of mixing at breakthrough can be as much as 88% higher compared to injecting a tracer with R = 0 (Figure 2.8). This is due to faster tip velocities of the fingers, which lead to faster spreading of the slug. As fingers develop, dispersion or spreading of the solute becomes stronger, resulting in sharper concentration gradients and a longer interface, which in turn augments the diffusive flux and the degree of mixing. 34 2.9. Summary and Conclusions We investigated the coupling between poroelastic deformation and solute spreading and mixing during viscous fingering using high-resolution numerical simulations. Motivated by the difference in the mathematical nature of the governing equations, I proposed a numerical formulation of a monolithicflow-geomechanicssolversequentiallycoupledwithatransportsolverviathevelocityfield and the concentration-dependent mobility model. Transport-deformation coupling mechanisms and their sensitivity to poromechanical parameters were identified and analyzed. A sensitivity analysis was performed to obtain the extent of poroelastic influence on the macroscopic transport quantities: the breakthrough time, the breakthrough curve, and the degree of mixing. The sensitivity analysis considered the heterogeneity in permeability field, the Biot coefficient, log viscosity ratio R and Peclet number. The parameters and their values were chosen to represent injection in real-world confined aquifers. The results show that the poroelastic processes influence solute transport. From Fig. 2.8, I observe that spreading of the injected tracer slug is accelerated, i.e. breakthrough occurs earlier, when the rock is less consolidated or more loosely packed, viscosity contrast increases, or diffusion is weaker compared to advection. Mixing of the injected slug is enhanced when rock becomes less consolidated and more deformable. Failure to account for this coupling will result in an inaccurate prediction of the breakthrough time, solute recovery, and the degree of mixing, thus an inaccurate determination of the efficiency of the groundwater remediation or enhanced oil recovery operation. The 2D study presented here can be extended to 3D. The physical mechanisms of advection, diffu- sion, and pressure-induced deformation will carry over from 2D to 3D, however, I may observe the following changes. The flow field around the wells will change from linear to radial resulting in a 3D stress field, which will not satisfy the plane stress approximation. Including the vertical dimension to the domain will activate gravity, which appears as a body force in the mechanics problem and as a velocity potential in the flow problem. The effect of gravity as a body force can be absorbed into the initial stresses because I assume the bulk density to remain constant in time. Also, since I assume a slightly compressible fluid with a density that is independent of the tracer concentration, 35 spatial variations of the fluid density can be neglected. Therefore, I believe my main result related to transport-deformation mechanisms are generalizable. Our workflow employs different discretization methods for different sub-problems: finite element method for the mechanics problem, finite volume method for the flow problem, and compact finite difference method for the transport problem. Finite element and finite volume methods can be applied easily to unstructured grids, because the volume and surface integrals required for terms in Eq. (2.25) can be evaluated for arbitrary shapes and orientations. However, the compact finite difference method is more appropriate for a structured grid because I assume that the grid lines are oriented parallel to the two Cartesian directions of my 2D domain with a fixed number of cells in each grid direction. In other words, fixed values of N x and N y in x and y directions are required to build the differentiation matrices and discretize the spatial derivatives in the transport solver. To extend my work to unstructured grids, one may consider using a compact high order finite volume method Ollivier-Gooch and Altena (2002) to solve the transport problem. 36 CHAPTER 3 Effect of Poroelastic Coupling and Fracture Dynamics on Solute Transport and Geomechanical Stability 3.1. Abstract The coupling between solute spreading and mixing mechanisms and geomechanical processes, e.g. stress-induced fracture activation, has emerged as an important hydrogeological challenge due to its role in applications such as underground waste disposal, carbon sequestration, contaminant reme- diation, enhanced oil recovery, and tracer-based reservoir surveillance. Despite recent advances in modelingflowandgeomechanicscoupling, aholisticapproachtocapturingthesynergybetweenfluid flow, solute transport, induced stresses, and fracture mechanics is lacking. This study investigates the rich interplay between these processes within a novel computational framework that is proposed to solve the coupled flow, transport, geomechanics, and fracture mechanics problem. Embedded Discrete Fracture Modeling (EDFM) is used to model the flow and transport processes in fractured porous media while an improved Bandis model is employed to capture the fracture-mechanical response to flow-induced stress perturbations. The role of transport-geomechanics coupling in mod- ulating the spreading and miscibility of a solute slug during viscously unstable flows is examined. I investigate how flow-transport coupling, parameterized through the solute viscosity contrast and the fracture permeability, influences the stress state and fracture stability in the domain. A case study, inspired by a huff-n-puff tracer flowback study, is conducted to investigate the applicability of the proposed framework in the field. A sensitivity analysis is performed to evaluate the dependence of global transport characteristics, permeability evolution, and fracture stability on parameters dictating the strength of coupling between geomechanics, flow, and transport. 3.2. Introduction Injectionoffluidsintonaturallyorhydraulicallyfracturedformationshasbeenanimportantresearch topicduetoitsrelevanceinavastnumberofengineeringapplications: wastedisposalWitherspoon et al. (1981);McCarthy and Zachara(1989);Cornaton et al.(2008),carbonsequestrationIding and Ringrose 37 (2010), contaminant transport Sahimi (2011), enhanced oil recovery Jiménez-Martínez et al. (2016), and tracer surveillance Hu and Moran (2005); Rugh and Burbey (2008); Warner et al. (2014). In many applications, a key objective of the numerical model built to simulate the processes is to quantify the mixing between the injected and resident fluids and quantify the extent of the asso- ciated mixing zone Dentz et al. (2011); Jha et al. (2011a); Cirpka and Valocchi (2007); Zhao et al. (2011); Bonazzi et al. (2020,0). In cyclic well operations, where the well alternates or cycles through injection and production/withdrawal stages, the evolution of fluid mixing is complicated by short- time scale variations in the pore pressure field. Often, the reservoir of interest is geomechanically sensitive and deformable due to the presence of pre-existing fractures, unconsolidated formations, or excessive flow-induced stresses. Resolving the geomechanical coupling effect in those reservoirs becomes critically dependent on modeling the effect of stress-modulated processes such as fracture opening and shear failure on mass transport processes such as advection and diffusion Matthäi et al. (2010); Fox et al. (2011b); Latham et al. (2013); Yan et al. (2019); Tran and Jha (2020).This is re- quired to answer questions such as how much change in the breakthrough time or degree of mixing corresponds to a given change in the average volumetric stress during groundwater withdrawal. The presence of hydraulically stimulated and natural fractures in a low permeability rock adds an extra dimension of complexity because the characteristic length and time scales of flow and defor- mation processes can be markedly different between the fracture and the host rock domains. Often understood as displacement discontinuities, fractures have been proven to make a first-order impact on the flow behavior of any rock Levison and Novakowski (2012); Hardebol et al. (2015) by alter- ing its permeability. While accounting for all individual fractures in a reservoir is computationally intractable and unnecessary, the contribution of a subset of fractures, at length scales relevant to the dominant physical processes in the problem, must be assessed. The role of flow-geomechanical coupling in modulating hydraulic properties such as stress-dependent permeability and porosity is increasinglyrecognizedasanecessityinreservoir-scalesimulationmodelsRutqvist and Stephansson (2003); Bubshait and Jha (2019); Meguerdijian and Jha (2021). Flow processes can also change the spatial and temporal evolution of stress state and fracture dynamics, which affect the mechanical stability of fractures and formation integrity Jin and Zoback (2018). Therefore, a tightly coupled 38 flow and geomechanics computational framework Tran and Jha (2020); Zhao and Jha (2019,0) is of- ten required to correctly encapsulate and explain the interaction among these underlying processes. This can help us find geomechanical signals, e.g. change in the magnitude and direction of principal stresses or activation/healing of fractures oriented along a certain direction, that I can monitor in the field to infer the plume size and dilution during injection. Modeling of discontinuities such as fractures usually falls under two popular categories of meth- ods: discrete fracture network modeling and continuum modeling (such as the continuum damage mechanics approach Pogacnik et al. (2016); Yan et al. (2019); Bubshait and Jha (2019)). Methods in each category have their own merits and drawbacks, most of which can be understood in terms of a method’s sensitivity to mesh size dependency, the capability to model localization of flow or deformation and the resulting anisotropic response, and the computational cost to achieve a certain order of numerical accuracy. The cost, which is related to the number of iterations required to solve the discretized linearized system, is often expressed in terms of the number of degrees of freedom and/or the condition number of the Jacobian matrix of the discrete system. The EDFM approach capitalizes on the best features of discrete and continuum approaches to model the effect of fractures while maintaining a relatively lower cost. The basic idea of EDFM is to represent fractures using a structured gridding scheme, separate from the host matrix grid, and use special non-neighboring connections to model the fluxes between fracture and matrix Lee et al. (2001a); Hajibeygi et al. (2011). EDFM allows a numerically fast and stable representation of individual fractures’ contribu- tion to flow and transport processes. Compared to traditional fracture modeling methods such as dual-continuum modeling, local grid refinement, and unstructured gridding, EDFM offers a balance between accuracy, flexibility, gridding, and computational efficiency. To improve cross-media flux transfer in media with high conductivity contrasts, e.g. in hydraulically fractured shales or geother- mal rocks, 1392017Tene et al.Tene, Bosma, Al Kobaisi, and Hajibeygi () proposed pEDFM with adjustmentsofmatrix-matrixandfracture-matrixtransmissibilitiesintheimmediatesurroundingof the explicit fractures. Compartmental EDFM (c-EDFM) allows sub-grid resolution by splitting the matrix grid when a fracture grid intersects a matrix cell Chai et al. (2018). 1282018Sangnimnuan etal.Sangnimnuan,Li,andWu()developedacoupledflow-geomechanicsmodelwithEDFMtochar- 39 acterize stress evolution during fluid production from ultra-low permeability reservoirs and to model stressredistributionandreorientationinducedbythedepletion. 1092016Norbecketal.Norbeck, Mc- Clure, Lo, and Horne () addresses how to include nonlinear evolution of nucleated tensile fracture and shear failure events. 1032013Moinfar et al.Moinfar, Sepehrnoori, Johns, and Varavei () models the change in aperture and associated permeability of an evolving fracture network via the effective normal stresses acting on fracture planes within an EDFM framework. Despite a recent shift of focus to coupled flow-geomechanical modeling, there have been limited attempts to model transport phenomena in a complex fracture network under the framework of poroelasticity. The effect of geomechanical coupling on solute spreading, mixing, and viscous fin- gering in stress-sensitive reservoirs has been explored recently Tran et al. (2018); Tran and Jha (2020); Sweeney and Hyman (2020). In addition, recent advances in understanding fracture dy- namics have shed new light on inherently strong nonlinearities that exist between the state of stress and closure and dilation of natural joints Olsson and Barton (2001); Li et al. (2020). Often, these nonlinearities lead to hysteresis in the mechanical response of the jointed rock when subjected to loading/unloading cycles over time scales of interest. Hysteresis in stress-strain behavior, post-peak shear strength, shear induced dilation are required to improve fracture modeling in aquifers deform- ingunderinjection/productionAsadollahi and Tonon(2010). Yet,manymodern-daygeomechanical simulators neglect this behavior because of the lack of a computationally efficient poroelastic model for fractured rocks. This paper attempts to bridge these technical gaps by proposing modifications in the classical Bandis model Bandis et al. (1983), implementing it within a coupled flow-transport- poroelastic framework, and successfully demonstrating its capability in predicting the following observables during tracer tests in aquifers with complex fracture networks: the tracer breakthrough profile, pore pressure field, and aquifer stress and strain fields around the well. The physical properties of a solute, which can be a contaminant or a remediation agent, can vary in space and time due to the dependence of solute density, viscosity, or diffusivity on the concen- trations of the dissolved species Flowers and Hunt (2007). Further, the contaminant viscosity may change with time due to phase separation, dissolution, and evaporation of lighter volatile com- 40 Chemical Viscosity (cp) R=ln µ 0 µ 1 Benzene 0.6 0.39 Toluene 0.56 0.46 Ethylbenzene 0.64 0.33 O-Xylene 0.76 0.16 Trichloroethylene (TCE) 0.55 0.48 No. 4 fuel oil 36 -3.7 Trans-1,2dichloroethylene 0.4 0.8 Chlorobenzene 0.8 0.11 m-Cresol (pesticides, antiseptics) 6.1 -1.92 Diesel 2.2 -0.9 MTBE (gasoline additive) 0.37 0.88 Propane 0.1 2.19 Table 3.1: Groundwater chemicals with viscosity µ 1 different from the water viscosity µ 0 . A water viscosity of 0.89 cp at 25 degC is assumed. ponents Mercer and Cohen (1990). Contaminant remediation processes often alter the dynamic properties of contaminant plumes. More specifically, remediation by injection of steam, surfactants, oxidizing agents, or polymer can produce a mobile, dissolved-phase plume with a viscosity that is different from the viscosity of the fluid bank ahead or behind it. Table 3.1 displays the viscosity and log-viscosity ratio R of some groundwater chemicals compared to water. Even for non-aqueous phase liquid (NAPL) contaminants, which are considered mostly immiscible with water, miscible plumes are created as components partition from the NAPL phase into the flowing groundwa- ter Poulsen and Kueper (1992). This has two major implications for contaminant transport. First, changes in the viscosity can affect contaminant mobility even in homogeneous permeability fields. Second, when a less viscous fluid is displacing a more viscous fluid, the transport becomes hydrody- namically unstable and culminates into viscous fingering Homsy (1987), which leads to stretching of the plume interface and enhanced mixing Jha et al. (2011b); Nicolaides et al. (2015), which then affects contaminant recovery. To answer questions such as whether solute fingering leads to more dilution in a stiffer rock with fracture activation or in a softer, intact rock with large porosity changes, I need to understand the role of transport-geomechanics coupling in advection-dominated transports. 41 Below I present my modeling framework in the context of a push-chase-pull tracer test and use the model to elucidate physical mechanisms that arise solely due to the two-way coupling between flow, transport, and poroelasticity. The novelty of my work lies in unraveling these coupling mechanisms and exploring their relationships during hydrodynamically unstable transports in fractured porous media. In particular, I focus on the mechanisms that affect the mixing and spreading of the in- jected tracer, because these processes control the effectiveness of contaminant removal/remediation in groundwater applications. Finally, I use the framework to analyze the sensitivity of the field ob- servables to changes in various rock-fluid properties, which is essential for the proposed framework to be applicable to subsurface problems where properties may vary over large and/or uncertain ranges. 3.3. Physical Model We consider a horizontal section of an aquifer with a horizontal well in the middle (Figure 3.1) that injects and produces alternatingly. The reservoir is fully saturated with a single-phase fluid (water) and is under hydromechanical equilibrium with external stresses S h and S H , which are the minimum and maximum principal horizontal compressions, respectively, from regional tectonics. The stochastic fracture network depicted there is inspired by a tracer flowback study Li et al. (2016) and is generated using the FracGen simulator based on fracture density, formation resistivity-based imaging logs, and the tracer breakthrough curve Myshakin et al. (2015); Li et al. (2016). This fracture network is designed to encompass complexities such as fracture-to-fracture intersections, well-to-fracture crossings, and heterogeneity in fracture orientation with respect to theS H direction. The cyclic well first acts as an injector, injecting a slug of tracer that is perfectly miscible with the in-situfluidbuthasalowerviscosity. Thetopandbottomboundariesareopentodrainageatafixed pressure that is lower than the initial aquifer pore pressure. As the less viscous slug moves through theaquiferanddisplacesthemoreviscousin-situfluid, theSaffman-Taylorhydrodynamicinstability initiates at the tracer front giving rise to viscous fingering. The strength of fingering depends on the viscosity contrast between tracer and water, the tracer injection rate, tracer diffusivity, and the characteristics of permeability field Tran and Jha (2020), which in this model is dominated by the fracture network. I stop the injection of tracer when it breaks through at either of the two 42 boundaries. Following breakthrough, the well switches to injection of the in-situ fluid (zero tracer concentration) to push the tracer further into the aquifer, and I call this duration the tracer soaking period. The soaking period lasts 30 days and is followed by the tracer extraction or withdrawal period, which is 90 days in the model. During withdrawal, the well produces both the tracer and the in-situ water. This process resembles either a typical huff-n-puff operation in enhanced oil recovery (EOR) or a push-and-pull tracer test in groundwater surveillance. In those applications, the observed pressure change during the injection period is often in the range of 10-30 MPa Li et al. (2016). The physical model is designed to investigate coupling between geomechanics and transport in fractured porous media during both injection and withdrawal episodes. Since the pore pressure is coupled to both the tracer velocity and the effective stress, I can use this model to probe the coupling between transport and deformation processes to extract possible physical mechanisms arising from such coupling. This will allow us to address questions such as how poroelasticity affects the mixing and spreading of solute transport and how flow-transport coupling, manifested through viscous fingering, impacts stress state and fracture stability. I use the same model to examine how stress changes dictate permeability evolution and how transport metrics are modulated by fracture mechanics. We make the following assumptions, which are appropriate for the physical problem at hand: (1) slightly compressible fluid, (2) small strain, (3) quasi-static equilibrium, (4) linear elasticity, (5) single-phase flow, (6) isotropic material behavior, (7) isothermal system, and (9) first contact mis- cibility between the fluids. Given my focus on groundwater transport problems with low values of induced stresses compared to, e.g., hydraulic fracturing problems, I make additional assumptions regarding the fracture behavior: (10) mechanical equilibrium in the domain is governed by the deformation of matrix only, and the fracture deformation is limited to stress-dependent changes in the fracture aperture. (11) the gradient of fracture pore pressure normal to the fracture is neg- ligible. Assumption (10) suggests that volumetric strain occurs only in the host matrix and not in the fractures. Thus, the fracture pressure is coupled only to the matrix pressure perturbation 43 Cyclic Well S H S h Figure 3.1: A physical model of tracer injection, tracer soaking or chase fluid injection, and tracer- water withdrawal in a reservoir populated with discrete natural fractures of different lengths and orientations. The horizontal well is intersected by some of the fractures. The initial state of stress in the reservoir is determined by boundary compressions S h and S H and zero normal displacement roller boundary conditions. The fractures are denoted by the green lines. The cyclic well is denoted in the thick blue line. The yellow-shaded region around the well indicates the mixture of injected tracer (viscosityµ 1 ) and resident water (viscosityµ 0 ), which forms a short time after tracer injection begins at the well. The remaining black domain is saturated with the resident fluid. and not coupled directly to the mechanical deformation of the matrix. Input parameters such as the well injection rate and the matrix permeability are selected such that the likelihood of fracture propagation or nucleation is negligible. Therefore, in this study, dynamic propagation of fracture is not considered and stress-dependent changes in fracture aperture are modeled via constitutive equations such as the Bandis model Bandis et al. (1983). As the fractures are generally very thin and highly permeable compared to the surrounding host rock, assumption (11) is generally valid. 3.4. Mathematical Model We consider single-phase flow and tracer transport through a poroelastic and deformable solid in two dimensions. Below, I present the governing equations in the strong form for both fracture and matrix domains. The equations will be converted into the respective weak forms and discretized. I adoptacontinuumrepresentationofamacroscopicmodeloftheporoelasticmedium, wherefluidand solid are viewed as overlapping continua Coussy (2004). Equations of quasi-static equilibrium, solid and fluid mass conservation, and advective-diffusive transport govern the system’s behavior in space 44 and time under the imposed initial and boundary conditions as well as the applicable constitutive laws for the solid, fluid, and tracer components. I will cast the problem into a formulation where the primary unknowns are as follows: the excess pressure in the matrix domain, the excess pressure in the fracture domain, the excess displacement in the matrix, the normalized concentration of tracer in the matrix, and the normalized concentration of tracer in the fracture domain. Here, ‘excess’ is defined as the difference between a quantity’s current value and its reference value; see section 3.4.6. 3.4.1. Fluid Mass Conservation The 2-D model in the Cartesian domain has an area Ω = ( L x × L y ), where L x and L y are the dimensions in the x and y directions. I denote the closed boundary of the domain by Γ . The initial matrix porosity is ϕ 0 , and the matrix permeability k m is assumed to be isotropic and constant in time. Subscript 0 indicates the initial condition and superscripts m, w and f denote matrix, well, andfracturedomainsrespectively. Theinitialfracturepermeabilityk f 0 isthesameforallfracturesso that the dependence of permeability on stress can be compared easily. The domain is saturated with a fluid of density ρ f , compressibilityC f , and viscosityµ f . The fluid mass conservation equations for single-phase slightly compressible flow under infinitesimal deformation in the matrix and fracture domains are 1 M m ∂p m ∂t +α m ∂ϵ v ∂t +∇· v m =Ψ mf +Ψ mw , (3.1) 1 M f ∂p f ∂t +∇ f · v f =Ψ fm +Ψ ff +Ψ fw , (3.2) where p is the pore fluid pressure ( p m in matrix and p f in fracture), t is time, 1/M =ϕ 0 C f +(α − ϕ 0 )/K s is the inverse Biot modulus, α =1− (K dr /K s ) is the Biot coefficient indicating the degree of cementation between grains, K dr is the drained bulk modulus of rock, K s is the bulk modulus of solid grains, C f =(1/ρ f )dρ f /dp is the fluid compressibility from the equation of state, v is the Darcy fluid velocity vector relative to the solid skeleton, and ϵ v is the volumetric strain. The Ψ terms indicate the mass flux transfer from fracture to matrix (superscript mf), matrix to fracture (superscript fm), fracture to fracture (superscript ff), well to matrix (superscript mw), and well to fracture (superscript fw). ∇ f is the nabla vector operator in the fracture domain, i.e. along the 1D fracture mesh that is embedded in a 2D matrix mesh. So, the tangential divergence of an 45 arbitrary vector quantity κ , i.e. ∇ f · κ , in the fracture domain can be related to the divergence ∇· κ in the matrix domain via the projection operator as follows: ∇ f · κ = I d − (n f ⊗ n f ) :∇· κ , wheren f is the fracture surface unit normal and I d is the identity tensor. The tangential gradient operator can be defined similarly. We assume that the Biot coefficient and porosity values are the same for matrix and fracture domains. Viscosity and density of the pore fluid are also assumed to be the same in the fracture and matrix domains. The main distinction between fracture and matrix domains lies in the permeability contrast between the two media. Ignoring gravity and buoyancy effects, the Darcy velocities in the matrix and fracture domains are v m = k m µ f (c m ) (−∇ p m ), (3.3) v f = k f µ f (c f ) (−∇ f p f ), (3.4) where µ f (c) is the concentration-dependent fluid viscosity. The fluid pressure experiences pertur- bations that stem from pressures at the domain boundaries, fluid flux (mass rate per unit area) across interfaces between the three domains (fracture, well, and matrix), fluid source/sink terms at the center well (mass rate per unit volume), and the volumetric strain rate which changes the fluid storage capacity of the medium. Below I discuss the closure relations for the Ψ functions. 3.4.2. EDFM Formulation Fracture-Matrix Intersection In EDFM formulation, pressures in the fracture and matrix domains are coupled to each other via the mass transfer functions. For closure purposes, these terms need to be expressed in terms of the pressures in the two domains. I write the matrix-to-fracture transfer function Ψ as Ψ fm =CI fm λ fm (p m − p f ), (3.5) where λ fm is the fluid mobility calculated as the harmonic average of the matrix cell mobility and the fracture segment mobility at their interface, CI fm is the connectivity index between a matrix cell and a fracture segment, and it represents the geometrical complexity of this intersection. The 46 fracture-to-matrix flux is assumed to be equal and opposite in sign: Ψ mf =− Ψ fm . After both the matrix and the fractures are meshed (see the numerical model section below), the discrete form of the connection mobility is λ fm i,k = 2λ m i λ f k λ m i +λ f k (3.6) and the connectivity index is CI fm i,k = D i,k ⟨d⟩ i,k , (3.7) where subscript i denotes a matrix cell, subscript k denotes a fracture segment, D i,k is the length fraction of fracture segment k inside matrix cell i,⟨d⟩ i,k is the calculated distance between matrix cell and fracture segment. This approach to couple fracture and matrix guarantees that the total flux between matrix and fracture is conserved. Fracture-Fracture Intersection Fracture-fracture interaction is captured via the transfer function Ψ ff . The flux between two fracture segments k andl, either belonging to the same fracture or different fractures that intersect, is defined as: Ψ ff kl =λ ff kl (p f l − p f k ) (3.8) λ ff kl is the harmonic average of the fluid mobilities in fracture segments k and l and it is calculated similar to Eq. (3.6). The fracture mobilities are defined similar to 1032013Moinfar et al.Moinfar, Sepehrnoori, Johns, and Varavei (): λ f k = k f w f µ f ∆ x f = k f √ 12k f µ f ∆ x f (3.9) Here, w f is the fracture aperture, which is calculated from the fracture permeability k f using the cubic law. ∆ x f is the length of the fracture segment along the fracture tangent direction. Fracture-Well Intersection The flux terms Ψ mw and Ψ fw represent fluxes from well-to-matrix and well-to-fracture domains, respectively. Well-intersecting fractures are expected to impact the flow significantly. I define Ψ fw 47 using the Peaceman well model Peaceman (1978) Ψ fw =WI fw (p w − p f ), (3.10) where the well index WI fw is calculated as in 1032013Moinfar et al.Moinfar, Sepehrnoori, Johns, and Varavei () WI fw = 2πk f w f µ f ln 0.14 q (L 2 f +h 2 f ) rw ! +s ! (3.11) In this equation, L f is the fracture length in the cell, w f is the fracture aperture, h f is the fracture height (assumed to be equal to the thickness of the horizontal aquifer in the z-direction, r w is the wellbore radius, and s is the wellbore skin due to near-wellbore effects such as precipitation, formation damage, non-Darcy pressure drop, etc. 3.4.3. Linear Momentum Balance The linear momentum balance equations for a fluid-saturated rock under quasi-static equilibrium and in absence of body forces are ∇· σ m =0, (3.12) where σ m = σ ′ ,m − α m p m I d is the Cauchy total stress tensor, and σ ′ ,m is the effective stress tensor. The bulk density of the matrix can be written as ρ b = ϕ m ρ f +(1− ϕ m )ρ s , where ρ f is the fluid density, ρ s is the solid grain density, and ϕ m is the Eulerian porosity of the matrix in the current configuration. A sign convention where the normal stresses are positive in tension and negative in compression is used. The 2D domain represents a thin layer of rock, thus a state of plane stress is assumed. In the compact engineering notation, σ = [σ xx ,σ yy ,σ xy ] T , ϵ = [ϵ xx ,ϵ yy ,2ϵ xy ] T , and I d = [1,1,0] T are the stress, strain and identity vectors in 2D, respectively. Substituting the effective stress into Equation (3.12), I obtain ∇· σ ′ ,m − α m ∇p m =0. (3.13) 48 Solving the flow and the geomechanics equations simultaneously implies a full monolithic coupling between the two sub-problems. In the mechanics sub-problem, the pore pressure gradient acts as an equivalent body force that drives changes to the rock deformation and stress state. In the flow sub-problem, the matrix volumetric strain rate acts as an equivalent fluid source, or fluid storage capacity change, causing changes in the pore pressure. Assuming a linear elastic behavior of the rock,σ ′ =Dϵ , whereϵ =(∇u+∇ T u)/2 is the infinitesimal strain tensor defined as the symmetric gradient of the displacement vector. The volumetric effective stress is related to the volumetric strain as σ ′ v =K dr ϵ v , where K dr =E/(2(1+ν )(1− 2ν )) in 2D. 3.4.4. Advective Dispersive Transport The Advection Dispersion Equations (ADE) for modeling tracer transport through matrix and fracture domains are ∂(ϕ m c m ) ∂t +∇· (v m c m − ϕ m D m d ∇c m )=χ mf +χ mw (3.14) ∂(ϕ f c f ) ∂t +∇ f · v f c f − ϕ f D f d ∇c f =χ ff +χ fm +χ fw (3.15) where c m and c f are the dimensionless normalized local concentrations of the tracer in the matrix and the fracture, and D d is the hydrodynamic dispersion coefficient assumed to be a constant. χ representsdifferenttransportfluxesacrossinterfacesbetweenthefracture, matrix, andwelldomains. During the injection episode, the well acts as a tracer source. I use an exponential model for the concentration-dependent viscosity of the mixture, µ f (c)=µ 1 e R(1− c) , (3.16) where c = c m in the matrix and c = c f in a fracture, and R = log e µ 0 µ 1 is the log viscosity ratio of the in-situ fluid viscosity µ 0 and the injected tracer viscosity µ 1 . Whenc m orc f is between 0 and 1, the viscosity of the mixture is higher than the viscosity of injected tracer and lower than the in-situ fluid viscosity, for positive values of R. Here, I assume that the mixture density is independent of the tracer concentration. 49 3.4.5. Fracture Mechanics Model The Bandis model is a widely used empirical model to provide a constitutive relationship between the state of stress and rock joint closure, which can be related to fracture hydraulic properties, i.e. permeability Barton et al. (1985). The Bandis model is based on a comprehensive laboratory investigation of the deformation behavior of a wide variety of natural unfilled joints under different loading/unloading and repeated load cycle conditions Bandis et al. (1983). The model is especially accurate for rough joints and fractures for which simpler models that assume smooth contacts often overestimatetheshearstrengthBarton and Choubey(1977). Foreachfracturesegmentwithnormal directionn f , the effective normal stress σ ′ n =σ ′ n f · n f , shear stress τ =|σ ′ n f − σ ′ n n f |, and shear slip displacement δ h are critical inputs to calculate the mechanical aperture of the segment. These equations are used to dynamically update the normal and shear stresses acting on the fracture segments using the matrix stress tensor output of the geomechanics solver and the known surface normal vectors of the segments. The slip displacement on a fracture segment is calculated using the displacement vectors of neighboring matrix nodes on either side of the fracture segment. The fracture permeability is updated based on the cubic law for flow between two parallel plates, which yields the following update relation k f =k f r a h a hr 2 (3.17) wherea h is the hydraulic aperture, and subscriptr indicates the reference value. The reference time is taken to be the end of the initialization period of a simulation (see initialization details below). Therefore, the reference permeabilityk f r is a function of the initial permeabilityk f i and the reference stresses at the end of initialization. Hydraulic aperture depends on the mechanical aperture and the ratio of shear displacement to peak shear displacement of each fracture segment Barton et al. (1985); Olsson and Barton (2001). Each fracture segment can cross multiple matrix cells, resulting in different aperture changes along its length because of the spatially changing stress tensor and slip vector corresponding to the matrix cells that the fracture segment intersects. The mechanical aperture of a fracture segment is the square root of the harmonic mean of the squared aperture of 50 its individual parts. The mechanical aperture of a fracture segment is measured by: a m =a 0 +∆ a n +∆ a s (3.18) where a 0 is the initial mechanical aperture under zero excess effective normal stress, ∆ a n is the change in mechanical aperture due to normal stress perturbation, and ∆ a s is the change in mechan- ical aperture due to shear-induced dilation. Hysteresis in the stress-deformation behavior has not been properly accounted for in existing models. Here, I implement a hysteretic stress-deformation model that captures the difference in fracture aperture change between the loading and unloading stages of rock under compression. This model also accounts for an expected fracture permeability decrease when the cyclic well changes from injection to withdrawal, shown by lines with circle and triangle markers respectively in Figure 3.2). Injection corresponds to decreased compression and henceunloadingoftherock, andproductioncorrespondstoloading. Basedonthefoundationalwork by 141985Barton et al.Barton, Bandis, and Bakhtar (), the dependence of mechanical aperture on normal stress change is calculated as: ∆ a n =                σ ′ n K ni if σ ′ n >0 (tensile state) 1 Vm− V i − K ni σ ′ n − 1 +V i if σ ′ n <0 (compressive state) and δσ ′ n >0 (unloading) 1 Vm − K ni σ ′ n − 1 if σ ′ n <0 (compressive state) and δσ ′ n <0 (loading) (3.19) where V i is the irrecoverable closure that contributes to the hysteretic behavior and δσ ′ n is the in- cremental normal effective stress between consecutive time steps. K ni ,V m , andV i are defined in the caption of Figure 3.2) and readily calculated using empirical correlations in 131982Barton (). my implementation has notable improvements compared to the original model in evaluating peak shear displacement (δ peak ), post-peak shear strength (B-C in the top right quadrant of Figure 3.2), aper- ture closure at small shear displacement, and fracture surface degradation. The model parameters are functions of the input parameters of the original Bandis model, which are the Joint Roughness Coefficient JRC, Joint Compressive Strength JCS, and Unconfined Compressive Strength UCS. JRC determines the shear strength of unfilled hard rock joints. It can be estimated by comparing the 51 appearance of a discontinuity surface with standard profiles published in 131982Barton (). JCS signifies the degree of weathering of the rock joint and can be estimated by the use of the Schmidt rebound hammer test. If the joints are completely unweathered, then JCS will be equal to the UCS of the unweathered rock. As the joints get weathered and weakened, JCS decreases accordingly. In addition to modeling the aperture closure dependency on the fracture normal stress, I also incor- porate a shear-related physical mechanism that alters hydraulic aperture. I take into consideration the reduction of post-peak shear strength from fracture surface degradation and shear-induced di- lation. The shear-induced dilation is formulated as a function of δ h Asadollahi and Tonon (2010): ∆ a s =        2δ 2 h 3δ peak − δ h 3 tan h JRC· log JCS σ ′ n i if δ h <δ peak ∆ a v,peak +δ peak R δ h /δ peak 1 tan h JRC p · log JCS σ ′ n · x − 0.381 i dx if δ h >δ peak (3.20) where x and dx are dummy variables of integration, δ peak is the peak shear displacement, ∆ a v,peak is the peak shear-induced dilation, JRC p is the peak joint roughness coefficient, and log refers to base-10 logarithm. These parameters are readily calculated based on the correlations given in52010AsadollahiandTonon(). Theshear-inducedeffectisintegratedduringtheinjectionepisodes only, with the parameters tuned to yield a reasonable relative change in permeability based on the experimental work 1562018Ye and Ghassemi (). Fractures tend to slip during the injection period when the effective normal compression acting on a fracture decreases. One of the limitations of my model is a lack of feedback from fracture slip to stress relaxation, i.e. the coupling between fracture dynamics and mechanical equilibrium (Eq. (3.13)) is one-way only (Figure 3.3). This can be justified based on the computational expense of solving the nonlinear contact problem on each of the fractures, many of which intersect each other, coupled to the equilibrium equations. One-way coupling has been used successfully in other transport studies, e.g. Sweeney and Hyman (2020). The coupling from fracture slip to domain-wide stresses can be important for long fractures slipping across their entire length and for propagating fractures, where an accurate calculation of induced seismic event time and magnitude may require an estimation of stress relaxation Jha and Juanes (2014); Zhao and Jha (2021); Meguerdijian and Jha (2021). To resolve this partially, I imposed a 52 Figure 3.2: The modified Bandis model to model nonlinear fracture mechanics in response to dy- namically evolving normal and shear stresses Li et al. (2017). The top left quadrant shows the behavior of a fracture in tension with K ni denoting the initial normal fracture stiffness. The bot- tom right quadrant shows a fracture in compression, displaying the initial mechanical aperture a 0 under zero excess effective normal stress, the maximum aperture of the fracture under compression (V m ), and the change in mechanical aperture ∆ a n due to normal stress perturbation. Hysteresis within a loading-unloading cycle on the fracture is shown with the triangle marker (loading or withdrawal period) and the circle marker (unloading or injection period) lines. The upper right quadrant shows the dependence of shear stress (τ ) on shear displacement (δ h ). The mobilization of JRC is divided into the pre-peak and post-peak periods, separated by point B corresponding to peak shear displacement δ peak and peak shear strength τ peak . As the shear displacement increases to the residual shear displacement δ res , the residual shear strength τ res is slowly reached at points C and D (upper right quadrant). limit on the maximum number of times a fracture segment is capable of slipping, and I set this limit to be equal to the number of sub-segments this fracture is divided into when intersecting the matrix mesh cells. Figure 3.2 illustrates my fracture mechanics model. 3.4.6. Initial and Boundary Conditions The coupled flow, transport, and geomechanics problem requires suitable boundary conditions for thethreesub-problems. ReferringtothephysicalmodelinFigure3.1,Iapplythefollowingboundary conditions: zero normal displacement on the bottom and right boundaries, compressions on the left and top boundaries with a prescribed stress anisotropy due to regional tectonics, zero normal fluid flux on the left and right boundaries, fixed pressures on the top and bottom boundaries which 53 may represent constant pressure of producing wells at those locations, and the natural outflux boundary condition for the tracer transport problem. The cyclic well in the domain is controlled by a prescribed injection rate q inj until the tracer breaks through at the boundary. This is followed by a soaking period of 30 days when only the in-situ fluid is injected. Finally, the tracer withdrawal episode is modeled by a prescribed well withdrawal rate q prod for 90 days. Units of q inj and q prod are volumetric rates per unit length in the third dimension. Initializationofacoupledflow-transport-geomechanicsproblem, especiallyinlowpermeabilityrocks with highly conductive fractures, is a challenging task. This is also true during lab experiments involving such rocks, e.g. initialization of pressure and concentration fields in a shale core subjected to coreflooding can take days. The goal of initialization is to achieve mechanical and hydrostatic equilibria such that the initial strains and fluid fluxes inside the domain are zero in absence of sources/sinks (i.e. wells). Spatially uniform fields of pressure, concentration, stresses, and dis- placements, e.g. zero excess pressure and displacements combined with initial stresses equal to boundary tractions, are prescribed as initial conditions that are close to, but not exactly at, the equilibrium. During the initialization period of the simulation, I solve the governing equations with the prescribed boundary compressions, initial conditions, body forces (in case of a 3D model), and initial rock-fluid properties to achieve the state of poromechanical equilibrium. No wells or tracers are active during the initialization period. The state variables, e.g. pressure and stress, at the end of initialization become my reference state. Stress-dependent properties such as porosity and permeability fields also experience perturbations during initialization and change from their initial user-prescribed values. I define ‘excess’ quantities, e.g. excess pressure, at a time t as the difference between the value at t and the value at the reference time. The cyclic well starts operating at the end of the initialization period with a fixed tracer concentration prescribed at the well location. This completes the definition of the coupled problem. The governing equations to be solved are equations (3.1), (3.2), (3.13) (after substitution to express it in terms of displacements), (3.14), and (3.15) subject to the initial and boundary conditions mentioned above. The constitutive equations for the fluid compressibility, tracer viscosity, stress-strain relation, well connectivity, mass transfer 54 fluxes, and fracture mechanics provide the necessary closure of the system of equations. It is a non- linear system of coupled differential equations with multiple length and time scales that correspond to multiple physical mechanisms in the system. As mentioned earlier, I select pore fluid pressures p m and p f , displacement vector u, and concentrations c m and c f as the primary unknowns of the problem. The effective stress σ ′ , porosity ϕ , Darcy velocities v m and v f , fracture aperture a h , and permeability k f are the secondary variables which I update explicitly by post-processing the primary variables. 3.5. Numerical Model Solvingthecoupledproblemdefinedaboverequiresanumericaldiscretizationschemethatcanhonor the continuityrequirements of the primary and secondaryvariables while allowing for multiple scales in the discrete system. Computational efficiency and robustness are the other two requirements on the numerical model because I intend to apply the model to aquifers with geometrically complex fracture networks, such as the one shown in Figure 3.1 with more than 20 fractures of arbitrary lengths and orientations intersecting each other. To satisfy these requirements, I use the finite vol- ume method to discretize the flow and transport equations, and the Galerkin finite element method to discretize the mechanics equations Jha and Juanes (2014); Tran and Jha (2020). I discretize the domain using an n d -dimensional uniform Cartesian mesh (n d = 2 for 2D case). The fractures are discretized uniformly on a separate n d − 1 dimensional mesh. The weak forms of the flow and mechanics equations are obtained by following the standard variational recipe: multiply the strong form by an appropriately defined smooth test function, integrate the equation over the domain, apply integration by parts and reduce the order of differentiation on the primary unknown, apply the divergence theorem, and impose Neumann boundary conditions Tran and Jha (2020). I use cell- centered piecewise constant pressure and concentration degrees of freedom and bilinear functions with nodal displacement degrees of freedoms. This choice satisfies the inf-sup condition because I consider a compressible system with finite values of the solid grain compressibility 1/K s and the fluid compressibility C f , which are typical of most subsurface formations and deep aquifers. The Darcy flux term in the flow equation is linearized using the two-point flux approximation. The advective flux term in the ADE is discretized using upwinding and the central finite difference is 55 Transport v ( p, k f ) μ f ( c ) k f ( δ h , σ n ) Figure 3.3: Flow and geomechanics are coupled two-way through strain-dependent porosity ϕ (ϵ v ) in the fluid mass balance equation (Equation 3.1) and the fluid pressure dependent effective stress σ ′ (p) in the linear momentum balance equation (Equation 3.13). The flow and transport are coupled two-way through the concentration-dependent fluid viscosity µ f (c) and pressure-dependent velocity field v(p). The geomechanics is coupled one-way with fracture mechanics by providing the nodal displacement, stress, and strain tensors to the Bandis model for slip activation, slip magnitude calculation, and fracture aperture update. The fracture mechanics delivers updated stress-dependent permeability k f (δ h ,σ ′ n ) to the flow solver based on shear slip δ h and normal dilation, which is dependent on the effective normal stress σ ′ n . used to discretize the diffusive flux term. We use the first-order accurate Backward Euler method to implicitly integrate the flow equations in time. Finally, I obtain the linear system of equations to be solved at each time step. 3.5.1. Coupling strategy and solution scheme Our coupled computational framework consists of four parts: geomechanics, flow, transport, and fracture mechanics (Figure 3.3). Each component is connected to at least one other component via a one-way or two-way coupling link. The strength of coupling is dictated by different param- eters: α , E, C f for flow-geomechanics coupling; R, k f for flow-transport coupling; and JRC, JCS for fracture mechanics-geomechanics coupling. We utilize a fully-coupled simultaneous solution approach to solve for displacement and pressure fields at each time step Jha and Juanes (2007); Lewis and Sukirman (1993b); Tran and Jha (2020). The matrix-vector system of linear equations 56 becomes           K − Q 0 0 Q T S m +δt A m − δt T mf − δt T mw 0 δt T fm S f − δt A f δt T fw 0 T wm T wf − A w                     U m p m p f p w           n+1 =           f u δt f p,m +Q T U m,n +S m p m,n δt f p,f − S f p f,n q w           (3.21) where U is the vector of nodal displacements from the finite element interpolation u =N T U, N is the matrix of bilinear shape functions; p m , p f , and p w are the vectors of cell-centered matrix pressure, fracture segment pressure, and well segment pressure unknowns respectively, n is the previous time step, and n+1 is the current time step at which displacements and pressures are solved such that δt = t n+1 − t n . The global stiffness matrix is assembled from element stiffness matrices, K i = Z Ω i B T i D i B i dΩ (3.22) whereB is the strain-displacement matrix computed by applying the symmetric gradient operator to N, and i is the element or cell index. There are two kinds of transmissibility matrices in Eq. (3.21). The global transmissibility matrix A is assembled from interfaces between cells of a single medium, indicated by the superscript m, f, or w in the equation. The global transmissibility matrixT is assembled from the individual interface transmissibilities between two different media, which are indicated by the two superscripts mf, wf, and mw in the equation. See my earlier paper Tran and Jha (2020) for expressions of the transmissibility matrices. The global storativity or compressibility matrixS for either the matrix or the fracture is a diagonal matrix with components, S ii = R Ω i 1/M i dΩ . The flow-geomechanics coupling matrix at the element level is, Q i = Z Ω i α B T i I d dΩ . (3.23) 57 The right-hand side of the weak form of the mechanics problem f u consists of traction (and body forces in case of a 3D model), f u = Z Γ t N T σ dΓ+ Z Ω N T ρ g dΩ (3.24) where Γ t denotes Neumann boundaries of the domain with the prescribed tractions σ = [0,σ H ] T on the top boundary and σ = [σ h ,0] T on the left boundary. The right-hand side of the weak form of the flow problem f p consists of prescribed fluxes and well terms. I solve Equation (3.21) using a direct solver. The Darcy velocity field is computed from the pressure solution. Next, the ADE is solved for the unknown concentrations. The flow-geomechanics problem is solved again at the same time step with the new concentration and fluid mobility fields. This process is repeated iteratively within an inner loop inside the time loop until pressure, concentration, and displacement solutions converge at the specified time step. The coupling between the poroelastic solver and the transport solver is implemented using a sequential one-way coupled approach. For the range of R considered in this study, the transport-to-deformation coupling strength is expected to be weak because there is no direct effect of concentration on the mechanics equation (3.13), and the effect of concentration on pressure is limited to fluid mobility. For the range of parameter values considered, the inner loop takes less than seven iterations to achieve convergence at any time step during the course of the simulation. 3.6. Model Verification and Benchmarking In recent literature, EDFM has been successfully validated against traditional methods of fracture modeling, e.g. EDFM’s pressure solutions are compared to pressure solutions from discrete frac- ture network simulations of multiphase flow in complex fracture networks Hajibeygi et al. (2011); Shakiba (2014). Further validation studies have been conducted using both a fine-grid simulation and a semi-analytical solution Shakiba (2014); Jansen et al. (2018). We conducted two model verifi- cation studies, shown in the supporting information Figures S1 and S2, where I compare my EDFM pressure solution with a commercial multiphysics simulator solution and an analytical solution, respectively. Below I present results from benchmarking on the Mandel consolidation problem to 58 verify poroelastic coupling in my framework. 3.6.1. Mandel-Cryer Effect Totesttheaccuracyofporoelasticcouplinginmycomputationalframework,Iappliedtheframework to the classic Mandel problem. An elastic, isotropic, homogeneous rock slab that is infinitely long in the z-direction is constrained between two frictionless, impermeable rigid plates (Figure 3.4). Due to the symmetry of the problem, the computational domain can be limited to the upper right quadrant (shaded gray in the panel (a) of Figure 3.4) with size L x × L y . The sudden application of compression at the boundary results in a positive excess pressure and linear displacements in the domain, also known as the undrained deformation, because the pore fluid does not have sufficient time to drain yet. The vertical left and right boundaries at x=± L x (or dimensionless x-coordinate x d =± 1) are free to deform at ambient pressure, which allows the fluid to drain over time. This implies that x = 0 can be modeled as a roller boundary condition. The analytical expressions of pressure, displacement, stressaregivenintheliteratureAbousleiman et al.(1996);Castelletto et al. (2015b). Figure 3.4(b) shows the dimensionless pressure profile along the horizontal distance in the shaded domain at four dimensionless time steps. The dimensionless pressure p d is obtained by normalizing the pressure by p 0 =FB(1+ν u )/(3L x ) where B =ν/ (ν 2 +(K dr /M)) is the Skempton coefficient and ν u =(3ν +Bν (1− 2ν ))/(3− Bν (1− 2ν )) is the undrained Poisson’s ratio. At early times (t d <0.1), the x=L x boundary becomes softer due to pressure dropping from fluid drainage at that boundary. The stress is transferred inwards into the domain, which is stiffer. This causes the pressure near thex=0 boundary to increase above the undrained pressure value, i.e. p d >1, which is known as the Mandel-Cryer effect. Next, I discuss results from a few representative numerical simulations. 3.7. Results All simulation cases below follow the same operational sequence for the cyclic well: the well starts injecting tracer until the tracer breaks through at either the top or bottom boundary, then the well injects the in-situ fluid for the 30-day soaking period, and finally, the well withdraws fluids for the 120-day withdrawal period. The withdrawal rate is set at 1.5 times the injection rate, which is 59 x Y b a 2F 2F x/L x L x L y p/p 0 0.0 0.2 0.4 0.8 1.0 0.6 0.0 0.2 0.4 0.8 1.0 0.6 1.2 (a) (b) Figure 3.4: (a) The Mandel problem setup with the origin of the x-y coordinate system positioned at the center of a plane strain 2D slab of dimension 2L x × 2L y . Symmetry allows modeling of a quarter of the slab, producing a computational domain of size L x × L y . The system is under a compressive vertical force denoted as 2F applied instantaneously at t = 0. (b) The dimensionless pressurep/p 0 vs. the dimensionless x-coordinatex d =x/L x at different dimensionless times t d . The lines show the simulation results and the dots are the analytical solutions. The simulation agrees with the analytical result and captures the Mandel-Cryer effect at the early time step of t d =0.05. typical in a push-and-pull test. 3.7.1. Base Case We define a base case using the values of the mesh, initial condition, and boundary condition pa- rameters shown in Table 3.2, together with the material parameter values shown in Table 3.3. Some of these values are taken from a published tracer flowback case study Li et al. (2016). Utilizing the proposed computational framework, I analyze the base case results to identify the key characteristics of flow, transport, and geomechanical behavior. Transport response in a coupled flow-geomechanics framework The tracer distribution in the matrix (colored contour lines) and fracture (colored straight lines) are monitored at three selected time steps: end of injection which is also the tracer breakthrough time, end of soaking, and end of 30-day withdrawal (Figure 3.5). During injection, the effect of the fracture network becomes evident as the tracer advects along the fractures that intersect the wellbore. Some amount of tracer directly enters the matrix due to diffusion. Tracer diffusion is also observed across the fracture-matrix interfaces, especially the ones situated within the tracer fingers that protrude outward from near-well fractures. Among the fractures that intersect the well, there 60 Parameter Symbol Value Number of cells in the x-direction N x 100 Number of cells in the y-direction N y 100 Number of cells in the z-direction N z 1 Domain length in the x-direction L x 200 meter Domain length in the y-direction L y 200 meter Domain length in the z-direction L z 75 meter Top boundary compression S H 82 MPa Left boundary compression S h 50 MPa Initial excess pressure ∆ p 0 0 MPa Initial pressure p 0 35 MPa Prescribed top and bottom boundary pressure p out 35 MPa Well injection rate q inj 7.4× 10 − 5 m 3 /sec/m Well production rate q prod 11.1× 10 − 5 m 3 /sec/m Poisson’s ratio ν 0.27 Solid grain density ρ s 2650 kg/m 3 Drained bulk modulus (matrix and fracture) K dr 21.9 GPa In-situ fluid viscosity µ 0 1.2× 10 − 3 Pa-sec Injected tracer viscosity µ 1 0.16× 10 − 3 Pa-sec Fluid density ρ f 1000 kg/m 3 Matrix permeability k m 5 millidarcy Initial fracture permeability k f 0 1256 darcy Initial matrix porosity ϕ m 0 0.1 Fracture porosity ϕ f 0.1 Unconfined compressive strength UCS 150 MPa Table 3.2: Simulation parameters Parameter Symbol Value Drained Young ′ s modulus E 25.55 GPa Biot modulus M 16.9 GPa Diffusion coefficient D d 2× 10 − 7 m 2 /sec Fluid compressibility C f 0.5× 10 − 9 Pa − 1 Biot coefficient α 0.5 Log viscosity ratio R 2 Joint compressive strength JCS 130 MPa Joint roughness coefficient JRC 8 Table 3.3: Base Case parameters 61 existsapreferentialdirectionforthemovementofthetracerbasedontheorientationofthefractures. Due to the coupling between the tracer viscosity and the pressure field, this preferential transport controlshowmuchthefracturepermeabilitychangeswiththein-situstressfield. Themorefavorably aligned a fracture is to the regional principal stress directions (in this case, fractures oriented quasi- parallel to S H , i.e. orthogonal to S h ), the more penetrating the fracture’s tracer distribution is. In fact, the tracer breaks through at the bottom boundary by flowing through the most vertically oriented fracture in the bottom right of the domain (panel (a) of Figure 3.5). During the soaking period, since there is no additional tracer injected, the tracer concentration field displays a more spread-out distribution in the aquifer with smaller tracer concentration values around the well. The main transport mechanism is diffusion across matrix-fracture and matrix-matrix interfaces, which contributes to less pronounced tracer fingers. The pressure field is continuous but not necessarily smooth across a fracture, demonstrated by the kinks in the pressure contours of bottom panels in Figure 3.5. During the tracer withdrawal period, the fracture network again shows its impact on the tracer distribution as displayed in the relatively granular tracer contours, despite a lower value of the concentration gradient inside the fracture compared to the injection period. Thetracerretrievaldataisavailableinmosttracersurveystudies, whichcanyieldusefulinformation if used correctly. The tracer retrieval curve can be used to deduce fracture morphology Li et al. (2016), especially if the fracture network is complex with non-uniform density, orientation, and permeability. The temporal characteristics of tracer retrieval curves at various locations along the horizontal well are shown in Figure 3.6. During the early period of injection (t < 5 day), concentration around the well rises sharply to reach a peak level of c w ≈ 0.5, and the fractures intersecting the well are saturated with the tracer. Then the tracer level increases at a much slower rate and the tracer slowly diffuses into the matrix through matrix-matrix or fracture-matrix interfaces. As expected, during the soaking period when the well does not inject any tracer, the tracer concentration level at the well drops sharply to zero. The soaking period yields minimal differences between different well locations due to a lack of tracer injection into the domain during this period. During withdrawal, a more gradual build-up is observed as the tracer returns to the well. Transport is non-uniform along the length of the well due to a non-uniform fracture density 62 Time 0.8 0.6 0.4 0.2 0 (a) push stage (c) pull stage breakthrough 20 m 0 4 -12 -9 -6 -18 2 6 12 8 10 -3 0 -15 MPa MPa Figure 3.5: Top row: Tracer concentrations in fractures and in the matrix at three successive time steps of the base case model: (a) at tracer breakthrough (t BT ) which is also the end of injection period, (b) at the end of the chase fluid injection period (soaking), and (c) 30 days into the withdrawal period. The well is denoted by the blue horizontal line. Matrix concentration c m is shown via contours drawn at an increment of 0.1. Fracture concentration c f is shown via colored straight lines drawn along the fracture length. Both follow the same color scale shown on right. The shape and size of the contours are controlled by the connectivity of the well to the fractures and the transport processes. Bottom row: corresponding excess pore pressure distribution at the three selected time steps. Color scales of the bottom panels are different, implying the sign reversal of excess pore pressure during push and pull periods. 63 (number of fractures intersecting the well per unit length of the well). This is true in my study as indicated by Figure 3.1: the left and right thirds of the wellbore have fewer fractures compared to the middle third of the well. Since fractures contribute significantly to the advective transport, which is away from the well during injection and toward the well during withdrawal, the middle point curve shows a higher peak concentration during injection and a lower peak concentration during withdrawal compared to the curves at the endpoints (Figure 3.6). An interesting observation is an increase in separation between the three curves during the with- drawal period compared to injection. This points to hysteresis in advective transport through a medium with discrete fractures. Here, hysteresis refers to the observation of different concentra- tions at a point as the direction of flow reverses, e.g. due to switching the well from injection to withdrawal. Hysteresis in diffusion-dominated transport is well-understood because the effect of diffusion on the concentration gradient at a point is monotonic with time, even if the flow direc- tion reverses. However, hysteresis in advection-dominated laminar flow (high Peclet number and low Reynolds number) is a topic of research and has many applications in the mixing industry, e.g. in microfluidics. my study set-up allows us to observe, quantify and control the magnitude of hysteresis in tracer transport through a fractured aquifer. Separation of the three concentration profile curves, which correspond to regions with three different fracture densities, also suggests that these observables have a diagnostic power. If such c vs. t data is acquired in the field by zonal fluid sampling along the wellbore at multiple time steps, the curves can be used to quantify the fracture density in the aquifer away from the well. This is valuable because electrical resistivity- based fracture logging tools sense only 1-6 inches away from the wellbore. The fracture density in the left and right segments of the well are similar, yet there exists a slight difference between the two tracer retrieval curves. This is due to the difference in the fracture orientations of these two segments: the orientations are more orthogonal to S H in the left segment than in the right segment (Figure 3.6). The variation in fracture orientation determines the action of the stress tensor on these discontinuities, resulting in different permeability changes during the simulation and hence a variation in the tracer retrieval curves. The impact of fracture permeability and orientation on the tracer retrieval curve will be revisited in later sections as they are strongly related to the stress 64 t (days) Left Endpoint Right Endpoint Midpoint End of Injection Start of withdrawal (a) (b) 50 m Figure 3.6: Evolution of the tracer concentration at selected points along the well. The curves at the right and left ends of the well are similar during the injection period while the middle point shows a higher concentration due to a higher fracture density near the center. The soaking period does not show any distinction among the three curves because the tracer is neither injected nor withdrawn. The withdrawal period curves are more separated from each other compared to the injection period because both advection and diffusion processes feel the fracture morphology deeper in the aquifer during later periods. field. Stress response to pressure coupling Next, I analyze how induced stresses in the aquifer vary with time and space during the episodic transport. The isotropic part of the deformation (volumetric expansion or contraction) is quantified by ϵ m v which is controlled by the effective volumetric stress σ ′,m v , which is related to the total volumetric stress σ m v = trace(σ m )/n d as σ ′,m v =σ m v +α m p m . The deviatoric part of deformation is controlled by the von Mises stress σ m VM = √ 3J 2 , where J 2 is the second invariant of the deviatoric stress tensor (σ m − σ m v I d ) with I d is the identity tensor Tran and Jha (2020). To characterize the state of stress in the domain, three different stress quantities are chosen: effective volumetric stress σ ′,m v for isotropic deformation, von Mises stress σ m VM , and poroelastic shear stress invariant I ′ MC for shear-induced rock failure and plasticity. With the effective stress σ ′,m denoted in Voigt engineering notation as [σ ′ xx ,σ ′ yy ,σ ′ xy ] T , the von Mises stress for the plane stress assumption are defined as follows: σ ′ VM = (σ ′ xx ) 2 +(σ ′ yy ) 2 − σ ′ xx σ ′ yy +3(σ ′ xy ) 2 1/2 (3.25) 65 The excess poroelastic shear stress I ′ MC is calculated as a function of the first and second stress invariants I ′ 1 and p J ′ 2 , representing the invariant form of the Mohr Coulomb (MC) yield func- tion Borja (2013). I ′ 1 =2σ ′,m v (3.26) J ′ 2 = 1 6 [(σ ′ xx − σ ′ yy ) 2 +σ ′2 xx +σ ′2 yy ]+σ ′2 xy (3.27) I ′ MC = q J ′ 2 − 1 3 I ′ 1 sin(arctan(µ s )), (3.28) whereµ s isthestaticfrictionalcoefficientofthehostrock. Thetemporalevolutionofexcesspressure and stress quantities are shown in Figure 3.7 at three locations: near the midpoint of the well (left diagram), the center of the upper half of the domain (middle diagram), and near the top boundary (right diagram). At the very beginning of the injection period, pressure rises due to the introduction of tracer fluid in a low permeability reservoir and the time not being long enough for fluids draining off in the top and bottom boundaries. The pressure then stabilizes and starts decreasing slightly because the injected tracer has a lower viscosity than the in-situ fluid it displaces in the domain. The pressure plateaus during much of the soaking period as no additional tracer is introduced. Then during the withdrawal period, the excess pressure turns negative, meaning the pressure is lower than the initial pressure as intended to allow draining of fluids in the cyclic well from the domain. The negative excess pore pressure continues to decrease throughout the withdrawal period as more tracer is displaced from the domain and replaced by the viscous resident fluid encroaching from the boundaries. Quantities characterizing the stress state evolve synergistically with the excess pore pressure in the matrix (Figure 3.7). Temporal evolution of the volumetric, deviatoric (von Mises), and Mohr- Coulomb first invariant stresses follow the excess pressure trend closely, confirming that the pressure field is driving the coupled dynamics of the system. During the injection period, pore pressure is higher than the initial pressure (positive excess pore pressure or ∆ p > 0), leading to a decrease of the effective volumetric stress, which means the host rock becomes more tensile. Because of the sign convention of the compressive state of the reservoir (σ ′,m v <0), this means that the excess effective 66 volumetric stress is positive (∆ σ ′,m v > 0) as shown in the lines with the diamond marker of Fig- ure 3.7. A similar argument can be made during the withdrawal period. As the rock becomes more volumetrically compact during production, the ∆ σ ′,m v behavior mirrors the trend of ∆ p. Investigat- ing the rich and complex spatial and temporal evolution of σ ′ v is important in practical applications. For example, in designing a new infill hydraulically fractured well, an optimal strategy is to choose drilling targets in a low compressive stress environment to facilitate the opening of stimulated frac- tures Guo et al. (2019). Besides the reversal of excess volumetric stress and excess pore pressure between injection and withdrawal periods, excess Von Mises stress and excess poroelastic shear stress are positive throughout these two episodes due to their definitions in Equation 3.28. Near the injection well, the Von Mises stress is very close to the excess effective volumetric stress because the excess shear stress is small compared to the excess effective normal stress ( ∆ σ ′ xy ≪ ∆ σ ′ xx ) and stress state is almost isotropic (∆ σ ′ xx ≈ ∆ σ ′ yy ) (panels (a) and (b) of Figure 3.7). Under these conditions, Equation 3.25 becomes σ ′ VM ≈ [(σ ′ xx ) 2 +3(σ ′ xy ) 2 ] 1/2 ≈ σ ′ xx . However, as I move away from the center of the domain, where the excess pressure perturbation is minimal, as the relative magnitude of ∆ σ ′ xy grows in comparison to ∆ σ ′ xx , ∆ σ ′ VM dominates ∆ σ ′ v . This suggests that shear-induced deformation becomes dominant away from the well as shown in the separation of the curves denoted in triangle and diamond markers. This can be explained by less tension (quantified by I ′ 1 ), or the increasing prominence of remote shear components (quantified by p J ′ 2 ), or both factors away from the injection well. The injection-induced rock expansion around the injection well is accompanied by slight compression at the boundary. Thedirectionalre-alignment ofmaximum andminimum horizontalprincipal stresses dueto pressure perturbation is of importance for the optimal design of a hydraulic fracturing job and well trajectory lay-out. When an existing well is operated cyclically as an injector and producer, the operator has to estimate the change in the principal stress direction, so a future drilling well can be placed and oriented optimally in the field. For example, in ultra-low permeability shale or sandstone formations that require hydraulic fracturing to boost wellbore injectivity/productivity, the objective is to align the well along the minimum principal stress direction, so that the fractures orthogonal to the wellbore can open up easily against the minimum compression. In these scenarios, a rigid approach 67 t (days) t (days) (MPa) ∆p ∆ σ t (days) ∆p ∆ σ or ∆I ′ MC ∆ σ ′ VM ∆p ∆ σ ′ V (a) (b) (c) Figure 3.7: The temporal evolution of excess pore pressure (∆ p) in the matrix, excess effective volumetric stress (∆ σ ′ v ), and excess Mohr-Coulomb stress Invariant (∆ I ′ MC ) at different locations in the domain, which are marked in the inset figures in the bottom left. The locations are (a) near the midpoint of the well, (b) center of the upper half of the domain, and (c) near the top boundary. Pressure (star marker), Mohr-Coulomb stress invariant (black marker), Von Mises stress (blue marker), and effective volumetric stress (diamond marker) evolve with time as the well goes through tracer injection, soaking (or chase fluid injection), and withdrawal episodes. to drilling wells without any consideration to the dynamically evolving stress field is guaranteed to yield unfavorable stress-alignment conditions and sub-optimal performance of the wells. Figure 3.8 shows the rotation angle ∆ β S H of the maximum principal stress S H , describing a re-orientation of the local principal stress field. I observe regions on the two sides of the well with opposite changes in the stress orientation (red and blue shades), and I observe larger changes of ∆ β S H near the left boundary compared to those near the right boundary (lighter shades). This can be explained by the contribution of anisotropic fracture network orientation (predominantly oriented Northwest- Southeast) to the non-uniform movement of tracer, perturbing the pressure field correspondingly. MoretracerinvasionintheNorthwest(NW)andSoutheast(SE)cornersleadstoaporefluidmixture of lower viscosity and less pressure drop away from the injection well (Figure 3.5). This results in the reorientation of S H in the NW-SE regions to be more aligned with the predominant fracture orientation, coinciding with the direction of the slowest change of pore fluid pressure. In addition, fluid-solid poroelastic coupling results in the occurrence of small shear stress perturbation (which is also observed in Figures 3.16 and 3.17) during injection, further augmenting the reorientation of principal stresses locally. Flow properties response to fracture dynamics and stress coupling Basedontheprescribedinitialstressfieldandfractureorientations, aninitialdistributionoffracture permeabilities is calculated using empirical correlations from 141985Barton et al.Barton, Bandis, 68 1.5 -1 -1.5 1 0 0.5 -0.5 ∆ β S H ( o ) (a) (b) 50 m Figure 3.8: Change in direction, β S H (degrees), of the maximum horizontal stress S H , measured from its initial direction along the y-axis.A positive angle means the principal stress rotates counter- clockwise while a negative angle indicates clockwise rotation. The dashed lines indicate the new S H direction. Green thick lines are the embedded fracture network. The thick blue line is the well trajectory. Two time steps are shown: (a) at the end of injection, and (b) at the end of withdrawal. and Bakhtar (). In the initialization step, fracture orientation contributes to the fracture perme- ability heterogeneity observed in the line with the square marker in Figure 3.9. Depending on how favorable the fracture orientation is compared to the anisotropic stress field applied initially, each fracture segment responds differently to the stress field. The empirical cumulative distribu- tion function (ECDF) of the fracture permeability is shown in Figure 3.9 at selected time steps. In my method, each fracture is composed of multiple fracture sub-segments that undergo different stress states based on their respective surrounding nodal stress state in the matrix, hence generating assorted permeability changes for each segment. When injection starts, besides aperture opening due to decreasing effective compression, shear-induced dilation facilitates the additional opening of fracture aperture Asadollahi and Tonon (2010). Therefore, the rich complexity of the stress field explains the spatial heterogeneity of fracture permeability despite the tight range of reference start- ing values (Figure 3.9). Because of the reduction in acting normal effective stress and the effect of shear-induced dilation, fracture permeability in the injection period (star marker) is approximately 17% higher compared to the withdrawal period (diamond marker). During injection, shear stress accumulates until shear slip is activated on a particular fracture sub-segment based on the Mohr- Coulomb failure criterion (i.e. fracture shear stress exceeds the product of fracture effective normal stress and static frictional coefficient τ f ≥ µ s σ ′,f n ). The shear-induced dilation effect is apparent in 69 ECDF 1320 1220 1270 1050 1090 1070 k f (Darcy) Figure 3.9: Empirical cumulative distribution function (ECDF) of the fracture permeability k f after initialization (square marker), at the end of injection (star marker), and at the end of with- drawal (diamond marker). Fracture permeabilities increase during injection and decrease during withdrawal. The two inset figures show the spatial distributions of k f at the respective timesteps: injection in the right inset and withdrawal in the left inset. The insets use different color scales to highlight the different fractures activated during different episodes–injection and withdrawal. augmenting the high-end values of k f during injection. Therefore, the spread of permeability range is amplified in the injection period compared to the withdrawal period. The temporal evolution of fracture permeability is governed by both fracture orientation and excess effective normal stress acting on the fracture surfaces (Figure 3.10). On panel (a) of Figure 3.10, three representative fractures are chosen from the segments that have physically slipped during the courseofinjection. Allfracturesstartwithaconstantandprescribedvalueofk f thatchangesduring the model initialization step discussed above, appearing as the first drop in k f (Figure 3.10(a)). Then, during the injection period, when no shear slip occurs, k f increases very slowly because the effective normal compression decreases slowly. The permeability trend exhibits an abrupt jump due to a shear slip event with the timing and magnitude of the events depending on the fracture orientation and the stress state. Fractures that are more favorably aligned (or quasi-parallel to the vertical S H , i.e. smaller θ f ) experience a larger increase in their stress-dependent permeability because their apertures open up against the minimum horizontal principal stress S h . The larger the angle, the smaller the fracture opening is due to increasing effective normal stress and smaller associated shear-induced dilation. As a result, the evolution of k f can be split into two periods 70 Slipped Unslipped ∆ (a) (b) k f (Darcy) k f (Darcy) (MPa) Figure 3.10: The temporal evolution of fracture permeability and its dependence on fracture orien- tation θ f . The fracture orientation θ f (degrees) is defined as the angle between the fracture plane and the initial S H direction (i.e. true North). It ranges from − 90 ◦ (clockwise from true North) to 90 ◦ (counter-clockwise from true North). The true North is the direction of S H in Figure 3.1. Panel (a) plots the permeability of three representative fracture segments that slipped and experienced shear-induced dilation during the injection period. The legend shows different orientations of each fracture segments. The horizontal time axis is plotted on log scale to illustrate the early injection period better. The mini-panel shows the locations of the fracture segments that are plotted. Panel (b) plots fracture permeability and orientation for all fracture segments where circle markers repre- sent unslipped segments and square markers denote slipped segments. The points are color-coded by the absolute magnitude of excess normal effective stress acting on the fracture segments ∆ σ ′,f n . with the corresponding mechanisms: the inter-slip period when the effective normal stress dictates k f and the co-slip period when the shear slip dictates k f . Interestingly, these shear slip events can be physically connected to microseismic episodes, if additional related studies are pursued. Finally, a sharp drop of fracture permeability between injection and withdrawal is expected owing to the sign reversal of the excess pore pressure. To increase the confidence of my fracture dynamics model, the quantitative trends of fracture permeability versus stress magnitude, stress history, and slip displacement are checked against existing literature Lamur et al. (2017); Watanabe et al. (2008). As observed in panel (b) of Figure 3.10, only a narrow range of fracture orientation facilitates shear slip (between 25 ◦ and 50 ◦ ). Slipped segments show higher permeability change compared to unslipped segments mainly due to the shear-induced dilation effect. Among the slipped segments, permeability change is highest for fractures of lower orientation angles (≈ 25− 30 ◦ ) and higher values of the excess effective normal stress ∆ σ ′,f n (denoted as red filled squares in Figure 3.10(b)). 71 Parameter Low· Base· High Unit Biot coefficient α 0.2· 0.5· 0.9 Fluid compressibility C f 0.4· 0.5· 1 GPa − 1 Log viscosity ratio R 0· 2· 2.5 Young Modulus E 10· 25.55· 50 MPa Joint Compressive Strength JCS 120· 130· 150 MPa Joint Roughness Coefficient JRC 8· 8· 20 Diffusion coefficient D d 2× 10 − 8 · 2× 10 − 7 · 2× 10 − 5 m 2 s − 1 Table 3.4: Sensitivity parameter values. 3.8. Sensitivity Analysis To better understand the coupling mechanisms between transport and deformation processes in an aquifer, whose properties can vary over a wide and/or uncertain range, I perform a sensitivity analysis to investigate the influence of key rock-fluid properties on different outputs of the Base Case defined above. As the outputs of the coupled simulation can be classified into flow (e.g. pressure), transport (e.g. degree of mixing), and geomechanics (e.g. stresses) categories, so the rock-fluid properties can also be classified into the corresponding categories that control the flow-geomechanics coupling strength (Biot coefficient α , fluid compressibility C f , rock’s Young modulus E), flow- transport coupling strength (log viscosity ratio R, diffusion coefficient D d ), and geomechanics- fracture mechanics coupling strength (JRC, JCS). Table 3.4 lists the range of investigation of these sensitivity parameters. The base case values are in the middle of the sensitivity range except for JRC, because JRC is empirically constrained by the original Bandis model based on the magnitude of in-situ stresses, which also narrows the sensitivity range of JCS. 3.8.1. Poroelasticity and fracture dynamics affect mixing and spreading Here, I investigate how different types of coupling lead to changes in the solute transport be- havior. Spreading and mixing are critical transport phenomena, which can be quantified using the tracer breakthrough time t BT and the domain-averaged degree of mixing at breakthrough, ζ BT Tran and Jha (2020); Jha et al. (2011a,0,0). The degree of mixing ζ can be defined in terms of the matrix concentration field: ζ (t)=4(c m ϕ d ϕ d (1− c m )), (3.29) 72 where ϕ d = ϕ/ϕ 0 and the overbar indicates volume averaging. ζ ranges from 0 to 1, with ζ = 0 representing an unmixed state and ζ = 1 representing a completely mixed state between injectant and resident fluids. The tornado plot in Figure 3.11 shows the range of magnitude change (in percentage) of t BT and ζ BT compared to the base case for all sensitivity parameters. I observe different levels of impact on t BT andζ BT for parameters belonging to different types of the coupling mechanism. For all cases except in one scenario, I observe a monotonous relation between each sensitivity parameter listed in Table 3.4 and the macroscopic transport metrics t BT andζ BT . There exists a non-monotonous trend in the relationship between Biot coefficient α and t BT , which means tracer breakthrough is fastest in a moderately consolidated reservoir. This nonlinearity may partly stem from the quadratic relationship between α and incremental porosity change δϕ m expressed in equation 2.8. For an incompressible system, or in the drained deformation case, the Biot coefficient can be interpreted as the change in porosity for a unit change in the volumetric strain Tran and Jha (2020). This phenomenon is also likely because α is present in both the fluid mass balance equation and the mechanics linear momentum balance equation, which complicates its effect on the pressure fieldandresultsinanonlinearrelationshipwiththetransportmetrics. Fortherangeof α Iexamined (from 0.2 to 0.9), changes in t BT and ζ BT are apparently small (approximately 6%) but could be significant for a real site where the absolute breakthrough time is in the order of months, meaning a 6% error margin could translate into an error in the estimated t BT of the order of a week. In those cases, it will be costly to ignore the coupling in the model. Besides, here, the degree of mixing is taken to be the average value for the whole domain, which may provide only low resolution information in a fractured aquifer where tracer generally accumulates around the major fractures. I will have a better look at the mechanical coupling impact on the spatial distribution of the tracer in the next section. Regarding the fluid compressibility, higher C f (more compressible fluid, weaker flow-geomechanical coupling) leads to a faster breakthrough (Figure 3.11). Higher C f leads to lower Biot modulus and higher values of the pressure time derivative, i.e. the fluid accumulation term. The Darcy velocity increases, especially in high permeability channels, which effectively reduces t BT . Regarding rock’s modulus, lowerE (softer rock, stronger flow-geomechanical coupling) leads to a faster breakthrough 73 of the tracer. Compared to flow-geomechanics coupling parameters ( C f and E), fracture dynamics parameters (JRC, JCS) have a larger impact on transport: tracer breakthrough is retarded (t BT increases) and the fluid mixing is enhanced ( ζ BT increases). This impact is physically intuitive in a fracture-dominated flow setting with a large permeability contrast between fracture and matrix (here the ratio is approximately 2× 10 5 ). Less smooth fracture surfaces (i.e. higher JRC) and less-weathered joints (i.e. higher JCS) result in a slower breakthrough, mainly as a consequence of lower fracture permeability from less shear-induced dilation. Interestingly, I find out that joint strength due to weathering exerts a more prominent impact than joint asperity in the macroscopic transport behavior in fractured rocks. Flow-transport coupling strength parameterR plays a less impactful role in controlling the variation oft BT andζ BT becauseofitslimitedpresenceinthemathematicalmodelofthecoupledproblem; see concentration-dependent viscosity equation (3.16). The diffusion coefficient D d , which is varying over three orders of magnitudes, has the most significant impact on the transport behavior. A more diffusive tracer propagates slower in the fractures, breaks through later, and creates a longer mixing zone as the tracer diffuses across matrix-fracture interfaces and spreads through the aquifer. I observe a positive correlation between t BT and ζ BT because the longer the tracer stays in the domain before the breakthrough, the longer is the contact duration between the tracer and in-situ fluids, which enhances mixing. We further probe the effect of poroelastic coupling on transport metrics in terms of spatial distribu- tions of the tracer cloud. In practical applications, there are several stress-altering processes such as drilling of a well, hydraulic fracturing, re-fracturing of a fractured well, fluid injection or production, heat injection or production during geothermal operations and steam injection, and natural or an- thropogenic seismicity. The in-situ natural fracture network may also evolve (opening or healing of fractures) due to geochemical processes involving mineral dissolution and precipitation, which adds extra complexities to the evolution of the stress field. As I have shown, consideration of different degrees of geomechanical coupling strength is going to impact the extent of the injectant mixing zone as well as how fast the injectant breaks through at neighboring wells. The extent of the mixing 74 0 D c JRC JCS E R 100% 50% −50% 0 100% 50% −50% 150% 150% D c JRC JCS E R Base Case Base Case 0.5 0.2 0.4 50 0 20 150 1 8 120 2.5 10 0.9 1 8 120 2.5 10 0.2 0.4 50 0 20 150 D d D d t BT ζ BT (a) (b) Figure 3.11: Tornado plot of the percentage change of breakthrough timet BT (panel (a)) and degree of mixing ζ BT (panel (b)) compared to the base case of different sensitivity parameters listed on the y-axis. From bottom to top on the vertical axis: Biot Coefficient α , Fluid Compressibility c f , Young Modulus E, Log Viscosity Ratio R, Joint Roughness Coefficient JRC, Joint Compressive Strength JCS, Diffusion Coefficient D d . Red columns to the left show reduction while blue columns to the right show increase of either breakthrough time or degree of mixing compared to the base case. The end values in each sensitivity bar indicate the corresponding values of each sensitivity parameter, with units shown in Table 3.4. zone and the degree of mixing between the tracer and in-situ fluids determine the efficiency of the injection operation and the quality of model calibration in tracer surveillance studies. In fractured aquifers, when the tracer breaks through at a well, changing the extent or degree of mixing becomes challenging because of the creation of a persistent pathway between the injection and pumping wells. Thus, a common objective in many operations is to maximize the degree of mixing or extend the mixing zone before any breakthrough occurs. Figure 3.12 illustrates how the Biot coefficient, which is a proxy for the degree of consolidation of the rock and also indicates the strength of flow- geomechanical coupling, modulates tracer distribution in the domain. When the aquifer is more consolidated, either naturally or becomes consolidated over time due to production-induced stresses, I observe a reduction of tracer concentration of up to 35% in the central mixing zone around the injector and along major fracture-matrix interfaces. 35% could be significant considering the cost of tracer injection. This observation supports the need to include proper geomechanical coupling while modeling transport in deformable aquifers. Next, I revisit the behavior of the tracer retrieval curve and study how it is regulated by not only the fracture morphology (Figure 3.6) but also the strength of flow-geomechanical coupling 75 0 -20 -30 -35 -10 -25 -15 -5 % (a) (b) 50 m Figure3.12: EffectofBiotcoefficient α onthetracerconcentrationdistribution. Thethickhorizontal line is the cyclic well and the thin black lines are the fracture planes. Difference between two concentration fields for (a) α = 0.5 and α = 0.2, and (b) α = 0.9 and α = 0.2 are shown at the breakthrough time. The loss of tracer concentration is intensified as the flow-geomechanical coupling strength increases with increasing α . represented by E and C f (Figure 3.13). The actual field data is restricted to the tracer withdrawal period and is made dimensionless by normalizing with the peak concentration to allow comparison with my solution. More compressible fluid and more compressible rock suppress the peak of the tracer retrieval curve. After the peak, the decreasing slope portion of the curve (Figure 3.13) is modulated mainly by the fracture morphology (fracture density and orientation), because the flow- geomechanical coupling parameters (E and C f ) have a minor impact on the slope. These effects can provide critical guidance during the calibration of the aquifer model via tracer retrieval data. The middle section of the well, with a denser and more favorably oriented fracture set, records lower tracer peaks due to faster drainage of the tracer. For the ranges used in this sensitivity study, C f induces a larger relative change to the tracer peak than E. It is important to note that matching the simulated tracer profile to the observed profile is not the goal here. Instead, my emphasis is on understanding how different coupling parameters regulate the shape of the tracer retrieval curve. In practice, this understanding can be used to improve the qualitative agreement between modeled and observed tracer curves in terms of their peak concen- tration and slope values. We examine the effect of fluid compressibility and rock stiffness on the temporal evolution of the degree of mixing (see Figure S3 in the supporting information). These 76 ∆t wd c w Midpoint c f =0.4 1/GPa Right-end c f =0.4 1/GPa Midpoint c f =1.0 1/GPa Right-end c f =1.0 1/GPa ∆t wd Midpoint E=10 GPa Right-end E=50 GPa Midpoint E=10 GPa Right-end E=50 GPa Field data Right-end Midpoint (a) (b) (c) c w Figure 3.13: (a) Tracer retrieval curve data from a case study Li et al. (2016). The inset figure shows the location of the well (thick blue line), fracture (thinner green lines), and tracer monitoring stations (midpoint of well and midpoint of the right well segment). The horizontal axis is the time elapsed from the start of the withdrawal period, ∆ t wd , in days. (b) and (c) show tracer retrieval curves during tracer withdrawal periods at different well locations for different values of the fluid compressibility C f (b) and rock Modulus E (c). Lower panels demonstrate the effect of flow- geomechanical coupling (different colors but same line styles) and fracture morphology (same colors but different line styles). All curves are shifted toward a common starting point on the horizontal axis for better visual comparison. 77 observations suggest that the classical definition of the flow-geomechanics coupling strength param- eter Armero and Simo (1992); Tran and Jha (2020), η = α 2 M K dr = α 2 (1− α )(α − ϕ )+ϕC f K dr cannot properly quantify the strength of transport-geomechanics coupling in fractured porous media. This is be- cause macroscopic transport metrics exhibit a non-linear relationship with α . Also, different trends of ζ BT and t BT are observed if I increase η by changing different component parameters. For exam- ple, decreasing C f retards tracer breakthrough and decreasing E accelerates tracer breakthrough, although both lead to increasing η (Figure 3.11). Next, I make a deeper investigation into the effect of poroelastic coupling on fracture permeability via stress field because of the critical importance of fracture permeability in dictating flow and transport behavior in fractured porous media. Here, the effect of poroelastic coupling on the synergy between fracture dynamics and transport is demonstrated through the fracture permeability (k f ) evolution during tracer injection and withdrawal (Figure 3.14). For the sensitivity range of input parameters used, I observe a larger separation of permeability trends in the injection compared to the withdrawal period, mainly because of the shear-induced dilation effect during injection. By the end of the withdrawal period, E shows minimal impact on fracture permeability while C f and α sensitivity yield approximately 5% change in k f . In comparison, across the same range of input, by the end of the injection period or tracer breakthrough, E and α contribute to about 16% variation while C f yields about 25% difference in k f . Rock that is more unconsolidated (higher α ) and softer (lower E) elevates the high end of the permeability range while less compressible pore fluid (lower C f ) inhibits fracture permeability growth. The augmentation of fracture permeability at higher α , higher C f , and lower E explains the decrease in the tracer breakthrough time for these cases in Figure 3.11. Another factor that could be considered is the compliance contrast between the fractures and the host matrix. Fractures are often mechanically more compliant than the host rock. Thus, I perform a sensitivity study where the bulk modulus of fractures is lower than that of the matrix Murdoch and Germanovich (2006). Please refer to Text S4 and Figures S4-S7 in the supporting information for the result of this study. 78 injection withdrawal (a) (b) (c) C f =1(1/GPa) E=10 GPa k f (Darcy) k f (Darcy) k f (Darcy) injection withdrawal injection withdrawal 1000 1200 1400 16001000 1200 1400 1600 1000 1200 1400 1600 0.4 0.2 0 0.8 1 0.6 Figure 3.14: Empirical cumulative distribution function (ECDF) of fracture permeability (k f ) plot- ted for different values of the Biot coefficient α (a), fluid compressibility C f (b), and rock’s modulus E (c). Injection and withdrawal responses are separated by plotting them as solid and dash lines, respectively. 3.8.2. Fingering and fracturing parameters impact the stress state We examine and characterize the impact of flow-transport and transport-geomechanics coupling on the state of stress. This understanding is critical to improve the operational design of a cyclic injection/withdrawal operation, for example, to maintain the mechanical integrity of the formation and to avoid reactivation of natural fractures. Here, I specifically focus on the effect of fracture dynamics and transport parameters on geomechanical quantities. While flow-geomechanics coupling has a more significant first-order impact on the stress state, flow-transport coupling through R (i.e. viscous fingering) also exerts considerable influence on the stress field as shown in Figure 3.15. The effectoffingeringismoreevidenttowardstheendoftheinjectionperiodwhenthetracerhasinvaded parts of the matrix (bottom panels of Figure 3.15). In case of no coupling from transport to flow, i.e. R=0, I observe a flat trend line signifying there is no dynamic synergy between geomechanics and transport. As R increases, curves between the root-mean-square (rms) of stress invariants (I ′ MC,rms and p J 2,rms ) and the transport metrics deviate downward, compared to the R=0 case. For either cases, the discrepancy between R = 0 and R = 2.5 can amount up to 7% for p J 2,rms and5% forI ′ MC,rms . This emphasizes the need for properly considering this flow-transport coupling mechanism, especially in critically stressed rocks close to failure. Near the injection well, where a significant fraction of the tracer is present and the flow-transport coupling (viscous fingering) dominates, I observe lower values of the Mohr-Coulomb failure function. This suggests that shear failure is less likely to happen in regions with stronger flow-transport coupling. 79 During the withdrawal period (dashed curves in Figure 3.15), different R curves start from different points near the bottom right corner of the plot. This horizontal offset between different curves for different R is due to the difference in the degree of mixing ( ζ ) and the average matrix concentration (c m ) at the start of the withdrawal period. Towards the end of the withdrawal period, the R = 0 curve flattens out, again suggesting negligible dependency between stress and transport parameters. However, with the advent of fingering at R>0, the curves deviate as they did during injection, and the impact of flow-transport coupling on the stress state becomes noticeable. We believe that the dependency of the stress state on transport metrics becomes stronger as the transport controlling parameter increases in value. I also believe that the above observation will hold true if viscous fingering is substituted with density-driven convective instability, which is commonly encountered during seawater intrusion into coastal aquifers Huyakorn et al. (1987) and CO 2 -brine mixing during carbon sequestration in saline aquifers Metz et al. (2005). In that case, the Rayleigh-Benard or Rayleigh-Taylor instability may substitute the Saffman-Taylor instability, and the density contrast will substitute the viscous contrast. Finally, I look at how fracture mechanical properties (joint weathering and joint strength) influence the stability of fractures during injection and withdrawal (Figure 3.16 and Figure 3.17). Because JCS and JRC values for an aquifer can be estimated from load testing experiments on core samples, my findings can be used to decide where and how the tracer should be injected to ensure the mechanical integrity of the subsurface. In Figures 3.16 and 3.17, the effective normal and shear stresses on a selected few fracture segments are plotted on the Mohr diagram to show their evolution with time and fracture properties. In a strike-slip faulting regime (i.e. S h < S v < S H where S v is the vertical principal stress), the difference between S h and S H determines the largest Mohr circle radius. Initially, in absence of any pore pressure perturbation, the initial stresses inside the domain balance the tectonic compressions on the boundary, which can be used to draw the initial Mohr circle. The initial stress state of the fracture segments (denoted by black squares in Figures 3.16 and 3.17) lie on the initial Mohr circle. The three squares are at different locations on the circle due to the difference in the fracture orientation. During injection, the normal effective compression σ ′ n acting on the fracture segments decreases, which causes the points to move left toward the static 80 Injection Injection Withdrawal Withdrawal ζ ζ (MPa) (MPa) (a) (b) (c) (d) (MPa) (MPa) c m c m R=0 1 2 2.5 Figure 3.15: Effect of the log viscosity ratio R, a transport parameter, on the stress state. (a) the root mean square of the square root of the second invariant of deviatoric stress tensor p J 2,rms is plotted against the degree of mixing ζ . (b) The root mean square of Mohr-Coulomb stress invariant I ′ MC,rms is plotted against the domain-averaged concentration c m . Evolution during the injection and withdrawal periods are shown via solid and dotted lines, respectively. The direction of time for each period is shown via arrows. The shaded regions are magnified in (c) and (d) for better visualization. 81 frictional failure envelope. Minor movement along the shear stress axis is also observed because the injection-induced expansion of the aquifer applies shear traction on some fractures, the magnitude of which depends on their orientation. In the figure, some points cross the failure envelope because my fracture mechanics model is one-way coupled to the domain equilibrium (section 3.4.5) and does not impose a strong restriction on such crossings. The magnitude of change in the stress state determines the likelihood of slip on the fracture. The amount of shift varies depending on the values of fracture mechanics parameters, JRC and JCS. In a reservoir with over-pressured pore fluid or more pronounced stress heterogeneity, these 4-5 MPa stress changes could easily induce slip on fractures. In Figure 3.16, fractures with higher JCS, implying less weathered joints, have larger stress changes and thus are more likely to fail in the shear mode. During tracer withdrawal, the less weathered fractures move farther away from the failure envelope, suggesting that they are less likely to slip during this period. For fractures of different JRC, i.e. different resistances to shear slip, the magnitude of poroelastic stress shifts are equally pronounced although the sensitivity spread of stress output among JRC values investigated is smaller than that of JCS (Figure 3.17). This suggests that, compared to the compressive strength, the resistance to shear has a less pronounced effect in governing fracture stability. Fractures with smaller JRC tend to move closer to the failure envelope during injection and farther during withdrawal. This implies that a fracture with lower shear strength is more likely to fail during injection. 3.9. Conclusions We proposed a novel multiphysics modeling framework to investigate the rich physical interchange between flow, geomechanics, transport, and fracture dynamics in naturally fractured aquifers. We show that successful integration of transport-geomechanical coupling in field-scale modeling studies of flow through fractured rocks can improve the prediction of macroscopic transport parameters: the tracer breakthrough time, the average and maximum concentrations of the tracer at a control plane, and the evolution of the degree of mixing with time. Additional contributions of the study are as follows: 1. IutilizedEmbeddedDiscreteFractureModeling(EDFM)tomodelthefracturedporousmedia 82 x x x Injection Withdrawal Figure 3.16: Evolution in the stress states (shear stress τ versus effective normal stress σ ′ n ) of three representative fracture segments is shown using a Mohr circle diagram. Three timesteps are shown– initial (black square), end of injection (solid circles), and end of withdrawal (asterisks)–for three different JCS values. The static failure envelope is plotted in red with the static frictional coefficient of 0.35. The shaded region in the top left diagram is shown in the right diagram with a zoom-in. The arrows indicate the direction of change in the stress, from the initial state, during injection and withdrawal. x x x Injection Withdrawal Figure 3.17: Stress states of three representative fracture segments are plotted in a Mohr diagram at different times: initial (black square), end of injection (solid circles), and end of withdrawal (asterisks) for different JRC values. The static failure envelope is plotted in red with the static frictional coefficient of 0.35. The shaded region in the top left diagram is shown in the right diagram with a zoom-in. The arrows indicate the directional change of stress state from the initial state during injection and withdrawal. 83 and implemented an improved Bandis model to simulate fracture mechanics. my EDFM- based framework allows modeling multiple randomly oriented intersecting fractures at a lower computational cost than many discrete fracture methods that require a conforming mesh. 2. I characterized the effect of geomechanical coupling on key transport metrics and quantified theimpactofeachcouplingparameterontracerbreakthroughanddegreeofmixing. Ifindthat parametersgoverningthefracturedynamics(JCS,JRC)haveahigherimpactonthespreading and mixing of tracer compared to parameters dictating the strength of flow-geomechanical coupling (Biot coefficient, fluid compressibility, rock’s bulk modulus). 3. Tracer retrieval curve characteristics e.g., the peak concentration, are modulated by both fracture morphology and the strength of geomechanical coupling. 4. Flow-transport coupling, via viscous fingering and fracture dynamics, alter the stress state and fracture stability in the domain. This is important to consider in critically stressed formations close to failure. 5. I quantified the impact of flow-geomechanics coupling strength on the evolution of fracture permeability. Our results suggest that the coupling between flow, transport, and geomechanics will be necessary to accurately model the fate of an active tracer (which dynamically modifies the viscosity or density of the flow) during reservoir surveillance, contaminant remediation, or fracture characterization studies in deformable rocks. For shallower reservoirs (for example, in groundwater applications), I hypothesize that the coupling between geomechanics and transport can become stronger when the rock is more stress-sensitive due to a lower degree of consolidation (higher Biot coefficient) and lower bulk modulus. The findings of this study could help motivate and inform the future design and deployment of sensors that can monitor small (<10%) changes in stress/strain fields that may arise from coupling with the flow and transport processes. An example of such high fidelity and precision measurement technique is Distributed Strain Sensing (DSS), which uses fiber optics technology to reliably detect small changes in stress/strain fields around a wellbore. 84 Our modeling framework assumes that the feedback from fracture slip to the matrix stress is via changes in fracture properties only. I neglect fracture slip-induced stress relaxation, which can become important in case of fractures/faults that rupture over a long distance or grow in length from their tips. A future extension of my work includes modeling the propagation and healing of fractures to capture the feedback from fracture deformation to stress. Extending the framework to enable the generation of synthetic microseismic events can also create new possibilities for assimi- lating microseismic monitoring data and calibrating the proposed multiphysics model. Finally, the proposed modeling framework can support multiphysics joint inversion endeavors Jha et al. (2015) at monitoring sites with multiphysics data from InSAR/GPS satellites, ground penetrating radar, and wells to estimate uncertain hydraulic and poroelastic properties of the subsurface. 85 CHAPTER 4 An Embedded Discrete Fracture–Phase Field Model to Quantify Fracturing-induced Pore Pressure 4.1. Abstract Propagation of fractures in fluid-saturated rocks can induce changes in the pore pressure of the surrounding rock matrix due to the poroelastic coupling. Pre-existing fractures can amplify these poroelastic signals, which can then be used to infer the geometry and hydraulic properties of the propagating fracture in real time if I can accurately model the processes at the appropriate field scale. By integrating Embedded Discrete Fracture Modeling (EDFM) and Phase Field Modeling (PFM) in a unified framework, I propose a method to simulate the dynamic interaction of fractures and the pressure and stress change signals associated with individual interactions. We quantify the impact of poroelastic coupling on fracture dynamics, i.e., nucleation, propagation, and interaction between fractures. I characterize the evolution in fracture length, propagation speed, and pressure for different fracture morphology (orientation and spacing) and poroelastic moduli. To establish a novel pressure-based fracture monitoring method for field applications, I apply the proposed framework to quantify the induced pressure signals of multiple fractures propagating in succession, where the signals are observed in monitoring wells in the vicinity of the well undergoing hydraulic fracturing. The proposed fracture modeling framework unifies and synthesizes the best features of EDFM and PFM. 4.2. Introduction The investigation of fracture-dominated flow processes in deformable porous media has remained a complex research topic, especially from the simulation viewpoint. This is due to the intri- cate interplay between multiple underlying processes: fluid flow in the rock matrix and fracture, fluid leak-off into fractures, poroelasticity, fracture dynamics, and damage to the rock’s structural frame Bai and Elsworth (2000a). During fluid-driven fracture propagation, the inter-process cou- pling is more evident because the fracture permeability dictates the fluid pressure distribution in the 86 domain, which impacts the stress state that determines the possible areas for nucleation of new frac- tures and growth of existing fractures Bai et al. (1993); Gangi (1978); Baghbanan and Jing (2008); Svenson et al. (2008); Rong et al. (2016); Norbeck et al. (2016); Manjunath et al. (2023). Frac- ture propagation perturbs the geomechanical state of the surrounding region and nearby fractures, e.g., causing shear slip and dilation of existing fractures Baghbanan and Jing (2008); Tran and Jha (2021) or closure of the fractures that fall in the stress shadow of a growing fracture. These effects introduce spatial heterogeneity and anisotropy in the fluid flow and deformation response of the rock during injection or production of fluids Bai and Elsworth (2000a); Bubshait and Jha (2019). Pre- diction and control of those effects require a robust numerical model to simulate fracture dynamics, from fracture nucleation and propagation to interaction with other fractures, in a poroelastic frame- work. Specific challenges include determination of the fracture path and simulating fracture kinking, branching, or arrest when a propagating fracture encounters existing geomechanical discontinuities or other active fractures in the rock. Measuring and analyzing the fracturing-induced signals has beenthecenterpieceofreal-timefracturemonitoringtechnologytounderstandandcontrolhydraulic fracturing in the field Hull et al. (2017); Richter et al. (2019); Liu et al. (2020). Examples of such technology include Distributed Acoustic Sensing (DAS), Distributed Temperature Sensing (DTS), microseismic monitoring, chemical tracer test, and pressure transient analysis. These signals are processed to deduce fracture characteristics, such as fracture length and height, orientation relative to principal stresses, and fracture conductivity, for quantifying the stimulated rock volume (SRV). Recently, analysis of pressure jumps measured in an offset well, different from the well undergoing fracturing, has been proposed as a method to diagnose fracture morphology Kampfer and Dawson (2016b); Roussel and Agrawal (2017b); Manjunath et al. (2023). The idea is that a fluid-saturated porous medium, subjected to a sudden external load due to opening of a fracture, experiences a pore pressure increase due to undrained deformation Wang (2000), which then progressively diffuses through the domain. In a low-permeability formation, pressure diffusion is slow enough such that the pressure jump induced by fracture growth can be observed at a point far from the growing fracture. Inthisstudy, mygoalistodevelopanumericalframeworktomodeltheeffectoffracturepropagation 87 and interaction on pressure-based monitoring of hydraulic fractures in fluid-saturated poroelastic media. One of the novelties of the framework is that it synthesizes the best features of two powerful methods for modeling fractures–embedded discrete fracture method and phase field method–as discussed below. The proposed framework can help us address questions such as how poroelasticity impacts the dynamic evolution of fracture and how to quantify the effect of fracture propagation on the induced pressure signal. 4.2.1. Brief review of discrete fracture modeling methods There are two popular categories of numerical methods for fracture modeling: discrete fracture network modeling and continuum damage modeling. Across the fracture surface, displacement is allowed to be discontinuous in the former approach while in the latter approach, the displace- ment field remains continuous everywhere. Discrete fracture methods are further categorized into conforming mesh methods and embedded mesh methods. Conforming mesh methods require an unstructured mesh with fractures along mesh element boundaries to accommodate fractures of dif- ferent orientations and intersecting fractures Goodman et al. (1968); Karimi-Fard et al. (2004). In embedded mesh methods, fractures are explicitly represented in a matrix mesh that remains static during the fracture growth. The idea is that fractures can grow arbitrarily and independent of the background mesh such that a complex, tortuous fracture can be modeled using a simple structured mesh. The extended finite element method (XFEM) Moës et al. (1999); Schwenck et al. (2015) and the embedded discrete fracture method (EDFM) Lee et al. (2001b); Hajibeygi et al. (2011) are two examples. In EDFM, fractures are represented using a structured gridding scheme, sepa- rate from the host matrix grid. Special non-neighboring connections are used to model the fluxes between fracture and matrix Tran and Jha (2021). Another class of discrete fracture method is the phase field method (PFM) Bourdin et al. (2008); Miehe and Mauthe (2016); Wu et al. (2020) where fractures are approximated as bands of finite thickness characterized by a phase field param- eter φ ranging from 0 (intact rock) to 1 (fracture). The sharp discontinuities are regularized by a diffused phase field variable, easily enabling this approach to deal with crack nucleation, fracture interaction (like merging and branching), and multiple crack fronts without additional complex- ity in the numerical treatment. Many applications of phase field modeling in hydraulic fracturing 88 have been proposed, e.g., see Ni et al. (2020); Yi et al. (2020); Chukwudozie et al. (2019). The im- pact of fracture spacing, reservoir permeability, and geomechanical heterogeneity on the interaction and propagation of multiple fractures has been studied Chukwudozie et al. (2019). PFM has been used to model dry cracks in functionally graded materials Natarajan et al. (2019). Coupling of the phase field variable with fluid flow and geomechanical processes has been implemented Wick et al. (2016); Mikelić et al. (2015), with a detailed formulation presented in Heister et al. (2015). Other studies captured fracture interaction complexity on variably saturated porous media Santillán et al. (2018); Cajuhi et al. (2018). One study Ni et al. (2020) improved the approach by resolving the issue of fracture reversibility and maintaining only tension-induced fracture with the inclusion of a monotonously increasing maximum energy density term in the phase field equation. A method was proposed to use the Darcy flow model for both matrix and fracture domains Yi et al. (2020), removing the need to convert the phase field value into a fracture aperture value for a lubrica- tion/Poiseuille flow model in the fracture. Interfacial damage between the propagating cracks and the rock matrix was considered in a coupled flow-geomechanics framework Xia et al. (2017), open- ing the possibility of modeling complex crack geometry in regular grids of voxels obtained from experimental imaging data on rock samples. Each method has its own advantages and disadvantages when compared to other methods. The comparison is often expressed in terms of the trade-off among practical concerns such as meshing requirement, capacity to incorporate fracture geometry and propagation related data, accuracy of stress and pressure solution, and the simulation runtime. Different flavors of one method or combi- nations of two methods exist to address these concerns in real-world applications. For example, a mixed EDFM-XFEM framework with discrete fractures has been proposed recently Ren and Younis (2021) to simulate fluid-driven hydraulic fracturing under prescribed injection and depletion scenar- ios relevant to petroleum recovery applications. my EDFM-PFM framework falls in that category of mixed methods and is unique because of its strong emphasis on modeling the fracturing-induced pressure and change in poromechanical parameters that evolve with the damage caused by fracture propagation. One of the challenges in PFM is that to obtain numerically robust and stable crack propagation, the mesh element size needs to be sufficiently refined to satisfactorily resolve the dam- 89 aged zone around the fracture, leading to potentially exorbitant computational demand. To reduce this demand and improve the convergence rate and the solution accuracy, I use a fourth-order crack energy density function. I use the isogeometric analysis (IGA) and adaptive h-refinement schemes to significantly save on the computational cost Goswami et al. (2020). 4.3. Mathematical Model We make the following assumptions, which are appropriate for the physical problem at hand: (1) slightly compressible pore fluid, (2) small strain deformation, (3) quasi-static equilibrium, (4) linear elasticity, (5) single-phase fluid flow, (6) isotropic material behavior, (7) isothermal conditions, and (8) no body forces, e.g., gravity. Equations of quasi-static equilibrium, solid and fluid mass conser- vation, and phase field evolution govern the system’s behavior in space and time under the imposed initial and boundary conditions. I formulate the problem with the following primary unknowns: fluid pressure in the matrix domain p m , fluid pressure in the fracture domain p f , displacement in the matrix u, and the phase field variable φ. The superscripts m and f are used for the matrix and fracture, respectively. The applicable constitutive laws for the solid, fluid, and the phase field variable provide closure for the system of equations. 4.3.1. Constitutive Models We approximate a fracture as a material region of finite thickness characterized by the phase field parameter φ, changing from 0 to 1. φ=0 represents undamaged rock while φ=1 represents fully cracked material and φ increasing from 0 to 1 indicates a progressively more damaged material. To accommodate the mechanical weakening effect of fracture presence and propagation on the rock frame, the magnitudes of key poroelastic parameters are defined as functions of the extent of damage in the domain. Because the phase field variable is a representative proxy for the extent of fracture propagation,φ can play a role analogous to that of a damage variable in continuum damage mechanics-based fracture modeling methods. Therefore, in my formulation, the Biot coefficient α , which controls the degree of rock consolidation, depends on phase field φ: α (φ)=1− K dr (φ) K s =1− (1− φ) 2 K 0 dr K s (4.1) 90 where K 0 dr is the undamaged drained bulk modulus, K dr is the drained bulk modulus, and K s is the solid bulk modulus. If φ = 1 then α = 1 (in fracture) and if φ = 0 then α = 1− K 0 dr Ks (in porous medium or intact rock). Assuming that Poisson’s Ratio ν does not change with the extent of damage, I modified the drained Young’s modulus E to be a function of the phase field variable via its relationship with the drained bulk modulus: E(φ)=2K dr (1− 2ν )(1+ν )=2(1− φ) 2 K 0 dr (1− 2ν )(1+ν ) (4.2) The effective volumetric stress σ ′ v is related to the volumetric strain ϵ v as σ ′ v =K dr ϵ v . In my EDFM formulation, the effective fracture porosity is small (on the order of 0.001) and is defined as ϕ f = L f w f A m b where L f is the fracture length, w f is the fracture aperture, and A m b is the cross sectional area of the matrix mesh element assigned to the fracture segment Xu (2015). To account for the fracture propagation effect, I model the change of fracture porosity as a function of the phase field variable Yi et al. (2020): ϕ f (φ)=ϕ m,0 +(ϕ f,0 − ϕ m,0 ) 1− (1− φ) 2 (4.3) Here, ϕ f,0 andϕ m,0 are the initial fracture and matrix porosities, respectively. The end point values are reasonable: φ = 1 implies ϕ f = ϕ f,0 and φ = 0 implies ϕ f = ϕ m,0 . In my staggered solution scheme, discussed in section 4.4.2, the matrix porosity is updated outside the inner loop with a converged value of the phase field variable. The Biot modulus M becomes a function of the phase field variable because of the dependence of porosity and Biot’s coefficient on φ. 1 M m (φ) =ϕ m (φ)c f + α (φ)− ϕ m (φ) K s , 1 M f (φ) =ϕ f (φ)c f (4.4) where c f is the fluid compressibility from fluid’s equation-of-state. During fracture propagation, the fracture permeability needs special treatment because fracture morphology (i.e., spatial density and direction) and fracture’s hydraulic properties (i.e., aperture) need to be updated based on the 91 evolving stress field and the extent of damage. To update the fracture permeability ( k f ) for newly formed as well as existing fracture segments, the well-established Bandis model is used to update the change of fracture aperture; see Tran and Jha (2020) for details. 4.3.2. Energy Functional Formulation Deformation of an elastic body under load increases the elastic energy. When this value approaches a critical value in a region, it is energetically favorable for the system to decrease the elastic energy throughthecreationofnewfracturesurfaces,increasingthephasefieldvaluetowardsunity. Fracture propagation results from the competition between the bulk energy in the matrix, away from the fracture, and the fracture surface energy. For a coupled poroelasticity problem, the total free energy of a volume domain V (with boundary ∂V) in the presence of a fracture (with Γ f as the fracture surface) consists of four terms: stored elastic strain energy in the solid Ψ st , surface fracture energy Ψ f su , fluid energy in the matrix Ψ m fl , and fluid energy in the fracture Ψ f fl . The energy functional is a function of the primary unknowns: displacement field u, pressure in matrix p m , pressure in fracture p f , and phase field variable φ. Ψ u,φ,p m ,p f =Ψ st (u,φ)+Ψ f su (u,φ)+Ψ m fl (u,p m ,φ)+Ψ f fl (p f ,φ) (4.5) The stored elastic strain energy in the solid is expressed as a function of displacement and phase field variable: Ψ st (u,φ)= Z V σ ′ :ϵ dV (4.6) Here,ϵ =(∇u+∇u T )/2 is the infinitesimal strain tensor, and σ ′ =2Gϵ +λ tr(ϵ )I is the effective stress tensor, with I being the identity tensor. Here, tr is the trace operator. The two Lame constants are defined as λ = νE (1+ν )(1− 2ν ) and G = E 2(1+ν ) . For an isotropic elastic body, the elastic energy term can also be defined as: Ψ st = λ 2 [tr(ϵ )] 2 +G tr(ϵ 2 ). As fracture propagates and the damaged region expands, matrix stiffness deteriorates, requiring us to modify equation (4.6) with the inclusion of the stress degradation function g(φ) and strain energy split: Ψ st (u,φ)= Z V g(φ)Ψ + st +Ψ − st dV = Z V [g(φ)(σ ′ :ϵ ) + +(σ ′ :ϵ ) − ]dV (4.7) 92 The stress degradation function accounts for the change in stored strain energy due to crack evolu- tion, i.e., how energy decays while generating a new crack surface. This function needs to satisfy the following conditions: g(0) = 1 to represent no damage in the intact zone, g(1) = 0 to describe the total damage in the fractured zone, g(φ)>0 if φ̸=1 to avoid a fully damaged state, g ′ (0)̸=0 be- cause crack nucleation initiates from a uniform initial state, and g ′ (1)=0 because energy converges to a finite value for a newly formed fracture. Many forms of the stress degradation function have been suggested in the literature. The quadratic form isg(φ)=(1− κ )(1− φ) 2 +κ , whereκ is a small dimensionless numerical parameter (O(10 − 6 )) to satisfy the well-posedness condition by avoiding the singularity of the stiffness matrix in the fully damaged region. A more popular choice is the cubic form of the stress degradation function, g(φ)=κ [(1− φ) 3 − (1− φ) 2 ]+3(1− φ) 2 − 2(1− φ) 3 . This choice modulates the degree of stiffness degradation prior to the start of cracking, mitigating the unrealistic nonlinearity of the stress-strain curve Goswami et al. (2020). In this study, I chose g(φ)=(1− φ) 2 +κ, g ′ (φ)=− 2(1− φ) (4.8) which ensures the convergence of the function to the functional of the discrete crack as the crack becomes infinitesimally thinner. With this form of the degradation function, φ=0 implies g(φ)= 1+κ ≈ 1 (intact rock) and φ = 1 implies g(φ) = κ ≈ 0 (fracture). The derivative conditions are also satisfied: g ′ (0)̸=0 and g ′ (1)=0 Thetensileandcompressivecomponentsofthestoredelasticenergyarecalculatedfromtheprincipal strain tensor ϵ Yi et al. (2020). In 2D, Ψ ± st (ϵ )= λ 2 ⟨ϵ 1 +ϵ 2 ⟩ 2 ± +G ⟨ϵ 1 ⟩ 2 ± +⟨ϵ 2 ⟩ 2 ± , (4.9) where the Macaulay bracket is defined as ⟨x⟩ ± = |x|± x 2 . I utilize an isotropic model for the elastic equation, meaningthereisnoenergysplitconsideredintheelasticequation(Ψ − st =0). Forthephase field governing equation, an anisotropic formulation is used such that only the tensile component of the stress field drives the change in the phase field variable (tensile fracture propagation). Thus 93 the stress degradation function is applied only to the tensile component of the elastic stress tensor (Ψ + st ) Goswami et al. (2020). Equation (4.7) then becomes: Ψ st (u,φ)≈ Z V g(φ)Ψ + st dV (4.10) Next, thefracturesurfaceenergyisexpressedasanarealintegrationofthecriticalenergyreleaserate G c (i.e., fracture toughness) over the fracture surface Γ f . This integral is numerically intractable, thus it is modified into a volume integral by using the fracture surface density function γ , a function of the phase field and the phase field gradient: Ψ f su (u,φ)= Z Γ f G c dΓ f ≈ Z V G c γ (φ,∇φ)dV (4.11) where γ is the second-order fracture surface energy density function for which I select the following form: γ = φ 2 2l 0 + l 0 2 (∇φ) 2 , (4.12) where l 0 is the length scale parameter (a material parameter) that acts as a proxy for the physical thickness of the diffused interface between a fracture and the intact rock surrounding the fracture. The gradient term in equation (4.12) regulates the amount of diffusion of the fracture width and, thereby, the numerical smoothness across the fracture interfaces. In the limiting case, l 0 → 0, and the gradient term disappears resulting in a sharp interface. For an intact rock with no fracture surface energy, φ=0, which implies γ =0. The phase field evolution is described by the microscale force balance equation Shen et al. (2019). Ignoring the effect of local inertia at the crack tip and external micro forces (such as chemical bond breaking, surface adsorption), the microscale traction vectorH m is related to the internal microscale force K m Arriaga and Waisman (2018). ∇· H m =K m . (4.13) The contribution of the fluid energy to the matrix domain or the fracture domain is a function of 94 displacementu and the incremental fluid volume content ζ , which can be written as in Mikelić et al. (2015): Ψ m fl (u,ζ m )= (p m ) 2 2M m = M m 2 (ζ m − αϵ v ) 2 (4.14) Ψ f fl ζ f = p f 2 2M f = M f 2 ζ f 2 (4.15) From equation (4.5), we rewrite the time derivative of the total energy functional in terms of its constituents and use the chain rule: ˙ Ψ= ˙ Ψ st + ˙ Ψ f su + ˙ Ψ m fl + ˙ Ψ f fl = ∂Ψ st ∂ϵ : ˙ ϵ + ∂Ψ st ∂φ : ˙ φ + ( ∂Ψ f su ∂∇φ :∇˙ φ+ ∂Ψ f su ∂φ : ˙ φ ) + ∂Ψ m fl ∂ϵ : ˙ ϵ + ∂Ψ m fl ∂φ : ˙ φ+ ∂Ψ m fl ∂ζ m : ˙ ζ m + ( ∂Ψ f fl ∂ϵ : ˙ ϵ + ∂Ψ f fl ∂φ : ˙ φ+ ∂Ψ f fl ∂ζ f : ˙ ζ f ) (4.16) where the overdot symbolizes the time derivative. 4.3.3. Clausius-Duhem Inequality As per the second law of thermodynamics, the total mechanical dissipation of an isothermal pro- cess should be non-negative, also known as the Clausius–Duhem inequality. In terms of the time derivative of entropy (˙ e) and the total free energy functional Ψ : ρ m ˙ e− ˙ Ψ ≥ 0 (4.17) Based on the internal energy formulation Shen et al. (2019) and combining with equation (4.16), I can expand this inequality by collecting the terms: h σ : ˙ ϵ +p m ˙ ζ m +p f ˙ ζ f +H m ·∇ ˙ φ+K m ˙ φ i − h ˙ Ψ st + ˙ Ψ f su + ˙ Ψ m fl + ˙ Ψ f fl i ≥ 0 (4.18) 95 Collecting terms that share a common factor, σ − ∂Ψ st ∂ϵ − ∂Ψ m fl ∂ϵ : ˙ ϵ + p m − ∂Ψ m fl ∂ζ m ˙ ζ m + " p f − ∂Ψ f fl ∂ζ f # ˙ ζ f + " H m − ∂Ψ f su ∂∇φ # ·∇ ˙ φ + " K m − ∂Ψ st ∂φ − ∂Ψ f su ∂φ − ∂Ψ f fl ∂φ − ∂Ψ m fl ∂φ # ˙ φ≥ 0 (4.19) The time derivatives of variables can take arbitrary values for different thermodynamic processes. Hence, the coefficients of the time derivative terms must be zero, which yields the following set of equations: σ − ∂Ψ st ∂ϵ − ∂Ψ m fl ∂ϵ =0 p m − ∂Ψ m fl ∂ζ m =0 p f − ∂Ψ f fl ∂ζ f =0 H m − ∂Ψ f su ∂∇φ =0 K m − ∂Ψ st ∂φ − ∂Ψ f su ∂φ − ∂Ψ f fl ∂φ − ∂Ψ m fl ∂φ =0 (4.20) Substituting from equations (4.10), (4.11), (4.14), and (4.15), I retrieve the constitutive equations for the coupled poroelasticity–phase field problem, i.e., the effective stress equation, poroelasticity equationsformatrixandfracture,andtherelationsbetweenmicroscaleforceandphasefieldvariable: σ = g(φ) ∂Ψ + st ∂ϵ +M m (ζ m − αϵ v )(− α )=g(φ)(σ ′ )− αp m I (4.21) p m = M m (ζ m − αϵ v ) (4.22) p f = M f ζ f (4.23) H m = G c l 0 ∇φ (4.24) K m = g ′ (φ)Ψ + st + G c φ l 0 + (p m ) 2 2 ∂(1/M m ) ∂φ + (p f ) 2 2 ∂(1/M f ) ∂φ (4.25) 96 To linearize equation (4.25), I recall the inverse Biot modulus expressions for matrix and fracture domains in equations (4.4). Substituting the stress degradation function in the matrix equation gives, 1 M m = ϕ m c f + α − ϕ m K s =ϕ m c f + 1− ϕ m K s − g(φ) K 0 dr K 2 s (4.26) The derivative of Biot modulus with respect to the phase field variable for matrix and fracture domains becomes ∂(1/M m ) ∂φ = c f ∂ϕ m ∂φ − 1 K s ∂ϕ m ∂φ − K 0 dr K 2 s ∂g ∂φ , (4.27) ∂(1/M f ) ∂φ = c f ∂ϕ f ∂φ (4.28) From equation (2.8), the matrix porosity is independent of the phase field variable. Substituting from equation (4.8) and equation (4.3), the expressions for matrix and fracture domains are, ∂(1/M m ) ∂φ = 2(1− φ) K 0 dr K 2 s , (4.29) ∂(1/M f ) ∂φ = 2(1− φ)c f (ϕ f,0 − ϕ m,0 ) (4.30) Substituting these into equation (4.25), the internal micro-force balance becomes K m =− 2(1− φ)Ψ + st + G c φ l 0 +(p m ) 2 (1− φ) K 0 dr K 2 s +(p f ) 2 (1− φ)c f (ϕ f,0 − ϕ m,0 ) (4.31) 4.3.4. Final Governing Equations Linear momentum balance and fluid mass continuity govern the behavior of the coupled fluid flow– geomechanical system: ∇· σ =0, ∂ζ m ∂t +∇· v m =0, ∂ζ f ∂t +∇· v f =0, (4.32) 97 where v m and v f are the Darcy fluid velocity vectors in matrix and fracture, respectively, and I express them in terms of p m and p f using Darcy’s law. I substitute equations (4.24) and (4.31) into equation (4.13) to obtain the phase field evolution equation. The set of final governing equa- tions is obtained by substituting equation (4.21) into the linear momentum balance equations and equations (4.22) and (4.23) into the fluid mass continuity equations: ∇· σ =∇· [g(φ)σ ′ − α (φ)Ip m ]=0 ∂ ∂t p m M m (φ) +α (φ)ϵ v −∇· k m µ ∇p m =ψ mf +ψ mw ∂ ∂t p f M f (φ) −∇ f · k f µ ∇p f =ψ fm +ψ ff +ψ fw ∇· (G c l 0 ∇φ)=− 2(1− φ)Ψ + st + G c φ l 0 +(p m ) 2 (1− φ)(K 0 dr /K 2 s ) +(p f ) 2 (1− φ)c f (ϕ f,0 − ϕ m,0 ) (4.33) where theψ terms indicate the mass flux transfer between subdomains, i.e., from fracture to matrix (superscript mf), matrix to fracture (superscript fm), fracture to fracture (superscript ff), well to matrix (superscript mw), and well to fracture (superscript fw) Tran and Jha (2020). The flux transfer functions are often expressed as a product of a connectivity index (CI) of the interface betweenthetwosubdomains, themobilityofthefluidcrossingtheinterface, andthedifferenceinthe pressure of the two subdomains, which resembles the two-point flux approximation model commonly used in the finite volume method for simulating fluid flow (LeVeque, 2002; Jha and Juanes, 2014). I follow the sign convention of normal stresses being positive under tension. Above, ∇ f is the nabla vector operator in the fracture subdomain, that is, divergence along the 1-D fracture mesh embedded in a 2-D matrix mesh. We define H + =2Ψ + st − (p m ) 2 K 0 dr K 2 s +c f (p f ) 2 (ϕ m,0 − ϕ f,0 ). (4.34) To realistically model the non-healing effect of the fracture creation, the phase field variable should not decrease over time. In other words, φ t+δt ≥ φ t . To enforce this irreversibility constraint in the constrained energy minimization approach, one option is to add a quadratic penalty term to the 98 energy functional, i.e., the exterior penalty method. I modify the history functional η + with the max condition to ensure the irreversibility constraint η + = max[⟨H + ⟩ + ]. For a rock matrix with spatially uniform fracture toughness, the governing equation for the phase field becomes − G c l 0 ∇ 2 φ+ G c φ l 0 =(1− φ)η + (4.35) 4.3.5. Hybrid EDFM-PFM Fracture Mechanics Model WeextendapreviousimplementationofEDFMTran and Jha(2020), whichwaslimitedtofractures offixedlength, topropagatingfractures. Iformulateatrackingalgorithmtodeducetheplacementof new fracture tips, which provides a basis to generate new fracture segments based on the converged solution of the phase field at each time step. In my EDFM formulation, the lower dimensional fracture domain consists of a set of fracture segments with existing endpoints. If the fracture propagation criteria are met, I assume that propagation of existing fracture happens only at the twooutermostendpoints, notatanyintermediatepointalongthebodyofthefracture. Newfracture tips are iteratively selected via a series of screening steps Santillán et al. (2018). This updates the fracture path Γ f , which is composed of a set of tip location vectors, i.e.,{x tip }. • Theconvergedphasefieldsolutionisscreenedtonarrowdownthepossibletipcandidates( ξ tip ) by selecting cell-centered elements in which φ exceeds a certain damage threshold (φ th =0.9). • Based on the shortest distance approach, these candidates are grouped corresponding to each existing tip i that tip i may propagate from (ξ i tip ). • For each existing fracture tip i that possesses a set of potential location to propagate to (ξ i tip ̸= ∅), I impose the trajectory of the new fracture segment such that the directional gradient of φ is minimum, ensuring that the new fracture path is formulated by minimizing the total energy functional around tip i: x i,new tip =argmin ξ i tip h ∇φ(x i,new tip )· (x i,new tip − x i,old tip ) i (4.36) 99 where∇φ(x i,new tip ) is the gradient of φ at the new tip. • The fracture permeability and fracture pressure of the newly-formed segments are assigned similarly to the segments they propagate from. • Following the completion of the tracking algorithm routine after each time step, I update the connectivity index (CI) of the non-neighboring connections (NNC) in the calculation of the fluid flux transfer functions ( ψ terms in equation (4.33)). This is a feedback from localized fracture propagation (PFM) to global pressure distribution in the domain (EDFM). Thus, the evolution of CI represents the geometrical complexity of the dynamic communication between the three subdomains–matrix, fracture, and well–during fracture propagation episodes. 4.3.6. Boundary and Initial Conditions Solution of the coupled problem requires suitable boundary conditions for the three sub-problems: flow, deformation, and fracture propagation. Referring to the physical model in Figure 4.4 and Figure 4.16, I apply the following boundary conditions: zero normal displacement (roller condition) and fixed pressure (zero excess pressure) on all boundaries. The wells located inside the domain are controlled by a prescribed injection rate q w over a prescribed time duration. The unit of q w is the fluid volume per unit time per unit length of the medium in the third dimension. Initialization of a coupled flow-geomechanics problem, especially in low permeability rocks with highly conductive fractures, is a challenging task. The goal of initialization is to achieve mechanical and hydrostatic equilibria such that the initial strains and fluid fluxes inside the domain are zero in absence of sources/sinks, i.e., wells. Spatially uniform fields of pressure, stresses, and displacements, e.g., zeroexcesspressure, zerodisplacements, andtheinitialstressesequaltotheboundarytractions, are prescribed as initial conditions that are close to, but not exactly at, the equilibrium. During the initialization period of the simulation, I solve the (discretized) governing equations with the prescribed boundary compressions and known initial values of rock-fluid properties to achieve the state of poromechanical equilibrium. The wells are inactive during the initialization period. The time step during the initialization period is based on the hydraulic diffusion time scale to attain 100 pressure equilibrium in the matrix and fracture domains. The state variables, e.g. pressure, stress and phase field, at the end of initialization become the reference state. Stress-dependent properties such as porosity and permeability experience perturbations during initialization and change from their initial values to the reference values. The injection well starts operating at the end of the initialization period. Initialization of the phase field sub-problem requires modeling of pre-existing cracks. Several op- tions are available to simulate the initial set of fractures: mesh-induced (via imposing discrete discontinuities in the mesh), phase field-dependent (via setting Dirichlet boundary conditions), or history functional (via establishing a local strain-history function) Goswami et al. (2020). The ini- tial strain-history function, η 0 , can be defined for all points in the domain that are finitely close to any line Γ f representing the initial discrete fracture set Borden et al. (2012): η 0 (x)= BG c 2l 0 1− d(x,Γ f ) n t l 0 , (4.37) where function d is defined as d(x,Γ f )≤ n t l 0 , parameter n t captures how thin the initial fracture set is, and parameter B captures how diffused the initial crack set is. I set n t =0.5 and B =1000. Points that are far away from any initial fracture line are assigned η 0 =0. This completes the definition of the coupled problem. The governing equations to be solved are equations (4.33) and (4.35) subject to the initial and boundary conditions mentioned above. The equation-of-state for the fluid compressibility, stress-strain relation, well connectivity relation, rela- tions for mass transfer fluxes between matrix, fracture and well, dependence of poroelastic parame- ters on the phase field variable, and fracture mechanics provide the necessary closure of the system of equations. As mentioned earlier, I consider the matrix and fracture fluid pressures ( p m and p f ), displacement vector u, and phase field variable φ as the primary unknowns of the problem. The effective stress σ ′ , porosity ϕ , matrix and fracture Darcy velocities (v m and v f ), fracture aperture a h , poroelastic parameters (α , M, and E), and matrix and fracture permeabilities (k m and k f ) are the secondary variables which I update explicitly by post-processing the primary variables. 101 4.4. Numerical Model Solvingthecoupledproblemrequiresanumericaldiscretizationschemethatcanhonorthecontinuity requirementsoftheprimaryandsecondaryvariableswhileallowingformultiplescalesinthediscrete systemandasufficientresolutionofthediscretefracturepropagationpath. Computationalefficiency and robustness are the other two requirements of the numerical model because I intend to apply the model to reservoirs with geometrically complex fractures/faults. To satisfy these requirements, I use the finite volume method to discretize the flow equations and the Galerkin finite element method to discretize the mechanics and phase field equations Jha and Juanes (2014); Tran and Jha (2020). 4.4.1. Weak Forms and Discretization We discretize the domain using an [n d ]-dimensional uniform Cartesian mesh (n d = 2 for the 2D case). The fractures are discretized uniformly on a separate [n d − 1] dimensional mesh. The weak forms of the flow and mechanics equations are obtained by following the standard variational recipe: multiply the strong form by an appropriately defined, smooth test function set [v p ,v u ,v φ ], integrate the equation over the physical domain, apply integration by parts and reduce the order of differentiation on the primary unknown, apply the divergence theorem, and impose Neumann boundary conditions Tran and Jha (2020). I use element-centered, piecewise constant pressure functions and bilinear functions for the nodal displacement and phase field variable. This choice satisfies the inf-sup condition because I consider a compressible system with finite values of the solid grain compressibility 1/K s and the fluid compressibility c f , which are representative of most subsurface rock formations. Below are the weak forms of the governing equations: Flow Following the standard recipe, I multiply the flow residual with the test function v p and integrate it over the whole domain Z Ω v p ∂ ∂t p m M m (φ) +α (φ)ϵ v −∇· k m µ ∇p m − ψ mf − ψ mw dΩ=0 . (4.38) 102 Integrating by parts and applying the Green’s theorem, we obtain Z Ω v p 1 M m ∂p m ∂t +p m ∂(1/M m ) ∂t +α ∂ϵ v ∂t +ϵ v ∂α ∂t dΩ+ Z Ω ∇v p · k m µ ∇p m dΩ − Z ∂Ω q v p k m µ ∇p m ∂Ω q dΓ − Z Ω (ψ mf +ψ mw )dΩ=0 , (4.39) Above, the time derivative of pressure and the time derivative of strain are the storage and coupling terms. The time derivatives of Biot modulus and Biot coefficient are new coupling terms between flow and geomechanics that arise due to fracturing-induced damage. The fifth term is the fluid mobility term while the sixth term indicates the Neumann boundary condition of the flow equation: k m µ ∇p m ∂Ω q =q w . The seventh term includes the leak-off term (or source term) arising due to flux transfer between matrix and fracture. The weak form of the fluid mass continuity in the fracture subdomain is Z Ω v p 1 M f ∂p f ∂t +p f ∂(1/M f ) ∂t dΩ+ Z Ω ∇v p · k f µ ∇p f dΩ − Z ∂Ω q v p k f µ ∇p f ∂Ω q dΓ − Z Ω (ψ fm +ψ ff +ψ fw )dΩ=0 , (4.40) To derive the weak form, I introduce the finite element approximations: u =N u U, p m =N p P m , and φ = N φ Φ . N u is the displacement shape function matrix, N p is the pressure interpolation function matrix, N φ is the phase field interpolation function matrix, U is the nodal displacement vector,P is the element-centered pressure vector, andΦ is the element-centered phase field variable vector. The test functions are interpolated similarly, i.e., v p =N p V p . The volumetric strain in the flow equation can be expressed in terms of the strain displacement matrix B u : ϵ v =∇· u=1L u u=1L u N u U =1B u U u=[u x ,u y ] T , L u =        ∂ ∂x 0 0 ∂ ∂y ∂ ∂y ∂ ∂x        , 1=[1,1,0] (4.41) 103 whereB u =L u N u . By substitution, I obtain the discretized form of the fluid flow equation: Z Ω N p 1 M m N p dP m dt +N p P m d(1/M m ) dt +α 1B u dU dt +1B u U dα dt dΩ + Z Ω ∇N p k m µ ∇N p P m dΩ − Z ∂Ω q N p q m dΓ − Z Ω (ψ mf +ψ mw )dΩ=0 , (4.42) We define the storativity matrix S m , the coupling matrixQ, and the mobility matrix for rock and fractureA m as follows: S m = Z Ω N p 1 M m N p dΩ Q= Z Ω 1N p α B u dΩ A m = Z Ω ∇N p k m µ ∇N p dΩ= Z Ω B p k m µ B p dΩ (4.43) Weusethefirst-orderbackwardEulermethodfortimeintegration: dP dt ≈ P n+1 − P n δt , dU dt ≈ U n+1 − U n δt , and δt = t n+1 − t n , where n and n+1 denote the current and future time steps, respectively. I rearrange terms such that the terms containing the unknown vectors, P m,n+1 , P f,n+1 , and U n+1 , are separate from other terms that are known at the beginning of t n+1 timestep. This yields the final discrete forms of the flow continuity equations for the matrix and fracture subdomains: (2S m,n+1 − S m,n +δt A m )P m,n+1 + 2Q n+1 − Q n U n+1 = δt f m p +S m,n P m,n +Q n U n (2S f,n+1 − S f,n − δt A f )P f,n+1 =δt f f p +S f,n P f,n (4.44) Here, the right-hand side terms f m p andf f p of equation (4.44) consist of prescribed fluxes and well terms in the matrix and fracture domains respectively. Mechanics Using a suitably smooth test function vector v u , we obtain the strong form of the mechanics equation: Z Ω v u :{∇· [g(φ)σ ′ − α (φ)Ip m ]}dΩ=0 , (4.45) 104 Applying integration by parts, using the divergence theorem, and substituting the finite element interpolation: v u =N u V u ,σ ′ =Dϵ =DLu=DLN u U =DB u U. Defining the stiffness matrix asK = R Ω B u DB u dΩ , I obtain the discretized form of the governing equation: g(φ)KU n+1 − Q n+1 P m,n+1 =f u (4.46) The right-hand side of the weak form of the mechanics problem consists of a traction term, f u = R Γ t N T σ dΓ where Γ t denotes the Neumann boundaries of the domain with the prescribed traction vectorσ . Phase Field Multiplying with a suitable test function v φ and integrating over the domain, I obtain the strong form for the phase field equation: Z Ω v φ − G c l 0 ∇ 2 φ+ G c φ l 0 − (1− φ)η + dΩ=0 (4.47) Applying integration by parts, expanding the first term using Green’s theorem, using the Neumann boundary condition on the phase field variable, and substituting v φ =N φ V φ , we obtain the discrete form of the phase field governing equation: Z Ω −∇ N φ G c l 0 ∇N φ Φ +N φ G c l 0 +η + N φ Φ − N φ η + dΩ=0 (4.48) We define W = Z Ω − B φ G c l 0 B φ dΩ , B φ =L φ N φ =∇N φ , L φ = ∂ ∂x , ∂ ∂y T Z = Z Ω N φ G c l 0 +η + N φ dΩ f φ = Z Ω N φ η + dΩ (4.49) 105 Finally, the discrete matrix-vector form of the phase field equation is (Z +W)Φ =f φ (4.50) 4.4.2. Coupling Strategy and Solution Scheme Our coupled physics computational framework consists of three main components: geomechanics, fluid flow, and fracture mechanics (Figure 4.1). Each component is connected to the other via a two-way coupling link. The primary unknowns are the matrix pressure, fracture pressure, matrix displacement, and the phase field variable. The strength of coupling is dictated by different param- eters: Biot coefficient α , Young’s Modulus E, fluid compressibility c f for flow–mechanics coupling; fracture toughness G c for phase field–mechanics coupling, and joint roughness coefficient (JRC) and joint compressive strength (JCS) for fracture mechanics–rock mechanics coupling. I employ a staggered scheme to solve the fully coupled system by performing a fixed point iteration from a simultaneous flow–mechanics solver to the phase field solver. The linearized matrix-vector system of equations to solve the coupled flow–geomechanics problem in the presence of wells with prescribed injection/production rates becomes:          R uu R up 0 0 R pu R m pp R mf pp R mw pp 0 R fm pp R f pp R fw pp 0 R wm pp R wf pp R w pp                   U P m P f P w          n+1 =          f u δt f m p +Q n U n +S m,n P m,n δt f f p +S f,n P f,n q w          (4.51) where I define R uu =[(1− φ) 2 +κ ]K n , R up =− Q n+1 , R pu =2Q n+1 − Q n R m pp =2S m,n+1 − S m,n +δt A m , R mf pp =− δt T mf , R mw pp =− δt T mw R fm pp =δt T fm , R f pp =2S f,n+1 − S f,n − δt A f , R fw pp =δt T fw R wm pp =T wm , R wf pp =T wf , R w pp =− A w (4.52) 106 At each time step in the time loop, because of the phase-field-dependent coefficients of equa- tion (4.51) and the impact of pressure and displacement on the strain energy (i.e., on η + ), I need to have an inner loop to perform necessary iterations to achieve convergence. The iteration starts with the initialized phase field variable, which is fixed to calculate the poroelastic parameters and the degraded stiffness matrix. Then the mechanical equilibrium and pressure equations are solved concurrently to generate new displacement and pressure fields, which are used to update the history functional η + . This new value of η + is used to solve equation (4.50) and obtain the new phase field for the next iteration in the inner loop. I continue until convergence is achieved. The simulation stops after a specified amount of injection time (Figure 4.2), denoted by t max . Within the inner loop, the flow transmissibility is calculated based on the updated geometric factors of the NNCs and the stress-dependent fracture aperture w f . Following the inner loop’s convergence, the fracture tip-tracking algorithm is activated to determine the new damage extent of the domain, which is followed by the calculation to update NNCs. 4.5. Model Verification and Benchmarking We perform a series of benchmark tests to demonstrate the capabilities of my EDFM-PFM frame- work in modeling fracture propagation under poroelasticity. First, to verify the decoupled phase field solver (i.e., no fluid flow coupling), I use the classic single-edge notched specimen subjected to tensile fracturing (Figure 4.3). Next, to benchmark the integration of phase field solver with poroelasticity, I use the well-known analytical solutions of Detournay Detournay (2004) to verify the dynamic evolution of fluid-driven fracture propagation. 4.5.1. Single-edge Notched Tension Test Wetestmyformulationintheabsenceoffluidflowwithaclassiccaseofasingle-edgenotchedtension test. I consider a square plate with a horizontal crack (notch) that connects the midpoint of the left boundary edge to the center of the domain (Figure 4.3). The nodes on the top and bottom boundary edges have prescribed displacements. The axial displacements of the bottom boundary nodes are constrained to be zero. The axial displacements of the nodes on the top edge are incremented at each time step, providing the applied loading to the system. The increment is δu y = 1× 10 − 4 107 φ( v ) Geomechanics Fracture mechanics Fluid flow ψ ∂σ v ∂t k f ( δ h , σ n ) M ( ϕ) α( ϕ) σ ( p) η + ( p) η + ( u) g ( ϕ) α( ϕ) σ , Figure 4.1: Flow and geomechanics are two-way coupled through strain-dependent porosity ϕ (ϵ v ) and the stress-dependent mass transfer function ψ in the fluid mass balance equation and the pressure-dependent effective stress σ ′ (p) in the linear momentum balance equation. The flow and fracture mechanics modules are coupled through the phase-field-dependent geomechanical param- eters M(φ) and α (φ) and the pressure-dependent strain energy functional η + (p). In addition, the fracture mechanics module delivers stress-dependent permeability k f (δ h ,σ ′ n ) to the flow module based on shear slip δ h and normal dilation, which is dependent on the effective normal stress σ ′ n . The geomechanics module is coupled to fracture mechanics by providing the displacement, stress, and strain to the Bandis model for slip activation, slip magnitude calculation, and fracture aperture update, as well as evaluating the strain energy functional η + (u). The fracture mechanics module updates the Biot coefficient α (φ) and the stress degradation function g(φ) in the geomechanics module. 108 Model Input and Initialization Update Time Loop: while No Phase Field Solver Flow-Geomechanics Solver Inner Loop: while Yes Post-processing No Yes t<t max t < t max ? K f , φ m , φ f , w f , x new tip Err p m , p f , u, ϕ > tol? Err p m , p f , u, ϕ > tol Figure 4.2: Process diagram of the coupled flow-geomechanics-fracture propagation framework. In the initialization step, the initial local strain history functional is specified as well as the initial pressure/displacement and permeability field. Inside the time loop, there is an inner loop to iterate through the solutions of the two solvers, flow-geomechanics solver and the phase field solver, to a specified error tolerance tol. Poroelastic parameters are updated for use in the next time step based on the converged solution of the phase field in the inner loop. Parameter Magnitude Unit Young’s Modulus E 210 GPa Poisson’s Ratio ν 0.3 m Critical energy release rate G c 2700 N/m Length scale parameter l 0 12.5× 10 − 6 m Table 4.1: Tension Test Parameters. mm up to u y = 4.5× 10 − 3 mm and δu y = 1× 10 − 6 mm up to the failure of the specimen. The external pull provides the tensile elastic energy to the domain for facilitating fracture propagation. To better capture the dynamic propagation process, the time step size is reduced as the onset time of propagation approaches by reducing the boundary displacement increment by 5000 times. Other parameters are specified in Table 4.1. The simulations are conducted under the plane strain condition. The resulting reaction forces at the bottom boundary are integrated periodically. They are output together with the current axial displacement value of the top boundary, to generate the force–displacement curve. The resulting crack patterns at different stages of boundary displacement are evaluated, and the force acting at the fixed bottom boundary is monitored and compared with published solutions Miehe et al. (2010); Natarajan et al. (2019). 109 Figure 4.3: Left: A square plate with a straight edge crack subjected to vertical pull on the top boundary to simulate a pure mode-I fracture. They-displacement on the top boundary is prescribed and increases monotonically with time. Right: Resultant force at the bottom boundary increases with the prescribed displacement until the onset of failure, which is marked with a sharp drop in the force. The proposed model agree with the literature. 4.5.2. Fluid-filled Fracture Propagation As a benchmark problem, I consider the fluid-driven propagation of a single horizontal fracture at the center of the domain. All sides of the domain are fixed roller boundaries with zero fluid pressure. To initiate fracture propagation, I activate fluid injection at a point at the center of the fracture. The key parameters are listed in the Table 4.2. Before I present my solution of the problem, it is worthwhile to discuss the difficulty in solving this problem, analytically or numerically. Analytical solution of hydraulic fracturing is challenging to establish because of the associated nonlinearity of the equations governing the flow of fluid in the fracture and the impact of the elastic response of the fracture on the domain. The nonlinearity originates from the dependence of the fracture permeability on the fracture aperture. Furthermore, the process is time-dependent; the fracture propagates with an unknown, time-dependent speed. Nonlocality and nonlinearity combine to result in a complex solution structure that involves processes at the micro-scale, near the tip of the fracture. These small-scale processes control the global response of the fluid-driven fracture. Numerical simulation of hydraulic fracturing is challenging because it involves the coupling between (i) mechanical deformation, (ii) fluid flow inside the fracture and the reservoir, and (iii) fracture propagation which changes the shape and size of the fracture subdomain (Γ f ) over which 110 Parameter Magnitude Unit Young’s Modulus E 210 GPa Poisson’s Ratio ν 0.3 m Critical energy release rate G c 2700 N/m Length scale parameter l 0 12.5× 10 − 6 m Table 4.2: 2D Fluid-filled fracture propagation. p = 0 p = 0 p = 0 p = 0 40 m ( ) p f ( x, t ) w f ( x, t ) q w σ σ Figure 4.4: Model setup to test the proposed framework on the fluid-filled fracture propagation problem. A hydraulic fracture propagating in an impermeable elastic medium along the direction orthogonal to the orientation of minimum compressive far-field stress σ . Fracture propagation is drivenbyinjectionofanincompressible, Newtonianfluidataconstantvolumetricrate q w Detournay (2004). fracture-relatedquantitiesarecomputed. Theframeworkproposedherecanaddressthesechallenges smoothly as evident from the discussion below. We observe that the fluid pressure is nearly constant within the crack due to the increased perme- ability along the main direction of the crack. Because the time step size is small compared with the characteristic time of pressure diffusion, the fluid pressure outside the crack is much lower than the pressure within the regularized crack. The contours of the elevated pressure area is approximately elliptical, which can be explained by the well-known “leak-off” phenomenon describing the fluid leaking from the hydraulic fracture into the surrounding rock matrix. The pressure rise caused by the injection of fluid, is offset by the pressure drop caused by fracture extension and fluid leak-off. As observed, fluid pressure inside the fracture elevates rapidly before reaching a peak value where 111 Phase Field X -displac ement (meter) Pressure (MPa) Y-displacement (meter) t = 3 sec t = 20 sec t = 50 sec Figure 4.5: From top to bottom: Spatial distributions of the phase field variable, fluid pressure, and horizontal and vertical displacements at three selected time steps. We observe a large pressure difference with low pressures (in blue color) in front of the fracture tips. Displacements are zero at the tips of the fracture, and the largest displacement is observed in the middle, a typical behavior of crack opening. the initial crack begins to propagate, causing a slight reduction of the fluid pressure. The frac- ture thickening effect observed during the simulation is typical of phase field-based approaches for simulating hydraulic fracturing. The difference in fracture length and pressure between my model and Detournay’s analytical model is due to the boundary effect in my simulation. Constant zero pressure at the boundaries makes the pressure gradient higher and fracture propagation faster. In Figure 4.6, I compare my numerical solution with a well-known analytical solution for fluid injection-driven fracture propagation Detournay (2004). Fluid pressure increases within the crack before reaching a peak point where the initially defined crack begins to propagate, resulting into a decrease in the fluid pressure. The solution agrees well with Detournay’s analytical approach for fluid-driven fracture growth in the viscosity-dominated regime, where the energy expended in the 112 Sim. Frac. Length (m) time (sec) time (sec) Pressure in Stim. Frac. (Pa) a a Figure 4.6: Evolution of the fracture pressure and fracture length with time. The results are compared with the analytical solution Detournay (2004). The slight mismatch stems from the difference in boundary condition which affects the solution when the fracture approaches the edge of the domain. creation of new fracture surfaces in the rock is relatively small compared to the energy dissipated in viscous flow. In the toughness-dominated regime, the viscous dissipation is small compared to the energy dissipated at the crack tip. 4.6. Stress Shadow and Fracture Interaction Building on the foundation of a single fracture model, I construct a mechanistic model to investigate the impact of stress shadow and fracture interaction on fracture propagation. From the undrained deformation theory Wang (2000), an external load applied to a fluid-filled porous material raises the fluid pressure inside the pores by an amount that is linearly proportional to the volumetric total stress applied with the Skempton’s coefficient as the proportionality factor. The excess pressure progressively dissipates through the medium until equilibrium is achieved. In a tight formation with low permeability, the dissipation process is slow enough that the excess pressure is main- tained throughout the hydraulic stimulation phase where the propagating hydraulic fracture creates the volumetric stress perturbation. These poroelastic responses can be captured by downhole and surface pressure gauges. Past research has shown that it is possible to utilize this induced pres- sure signal for the purposes of fracture monitoring and diagnostics (Kampfer and Dawson, 2016b; Manjunath et al., 2023) to decipher fracture morphology parameters such as fracture orientation, length, and spacing from other pre-existing fractures. Based on the framework proposed here, I build a mechanistic model to test this hypothesis. In my setup, there is a pre-existing fracture near the top boundary, acting as the monitoring point for pressure induced by the growing fracture near 113 Parameter Magnitude Unit Dimensions L x ,L y 40 m Number of grid cells N x , N y 120 Matrix porosity ϕ m 0.15 Fracture porosity ϕ f 0.001 GPa Matrix permeability k m 5 mD Fluid density ρ f 1000 kg/m 3 Fluid viscosity µ 1 cP Fluid compressibility c f 0.455 1/GPa Fluid diffusivity D f 2× 10 − 7 m 2 /s Young’s modulus E 16 GPa Poisson’s ratio ν 0.2 m Solid grain density ρ s 2.65 g/cc Biot coefficient matrix α (φ=0) 0.79 Biot coefficient fracture α (φ=1) 1 Critical energy release rate G c 80 N/m Critical cohesive tensile stress σ c 0.45 MPa Length scale parameter l 0 0.33 m Table 4.3: Model parameters for investigating stress shadow and fracture interaction Parameter Low· Medium· High Unit Approach angle β 0· 60· 90 degree Spacing between the two fractures 8· 16· 22 m Matrix Young’s Modulus 8· 16· 22 GPa Table 4.4: Sensitivity parameter values. the bottom boundary (Figure 4.7). During stimulation, pore pressure inside the bottom fracture elevates, prompting its propagation. I will investigate the impact of different fracture morphology to evaluate whether or not the induced pressure signal is sensitive to those parameters. The key parameters of this simulation are listed in Table 4.3. 4.6.1. Effect of Fracture Approach Angle In the subsurface, the direction of minimum principal stress may vary in space, influencing the path of hydraulic fracture propagation. It is critical to quantify the impact of stress shadow exerted by one fracture on the other, especially when multiple fractures are created side-by-side along a well or when the wells with fractures are closely spaced in the field. One of the parameters that can be tracked is the angle of approach of the growing or stimulated fracture to an existing observation 114 Phase Field Pressure (MPa) t = 30 sec t = 60 sec t = 120 sec Obs. Frac. Stim. Frac. Figure 4.7: Propagation of a hydraulically stimulated fracture (Stim. Frac.) monitored via a static observation fracture (Obs. Frac.). Spatial distribution of the phase field variable and the fluid pressure at three selected time steps (along columns) are shown. The thin zig-zag line inside the Stim. Frac. in the bottom row panels show the fracture propagation trajectory determined by my fracture tracking algorithm in section 4.3.5. fracture. The two fractures can be located in the same well or in two different wells. We evaluate how the fracturing-induced pressure signal in the observation fracture varies with the stimulated fracture’s approach angle. In Figure 4.8, I tested different approach angle values with β ranging from 0 to 90 degrees. We find that the damage zone around the stimulated fracture grows in width perpendicular to the fracture orientation due to the coupling between mechanical degradation and fluid leak-off at the fracture interfaces. Another interesting result is the change of direction of fracture propagation as the stimulated fracture approaches the domain boundary. Naturally, the fracture bends to find ‘the path of the least resistance’ and more space to grow. The left panel of Figure 4.9 shows the variation of pressure inside the stimulated fracture during the fracturing process. The fracture initiation pressure is independent of the approach angle. As the fracture starts to propagate, pore pressure declines from the peak initiation pressure, with fluctua- tions corresponding to individual failure episodes. The average decline rate is different for different approach angles. As the approach angle increases (fracture turns from horizontal to vertical), the pressure drops faster for a given growth in the length because the stress shadow effect from the 115 Time Increasing approach angle t = 3 sec t = 20 sec t = 10 sec β = 0 β = 30 β = 60 β= 90 t = 30 sec t = 40 sec β Obs. Frac. Stim. Frac. Figure 4.8: Snapshots of the spatial distribution of the phase field at five selected time steps for different approach angles β . The rows are for different approach angles which change the fracture orientation from horizontal to vertical in the domain. The columns from left to right indicate the evolution in time. 116 Stim. Frac. Length (m) Normalized Time Normalized Pressure in Obs. Frac. degree Pressure in Stim. Frac. (MPa) (a) (b) Figure 4.9: Effect of the approach angle β of the stimulated fracture to the observation fracture. (a) As the stimulated fracture’s orientation with respect to the observation fracture changes from parallel (β =0) to perpendicular (β =90), the stimulated fracture pressure gets smaller on average. For any given value of the angle, the pressure first drops and then rises as the fracture grows to sense boundaries in the domain. (b) Observation fracture pressure (normalized) versus the normalized time for different approach angles. observation fracture is weaker. For the 90 degree approach angle case (stimulated fracture ori- ented vertically to the observation fracture), the simulation duration is shorter because the fracture reaches the domain boundary sooner. The right panel of Figure 4.9 shows the temporal evolution of pore pressure in the observation fracture, normalized by the concurrent pressure at the stimu- lated fracture. The impact of fracture angle is apparent from the figure; a smaller β results in a larger observation pressure. This is desirable in the field because the signal-to-noise ratio improves, leading to robust estimates of the stimulated fracture’s length and orientation from interpreting the observation fracture’s pressure. When the stimulated fracture grows perpendicular to the observa- tion fracture, it produces almost no response in the observation fracture, which implies that the stimulated fracture properties cannot be estimated from the observation fracture’s data. The isotropic part of the deformation (volumetric expansion or contraction) is quantified by ϵ v which is controlled by the effective volumetric stress σ ′ v , which is related to the total volumetric stress σ v = tr(σ )/n d as σ ′ v =σ v +αp m . The deviatoric part of deformation is controlled by the von Mises stress σ VM = √ 3J 2 , where J 2 is the second invariant of the deviatoric stress tensor (σ − σ v I) withI is the identity tensor Tran and Jha (2020). To characterize the state of stress in the domain, three different stress quantities are chosen: effective volumetric stress σ ′ v for isotropic deformation and the von Mises stress σ VM for shear-induced failure. 117 Von Mises Stress (Pa) Eective Volumetric Stress (Pa) Von Mises Stress (Pa) Eective Volumetric Stress (Pa) A B Stress trajectory of point A Obs. Frac. Stim. Frac. Stress trajectory of point B β = 0 ◦ β Figure 4.10: Effect of the fracture approach angle β on stress trajectories in the invariant space of von Mises stress versus effective volumetric stress. Stress trajectory at a point also depends on the point’s position relative to the stimulated fracture. Left: Point A is at the center of the domain and falls on the path of the β =90 ◦ stim. frac. (stimulated fracture becoming perpendicular to the obs. frac.). Since stimulated fractures are predominantly tensile-mode failures, the effective volumetric stress is tensile (positive) for β = 90 ◦ and remains tensile even after the tip moves past point A. For other values of β , the effective volumetric stress is compressive (negative), which is the origin of the stress shadow effect. Right: Point B falls closest to the path of the β = 30 ◦ stim. frac. As the fracture tip approaches point B, the effective volumetric stress becomes tensile. Further growth of the stimulated fracture causes compression at point B (turning of the trajectory) due to the stress shadow effect. 118 Time Increasing approach angle t = 3 sec t = 20 sec t = 10 sec β = 0 β = 30 β = 60 β = 90 t = 30 sec t = 40 sec σ ( MP a ) Figure 4.11: Spatial distribution of the total volumetric stress at five selected time steps for different values of the fracture approach angle β . 119 Increasing approach angle t = 3 sec t = 20 sec t = 10 sec β = 0 β = 30 β = 60 β = 90 t = 30 sec t = 40 sec σ ( ) Time MPa Obs. Frac. Stim. Frac. β Figure 4.12: Spatial distribution of the von Mises stress at five selected time steps for different values of the approach angle β . Deviatoric stresses are affected by the approach angle. 120 4.6.2. Effect of Fracture Spacing Horizontal wells drilled in low-permeability rocks are often fractured at multiple positions along the horizontal section of the wellbore to enhance hydraulic conductivity of the surrounding rock and establish connectivity between the well and the surrounding rock. The idea is to create multiple quasi-vertical fractures that branch out to efficiently drain the stimulated region. This is called multi-stage hydraulic fracturing Manjunath et al. (2023), where ‘stage’ refers to a length interval along the wellbore subjected to high-pressure injection to initiate fracturing. The distance between neighboring fracture clusters or stages in the stimulated well is an important design parameter of a hydraulic fracturing job. Here, I investigate the sensitivity of the fracturing-induced pressure signal, measured at the observation fracture, to the spacing between the observation and stimulated fractures. I perform simulations for three values of the fracture spacing parameter, as shown in the first column of Figure 4.13. The pressure at the observation point increases with time as the stimu- lated fracture propagates, which is expected from (1) growing compression in the rock surrounding the stimulated fracture, and (2) diffusion of the injection pressure through the permeable medium between the two fractures. Both mechanisms are strongest at the smallest value of fracture spacing, which explains the highest observation pressure for the lowest spacing value. Figure 4.14 illustrates the impact of fracture spacing on the induced pressure signal observed in the observation fracture. An increase in the fracture spacing by a factor of 3 (from 8 meters to 22 meters) results in more than doubling of the induced pressure signal at the end of the simulation when the normalized time reaches 1. By the end of the simulation, the observation pressure in the 8-meter spacing case is equal to the injection pressure in the stimulating fracture, suggesting a complete interference or “short-circuit" between the two fractures, which may be undesirable in the field. 4.6.3. Effect of mechanical heterogeneity Subsurface geomechanical heterogeneity plays a central role in determining the characteristics of the hydraulicfracturenetworkaswellastheextentofstressshadoweffect. Toinvestigatethis, Iconsider a scenario where the domain is made up of two layers of different mechanical properties (Young’s Modulus and Biot’s Coefficient) and the two fractures (stimulated fracture and observation fracture) 121 Increasing fracture spacing Time 4.5 3.5 2.5 1.5 0.5 0 MPa 22 m 16 m 8 m Obs. Frac. Stim. Frac. Figure 4.13: Evolution of pore pressure at five selected time steps (along columns) for three values of the fracture spacing parameter (along rows). Normalized Time Normalized Pressure in Obs. Frac. Pressure in Stim. Frac. (MPa) Stim. Frac. Length (m) (a) (b) Figure 4.14: Observation fracture pressure (normalized) versus the normalized time for three dif- ferent fracture spacing values. Fracture spacing is the distance between the stimulated fracture and the observation fracture. At closer spacing, the stimulated fracture has lower pressure and shorter length because of the pressurization of the observation fracture. 122 Time (sec) Normalized time Obs. Frac. Stim. Frac. Normalized Pressure in Obs. Frac. E /2, α E , α 2E , α E , α/2 E /2, α E , α 2E , α E , α/2 Stim. Frac. Length (a) (b) Figure 4.15: Effect of change in mechanical properties (Young’s modulus E and Biot’s coefficient α ) of the rock surrounding the stimulated fracture region (colored box on the left). (a) Observation fracture pressure, and (b) length of the stimulated fracture. A 50% reduction in α causes almost 54% reduction in the observation fracture pressure are located in separate layers. I analyze the effect of such heterogeneity on the stimulated fracture length and the induced pressure signal in the observation fracture (Figure 4.15). A lower value of Biot’s coefficient weakens the fluid-to-solid coupling (recall the effective stress equation (4.21)), thus diminishing the induced pressure signal. At the same time, a lower value of Biot’s coefficient implies more consolidated rock and, therefore, faster fracture propagation and a longer fracture after a certain time duration. Also, rocks with higher tensile stiffness (larger Young’s modulus) show a tendency to facilitate faster fracture propagation. 4.7. Field Case Study We conduct a field case study of pressure-based monitoring of hydraulic fracturing. I consider a horizontal section of a reservoir with a triplet of horizontal wells, named W1, W2, and W3; see Figure 4.16. The sequence of injection operation is as follows. Well W1 and W3 were hydraulically stimulated to generate a near-orthogonal fracture set at the two wells, with the fracture orientation estimatedusingmicroseismicdataanalysisinthefield. Afterthat, thecenterwell, W2, isstimulated to generate a few near-orthogonal hydraulic fracture sets. Stimulation of W2 is monitored by recording pressures at W1 and W3. The pre-existing fractures in the monitoring wells contribute to the observation pressure signal generated by W2’s fracturing. The reservoir is fully saturated with water, and the injected fluid to create the fractures is ‘slickwater’, which has physical properties (viscosity and density) similar to the pore water. We build a model to simulate the fracturing of well W2 while monitoring pressure at W1 and W3. 123 Parameter Magnitude Unit Dimension L x 914.4 meter Dimension L y 304.8 meter Number of grid cells N x 300 Number of grid cells N y 100 Matrix porosity ϕ m 0.1 Young modulus E 41.47 MPa Critical energy release rate G c 17.3 kN/m Critical cohesive tensile stress σ c 3.5 MPa Length scale parameter l 0 3.1 meter Table 4.5: Field case study input parameters The model domain boundaries are set at constant pressure with roller conditions (fixed displace- ments). The origination points of the propped fractures along the lateral sections of the wells are assumed to be constrained by the locations of the perforation holes and hydraulic isolation packers in the wellbore. We assume single-phase flow inside the fluid-driven hydraulic fracture embedded in a deformable poroelastic solid. The parameters of the simulation of this field trial are listed in Table 4.5, with the rest of parameters are similar to Table 4.3. Opening of a hydraulic fracture in the stimulated well can be ‘seen’ as a pressure jump in an observation well that is sufficiently close in space and is properly instrumented with either downhole or wellhead pressure gauges. The pressure signal has two parts: (1) poroelastic part related to the undrained deformation of the rock in the vicinity of the propagating fracture and (2) hydraulic part due to movement of the fluid from the stimulated well to the observation well, also known as the ‘frac hit’ (Lawal et al., 2013; Manjunath et al., 2023). The time scales and the magnitudes of the two parts depend on the hydraulic properties (permeability, porosity, fluid compressibility, and fluid viscosity) and poroelastic properties (Skempton’s coefficient, drained bulk modulus and Poisson’s ratio, and Biot coefficient). Poroelastic signal travels faster than the hydraulic signal. In low-permeability rock considered here, the poroelastic signal is larger in magnitude than the hydraulic signal. Figure 4.17 shows the measured data in the field. Injection pressure at well W2 and observation pressures at wells W1 and W3 are plotted along with the fracturing fluid (water) injection rate in 124 p = const p = const p = const Monitor W1 Monitor W3 Stimulation W2 stage 2 p = const stage 1 Figure 4.16: Schematic of the field case study. W1, W2, and W3 are three horizontal wells in a horizontal section of a reservoir rock. Fractures are shown by oval shapes placed on the wells. Well W2 is hydraulically stimulated to create fracture stages 1 and 2 while pressure is monitored at wells W1 and W3. W1 and W3 contain one hydraulic fracture each from a prior stimulation job. W2. The pressure in well W2 increases significantly due to the pumping of the fracturing liquid. This causes the observed pressure at W3 to also increase, but gradually and with some delay. This is diagnostic of a poroelastic signal in this low permeability rock and rules out ‘frac hit’ or direct hydraulic communication between W2 and W3. Well W1 exhibits a more pronounced induced pressure signal (≈ 4 MPa) compared to well W3 in response to the fracturing at well W2 (≈ 0.7 MPa). Because of the uncertainties in the reservoir rock properties and fracturing job parameters, the scope of my study is limited to simulating the physical mechanisms of induced poroelastic signals and its qualitative evaluation. Quantitative ‘history matching’ of the simulation model to field measurements is outside the scope of this study. In the numerical model, the orientation and length of existing fractures are based on a prior analysis of microseismic data in the field. Fixed roller boundary and constant fluid pressure are prescribed on all boundaries as per Figure 4.16. The boundary pressure is prescribed to be 3500 psia = 24.1 MPa, which is approximately equal to the initial flowing bottom hole pressure of the stimulation well, to achieve initial equilibrium. Fluid is injected into well W2 to create stage-1 and stage-2 hydraulic fractures, sequentially, while pressure is monitored with time at observation fractures in wells W1 and W3. 125 9.65 13.10 16.55 20.00 23.44 26.89 30.34 Pressure (MPa) Injection Rate (cubic meter/min) 0 2.38 4.77 7.15 9.54 11.92 14.30 Elapsed Time (minute) 0 15 30 45 60 75 W1 Obs. Pres. W2 Inj. Pres. W2 Inj. Rate W3 Obs. Pres. Stage 1 Stage 2 Figure 4.17: Field case study data: Pressure and rate data during hydraulic fracturing of well W2 and its monitoring at wells W1 and W3. Two stages of fracturing are shown with similar injection rates for both stages (approximately 11.13 cubic meter per minute) and similar pressures (29.65 MPa). Stage 1 injection lasts about 30 minutes while stage 2 injection lasts about 20 minutes. 1 min Time 20 min 50 min Phase Field Pressure 70 min MPa W1 Obs. Frac. W3 Obs. Frac. stage-1 W2 Stim. Frac. Stage-1 fracturing Stage-2 fracturing stage-2 Figure 4.18: Field case study: Spatial distribution of fracture pressure and the phase field variable at four selected time steps during stage-1 and stage-2 fracturing at well W2. Left side panels show that injection-induced overpressure in stage-1 (red color) causes pressure diffusion and stress shadow that affects stage-2’s trajectory and length. Right side panels show that the stage-1 fracture grows in a direction that is predominantly parallel to the observation fractures. Stage-2 fracture senses the presence of relatively lower compression inside stage-1 and higher compression in the region between the two stages, which dictates stage-2’s growth in the orthogonal direction and the intersection of the two fractures. 126 1 minute Time 20 minutes 50 minutes 70 minutes Von Mises Stress (MPa) Effective Volumetric Stress (Mpa) σ ′ ( MPa ) σ ( MPa ) Figure 4.19: Field case study: Spatial distribution of effective volumetric stress (left column) and Von Mises stress (right column) at four selected time steps during stage-1 and stage-2 fracturing at well W2. The time steps are the same as the figure above. 127 Figure 4.18 depicts pressure evolution and fracture propagation. The first event is an increase in pressure at well W2’s stage-1 location, which extends the fracture in a direction approximately perpendicular to the wellbore direction; see Figure 4.16. This result agrees with the fracture ori- entation estimated in the field based on the prevailing state of stress (principal stress magnitudes and directions). The pore pressure field shows growing interference towards the location of stage-2 as stage-1 fracture propagates. Due to the opening of stage-1 fracture, the rock above and below stage-1 is under higher compression relative to the far-field stress. Also, the stress inside the stage-1 fracture is least compressive (Figure 4.19). These two aspects of the new stress state cast a stress shadow on the growth of stage-2 fracture, when stage-2 pumping begins. The stage-2 fracture starts to propagate in a direction orthogonal to stage-1 because of the stress shadow effect. The two fractures coalesce, and the stage-2 fracture continues growing towards the north boundary of the domain. Observation wells W1 and W3 capture the entire dynamics in terms of the pressures measured at the respective observation fractures, as shown in Figure 4.20. The poroelastic response at W1 is consistent at ≈ 0.2 MPa for both stage of fracturing in W2. On the other hand, the response at W3 to the two stages varies from 0.15 MPa (stage 1) to 0.45 MPa (stage 2). Overall, the simulated responses are relatively lower than what the field data suggests. As proven in the previous section, the magnitude of poroelastic response is a function of several variables most of which are uncertain, for example, the mechanical properties of the rock (Young’s modulus, Poisson’s ratio, Biot coefficient), the approach angle of the stimulated fracture, and the fracture spacing. Also, the merging of two stimulated fractures complicates the poroelastic responses at the two observation wells, W1 and W3. Due to a large number of uncertainties in terms of hydraulic, petrophysical, and mechanical heterogeneity of the reservoir, the quantitative history matching of this complex stimulation sequence is not part of the scope of this paper. 4.8. Conclusions Dynamic fracture propagation is a complex physical process involving many factors such as the evolving state of stress, flow–geomechanical coupling strength, fluid leak-off, and mechanical degra- 128 Injection pressure (MPa) Elapsed Time (minute) Observed Pressure (MPa) W2 stage-1 W2 stage-2 W1 Obs. Pres. W3 Obs. Pres. Figure 4.20: Field case study simulation: Observed poroelastic signal at wells W1 and W3 when well W2 undergoes stage-1 and stage-2 hydraulic fracturing stimulation. W2 stage-1 and W2 stage-2 curves are injection pressures monitored at two separate locations corresponding to the stage-1 and stage-2 locations; see Figure 4.16. The peak in an injection pressure curve corresponds to tensile failure of the rock. Post-peak evolution corresponds to the fracture propagation period for the respective fracture. dation of the rock’s structural frame. I proposed a mixed EDFM-PFM (embedded discrete fracture method–phase field method) within the coupled flow-geomechanics framework to model fracture in- teraction and fracturing-induced poroelastic signals during hydraulic fracturing of saturated porous media. Theproposedmodelinvestigatesthesynergisticcouplingbetweenporoelasticityandfracture dynamics in fluid-saturated porous media during the events of fracture propagation. The sharp dis- continuities at the fracture interfaces are regularized by a diffused phase field approximation within a continuum formulation, allowing it to flexibly model fracture nucleation, non-singular fracture tips, fracture coalescence, and branching. I showed that heterogeneity in poroelastic properties can affect the fracture propagation solution significantly. By analyzing the interaction between a grow- ing fracture and a static fracture, I demonstrated that deformation-induced poroelastic pressure signal can be diagnostic of fracture morphology and reservoir heterogeneity. 129 CHAPTER 5 Key Conclusions In the first component of my thesis, I conducted an in-depth investigation into the interplay be- tween poroelastic deformation and the spreading and mixing of solutes during viscous fingering, employing high-resolution numerical simulations. my motivation for this research was rooted in the distinct mathematical characteristics of the governing equations, which led us to propose a numer- ical framework involving a unified flow-geomechanics solver coupled sequentially with a transport solver, utilizing the velocity field and a concentration-dependent mobility model as connectors. I identified and analyzed the mechanisms underlying the coupling between transport and deforma- tion, examining their sensitivity to various poromechanical parameters. A sensitivity analysis was performed to gauge the extent of poroelastic influence on key macroscopic transport metrics, in- cluding breakthrough time, breakthrough curve, and mixing degree. The results demonstrate the impact of poroelastic processes on solute transport. Neglecting this coupling in modeling can lead to inaccurate predictions of breakthrough times, solute recovery rates, and mixing levels, thereby compromising the accuracy of assessments in groundwater remediation or enhanced oil recovery operations. In the second part of my thesis, I introduced an innovative multiphysics modeling framework to explore the complex interactions between flow, geomechanics, transport, and fracture dynamics in naturally fractured aquifers. My work demonstrates that integrating transport-geomechanical cou- pling in field-scale modeling of flow through fractured rocks can substantially enhance predictions of critical transport parameters, such as tracer breakthrough time, average and maximum tracer con- centrations at a control plane, and the evolution of mixing over time. Furthermore, I leveraged the Embedded Discrete Fracture Modeling (EDFM) approach to simulate fractured porous media and employed an improved Bandis model to simulate fracture mechanics. My EDFM-based framework allowedformodelingmultiplerandomlyorientedintersectingfracturesatalowercomputationalcost comparedtoothermethodsrequiringaconformingmesh. Iquantifiedtheinfluenceofgeomechanical 130 coupling on key transport metrics and identified the relative impact of various coupling parameters on tracer breakthrough and mixing. The characteristics of tracer retrieval curves, including peak concentration, were influenced by both fracture morphology and the strength of geomechanical coupling. The interplay between flow-transport coupling, involving viscous fingering and fracture dynamics, altered the stress state and fracture stability within the domain, a critical consideration in formations close to failure. I quantified the impact of the strength of flow-geomechanics coupling on the evolution of fracture permeability. My findings suggest that coupling flow, transport, and geomechanics is essential for accurately modeling scenarios involving active tracers (which mod- ify flow viscosity or density) during reservoir surveillance, contaminant remediation, or fracture characterization studies in deformable rocks. In the final part of my thesis, I introduced a mixed EDFM-PFM (embedded discrete fracture method-phase field method) within the coupled flow-geomechanics framework to model fracture interaction and poroelastic signals induced by fracturing in saturated porous media. This model explores the intricate interplay between poroelasticity and fracture dynamics during fracture prop- agation events. It employs a diffused phase field approximation to regularize sharp discontinuities at fracture interfaces within a continuum formulation, allowing for flexible modeling of fracture phenomena such as nucleation, non-singular fracture tips, coalescence, and branching. My research revealed that heterogeneity in poroelastic properties can significantly affect fracture propagation solutions. Additionally, through an analysis of the interaction between a growing fracture and a static fracture, I demonstrated that deformation-induced poroelastic pressure signals can serve as diagnostic indicators of fracture morphology and reservoir heterogeneity. 131 BIBLIOGRAPHY Y Abousleiman, AH-D Cheng, L Cui, E Detournay, and J-C Roegiers. Mandel’s problem revisited. Geotechnique, 46(2):187–195, 1996. A. N. Alshawabkeh, N. Rahbar, and T. Sheahan. A model for contaminant mass flux in capped sediment under consolidation. J. Contam. Hydrol., 78:147–165, 2005. F. Armero and J. C. Simo. A new unconditionally stable fractional step method for non-linear cou- pled thermomechanical problems. International Journal for Numerical Methods in Engineering, 35(4):737–766, 1992. https://doi.org/10.1002/nme.1620350408. Miguel Arriaga and Haim Waisman. Stability analysis of the phase-field method for fracture with a general degradation function and plasticity induced crack generation. Mechanics of Materials, 116:33–48, 2018. Pooyan Asadollahi and Fulvio Tonon. Constitutive model for rock fractures: Revisiting barton’s empirical model. Engineering Geology, 113(1-4):11–32, 2010. K. Aziz and A. Settari. Petroleum Reservoir Simulation. Elsevier, London, 1979. Alireza Baghbanan and Lanru Jing. Stress effects on permeability in a fractured rock mass with correlated fracture length and aperture. International journal of rock mechanics and mining sciences, 45(8):1320–1334, 2008. M. Bai and D. Elsworth. Coupled Processes in Subsurface Deformation, Flow and Transport. ASCE Press, 2000a. M. Bai and D. Elsworth. Coupled Processes in Subsurface Deformation, Flow, and Transport. Amer. Soc. Civil Eng., Reston,VA, 2000b. Mao Bai, Derek Elsworth, and Jean-Claude Roegiers. Multiporosity/multipermeability approach to the simulation of naturally fractured reservoirs. Water Resources Research, 29(6):1621–1633, 1993. https://doi.org/10.1029/92WR02746. SC Bandis, AC Lumsden, and NR Barton. Fundamentals of rock joint deformation. In International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts, volume 20, pages 249–268. Pergamon, 1983. N.BartonandV.Choubey. Theshearstrengthofrockjointsintheoryandpractice. Rock Mechanics, 10:1–54, 1977. Nick Barton. Modelling rock joint behavior from in situ block tests: implications for nuclear waste repository design, volume 308. Office of Nuclear Waste Isolation, Battelle Project Management 132 Division, 1982. Nick Barton, S Bandis, and K Bakhtar. Strength, deformation and conductivity coupling of rock joints. In International journal of rock mechanics and mining sciences & geomechanics abstracts, volume 22, pages 121–140. Elsevier, 1985. B. Berkowitz, A. Cortis, M. Dentz, and H. Scher. Modeling non-Fickian transport in geological formations as a continuous time random walk. Rev. Geophys., 44:RG2003, 2006. A. Bonazzi, B. Jha, and F. P. J. de Barros. Transport analysis in deformable porous media through integraltransforms. Int. J. Numer. Anal. Methods Geomech., pages1–18, 2020. 10.1002/nag.3150. A. Bonazzi, M. Morvillo, J. Im, B. Jha, and F. P. J. de Barros. Relative impacts of permeability heterogeneity and viscosity contrast on solute mixing. Phys. Rev. Fluids, 6:064501, Jun 2021. 10 .1103/PhysRevFluids.6.064501. URL https://link.aps.org/doi/10.1103/PhysRevFluids.6 .064501. Michael J Borden, Clemens V Verhoosel, Michael A Scott, Thomas JR Hughes, and Chad M Landis. A phase-field description of dynamic brittle fracture. Computer Methods in Applied Mechanics and Engineering, 217:77–95, 2012. Ronaldo I Borja. Plasticity: modeling & computation. Springer Science & Business Media, 2013. Blaise Bourdin, Gilles A Francfort, and Jean-Jacques Marigo. The variational approach to fracture. Journal of elasticity, 91(1-3):5–148, 2008. F. Brezzi. On the existence, uniqueness and approximation of saddle-point problems arising from lagrangian multipliers. RAIRO Analyse numerique, 8:129–151, 1974. A. Bubshait and B. Jha. Coupled poromechanics-damage mechanics modeling of fracturing during injection in brittle rocks. Int. J. Numer. Methods Eng., pages 1–21, 2019. 10.1002/nme.6208. M. Budhu. Earth fissure formation from the mechanics of groundwater pumping. Int. J. Geomech., 11:1–1, 2011. Tuanny Cajuhi, Lorenzo Sanavia, and Laura De Lorenzis. Phase-field modeling of fracture in variably saturated porous media. Computational Mechanics, 61(3):299–318, 2018. N. Castelletto, J. A. White, and H. A. Tchelepi. Accuracy and convergence properties of the fixed-stress iterative solution of two-way coupled poromechanics. Int. J. Numer. Anal. Methods Geomech., 39:1593–1618, 2015a. 10.1002/nag.2400. Nicola Castelletto, Joshua A White, and HA Tchelepi. Accuracy and convergence properties of the fixed-stress iterative solution of two-way coupled poromechanics. International Journal for Numerical and Analytical Methods in Geomechanics, 39(14):1593–1618, 2015b. 133 Zhi Chai, Hewei Tang, Youwei He, John Killough, and Yuhe Wang. Uncertainty quantification of the fracture network with a novel fractured reservoir forward model. In SPE annual technical conference and exhibition. Society of Petroleum Engineers, 2018. C. Y. Chen and E. Meiburg. Miscible porous media displacements in the quarter five-spot configu- ration. Part 2. Effect of heterogeneities. J. Fluid Mech., 371:269–299, 1998. Chukwudi Chukwudozie, Blaise Bourdin, and Keita Yoshioka. A variational phase-field model for hydraulic fracturing in porous media. Computer Methods in Applied Mechanics and Engineering, 347:957–982, 2019. Olaf A Cirpka and Albert J Valocchi. Two-dimensional concentration distribution for mixing- controlled bioreactive transport in steady state. Advances in water resources, 30(6-7):1668–1679, 2007. F. J. Cornaton, Y. J. Park, S. D. Normani, E. A. Sudicky, and J. F. Sykes. Use of groundwater lifetimeexpectancyfortheperformanceassessmentofadeepgeologicwasterepository: 1.Theory, illustrations, and implications. Water Resour. Res., 44, 2008. {10.1029/2007WR006208}. Olivier Coussy. Poromechanics. John Wiley & Sons, 2004. A. Daneshy. Dynamic active fracture interaction (DAFI) in horizontal wells. Hydraulic Fracturing J., 2(2):8–19, 2015. A. De Wit and G. M. Homsy. Viscous fingering in periodically heterogeneous porous media. II. Numerical simulations. J. Chem. Phys., 107:9619–9628, 1997. Marco Dentz, Tanguy Le Borgne, Andreas Englert, and Branko Bijeljic. Mixing, spreading and reaction in heterogeneous media: A brief review. Journal of contaminant hydrology, 120:1–17, 2011. Emmanuel Detournay. Propagation regimes of fluid-driven fractures in impermeable rocks. Inter- national Journal of Geomechanics, 4(1):35–45, 2004. V. Eswaran and S. B. Pope. Direct numerical simulations of the turbulent mixing of a passive scalar. Phys. Fluids, 31:506–520, 1988. 10.2118/21237-PA. M. Ferronato, N. Castelletto, and G. Gambolati. A fully coupled 3-D mixed finite element model of Biot consolidation. J. Comput. Phys., 229(12):4813–4830, 2010. Tracey C Flowers and James R Hunt. Viscous and gravitational contributions to mixing during vertical brine transport in water-saturated porous media. Water resources research, 43(1), 2007. P. J. Fox, J. Lee, and J. J. Lenhart. Coupled consolidation and contaminant transport in compress- ible porous media. Int. J. Geomech., 11:113–123, 2011a. 134 PatrickJFox, JangguenLee, andJohnJLenhart. Coupledconsolidationandcontaminanttransport in compressible porous media. International Journal of Geomechanics, 11(2):113–123, 2011b. Anthony F. Gangi. Variation of whole and fractured porous rock permeability with confining pres- sure. International Journal of Rock Mechanics and Mining Sciences Geomechanics Abstracts, 15: 249–257, 1978. https://doi.org/10.1016/0148-9062(78)90957-9. D. Gawin, P. Baggio, and B. A. Schrefler. Coupled heat, water and gas flow in deformable porous media. Int. J. Numer. Methods Fluids, 20:969–987, 1995. 10.1002/fld.1650200817. R. E. Goodman, R. L. Taylor, and T. L. Brekke. A model for the mechanics of jointed rock. J. Soil Mech., 94:637–659, 1968. Somdatta Goswami, Cosmin Anitescu, and Timon Rabczuk. Adaptive fourth-order phase field analysis for brittle fracture. Computer Methods in Applied Mechanics and Engineering, 361: 112808, 2020. Xuyang Guo, Kan Wu, Cheng An, Jizhou Tang, and John Killough. Numerical investigation of effects of subsequent parent-well injection on interwell fracturing interference using reservoir- geomechanics-fracturing modeling. SPE Journal, 24(04):1–884, 2019. Hadi Hajibeygi, Dimitris Karvounis, and Patrick Jenny. A hierarchical fracture model for the iterative multiscale finite volume method. Journal of Computational Physics, 230(24):8729–8743, 2011. Nico J Hardebol, Christine Maier, Hamid Nick, Sebastian Geiger, Giovanni Bertotti, and Herman Boro. Multiscale fracture network characterization and impact on flow: A case study on the latemar carbonate platform. Journal of Geophysical Research: Solid Earth, 120(12):8197–8222, 2015. Timo Heister, Mary F Wheeler, and Thomas Wick. A primal-dual active set method and predictor- corrector mesh adaptivity for computing fracture propagation using a phase-field approach. Com- puter Methods in Applied Mechanics and Engineering, 290:466–495, 2015. George M Homsy. Viscous fingering in porous media. Annual review of fluid mechanics , 19(1): 271–311, 1987. Qinhong Hu and Jean E. Moran. Simultaneous analyses and applications of multiple fluorobenzoate and halide tracers in hydrologic studies. Hydrological Processes, 19(14):2671–2687, 2005. https:// doi.org/10.1002/hyp.5780. RA Hull, Robert Meek, Hector Bello, and Douglas Miller. Case history of das fiber-based mi- croseismic and strain data, monitoring horizontal hydraulic stimulations using various tools to highlight physical deformation processes (part a). In Unconventional Resources Technology Con- ference, Austin, Texas, 24-26 July 2017, pages 3050–3062. Society of Exploration Geophysicists, 135 American Association of Petroleum ..., 2017. Peter S. Huyakorn, Peter F. Andersen, James W. Mercer, and Harold O. White Jr. Saltwater intrusion in aquifers: Development and testing of a three-dimensional finite element model. Water Res. Res., 23:293–312, 1987. 10.1029/WR023i002p00293. Martin Iding and Philip Ringrose. Evaluating the impact of fractures on the performance of the In Salah CO2 storage site. International Journal of Greenhouse Gas Control, 4(2):242–248, 2010. Gunnar Jansen, Benoît Valley, and Stephen A Miller. Thermaid-a matlab package for thermo- hydraulic modeling and fracture stability analysis in fractured reservoirs. arXiv preprint arXiv:1806.10942, 2018. B. Jha. Flow through porous media: from mixing of fluids to triggering of earthquakes . PhD thesis, Massachusetts Institute of Technology, 2013. B. Jha and R. Juanes. A locally conservative finite element framework for the simulation of coupled flow and reservoir geomechanics. Acta Geotech., 2:139–153, 2007. B. Jha, F. Bottazzi, R. Wojcik, M. Coccia, N. Bechor, D. Mclaughlin, T. A. Herring, B. H. Hager, S. Mantica, and R. Juanes. Reservoir characterization in an underground gas storage field using joint inversion of flow and geodetic data. Int. J. Numer. Anal. Methods Geomech., 39:1619–1638, 2015. 10.1002/2013WR015175. Birendra Jha and Ruben Juanes. Coupled multiphase flow and poromechanics: A computational model of pore pressure effects on fault slip and earthquake triggering. Water Resour. Res., 50(5): 3776–3808, 2014. Birendra Jha, Luis Cueto-Felgueroso, and Ruben Juanes. Fluid mixing from viscous fingering. Phys. Rev. Lett., 106(19):194502, 2011a. Birendra Jha, Luis Cueto-Felgueroso, and Ruben Juanes. Quantifying mixing in viscously unstable porous media flows. Physical Review E, 84(6):066312, 2011b. Birendra Jha, Luis Cueto-Felgueroso, and Ruben Juanes. Synergetic fluid mixing from viscous fingering and alternating injection. Phys. Rev. Lett., 111(14):144501, 2013. Joaquín Jiménez-Martínez, Mark L Porter, Jeffrey D Hyman, J William Carey, and Hari S Viswanathan. Mixing in a three-phase system: Enhanced production of oil-wet reservoirs by CO2 injection. Geophysical Research Letters, 43:196–205, 2016. Lei Jin and Mark D Zoback. Fully dynamic spontaneous rupture due to quasi-static pore pressure andporoelasticeffects: Animplicitnonlinearcomputationalmodeloffluid-inducedseismicevents. J. Geophys. Res. Solid Earth, 123:9430–9468, 2018. 136 W. P. Jones and B. E. Launder. The prediction of laminarization with a two-equation model of turbulence. Int. J. Heat Mass Trans., 15:301–304, 1972. G. Kampfer and M. Dawson. A novel approach to mapping hydraulic fractures using poromechanic principles. In 50th US Rock Mechanics/Geomechanics Symposium, pages ARMA 16–843. Am. rock Mech. Assoc., 2016a. Günther Kampfer and Matthew Dawson. A novel approach to mapping hydraulic fractures us- ing poromechanic principles. In ARMA US Rock Mechanics/Geomechanics Symposium, pages ARMA–2016. ARMA, 2016b. V. Kapoor and L. W. Gelhar. Transport in three-dimensionally heterogeneous aquifers 1. Dynamics of concentration fluctuations. Water Resour. Res., 30:1775, 1994. V. Kapoor and P. K. Kitanidis. Concentration fluctuations and dilution in aquifers. Water Resour. Res., 34:1181–1193, 1998. Mohammad Karimi-Fard, Louis J. Durlofsky, and Khalid Aziz. An efficient discrete-fracture model applicable for general-purpose reservoir simulators. Spe Journal, 9:227–236, 2004. URL https:// api.semanticscholar.org/CorpusID:131137515. E. J. Koval. A method for predicting the performance of unstable miscible displacement in hetero- geneous media. Soc. Pet. Eng. J., 3(02):145–154, 1963. L. W. Lake, R. Johns, B. Rossen, and G. Pope. Fundamentals of Enhanced Oil Recovery. Society of Petroleum Engineers, 2014. A Lamur, JE Kendrick, GH Eggertsson, RJ Wall, JD Ashworth, and Y Lavallée. The permeability of fractured rocks in pressurised volcanic and geothermal systems. Scientific reports , 7(1):1–9, 2017. John-Paul Latham, Jiansheng Xiang, Mandefro Belayneh, Hamidreza M. Nick, Chin-Fu Tsang, and Martin J. Blunt. Modelling stress-dependent permeability in fractured rock including effects of propagatingandbendingfractures. International Journal of Rock Mechanics and Mining Sciences, 57:100–112, 2013. 10.1016/j.ijrmms.2012.08.002. Hamed Lawal, Greg Jackson, Nnamdi Abolo, and Cecilia Flores. A novel approach to modeling and forecasting frac hits in shale gas wells. In SPE Europec featured at EAGE Conference and Exhibition?, pages SPE–164898. SPE, 2013. T. Le Borgne and P. Gouze. Non-Fickian dispersion in porous media: 2. Model validation from measurements at different scales. Water Resour. Res., 44:W06427, 2008. T. Le Borgne, M. Dentz, and E. Villermaux. Stretching, coalescence, and mixing in porous media. Phys. Rev. Lett., 110:1468–1475, 2013. 137 Tanguy Le Borgne, Marco Dentz, Diogo Bolster, Jesus Carrera, Jean-Raynald De Dreuzy, and Philippe Davy. Non-fickian mixing: Temporal evolution of the scalar dissipation rate in hetero- geneous porous media. Adv. Water Resour., 33(12):1468–1475, 2010. Seong H Lee, MF Lough, and CL Jensen. Hierarchical modeling of flow in naturally fractured formations with multiple length scales. Water resources research, 37(3):443–455, 2001a. Seong Hyuk Lee, Michael F. Lough, and Clair Jensen. Hierarchical modeling of flow in naturally fractured formations with multiple length scales. Water Resources Research, 37:443–455, 2001b. Sanjiva K Lele. Compact finite difference schemes with spectral-like resolution. J. Comput. Phys., 103(1):16–42, 1992. R. J. LeVeque. Finite Volume Methods for Hyperbolic Problems. Cambridge University Press, UK, 2002. Jana K. Levison and Kent S. Novakowski. Rapid transport from the surface to wells in fractured rock: A unique infiltration tracer experiment. Journal of Contaminant Hydrology, 131(1):29–38, 2012. ISSN 0169-7722. https://doi.org/10.1016/j.jconhyd.2012.01.001. R.W.LewisandB.A.Schrefler. The Finite Element Method in the Static and Dynamic Deformation and Consolidation of Porous Media. Wiley, 2nd edition, 1998. R. W. Lewis and Y. Sukirman. Finite element modelling of three-phase flow in deforming saturated oil reservoirs. Int. J. Numer. Anal. Methods Geomech., 17:577–598, 1993a. R. W. Lewis and Y. Sukirman. Finite element modeling for simulating the surface subsidence above a compacting hydrocarbon reservoir. Int. J. Numer. Anal. Methods Geomech., 18:619–639, 1993b. R. W. Lewis, A. Makurat, and W. K. S. Pao. Fully coupled modeling of seabed subsidence and reservoir compaction of North Sea oil fields. Hydrogeol. J., 11(1):142–161, 2003. Junjian Li, Yanli Pei, Hanqiao Jiang, Lin Zhao, Linkai Li, He Zhou, Yuyun Zhao, and Zhentao Zhang. Tracer flowback based fracture network characterization in hydraulic fracturing. In Abu Dhabi International Petroleum Exhibition & Conference. Society of Petroleum Engineers, 2016. S. Li, D. Zhang, and X. Li. A new approach to the modeling of hydraulic-fracturing treatments in naturally fractured reservoirs. SPE Journal, 22(04):1–64, 2017. X. Li, X. Li, D. Zhang, and R. Yu. A dual-grid, implicit, and sequentially coupled geomechanics- and-composition model for fractured reservoir simulation. SPE Journal, 2020. Yongzan Liu, Kan Wu, Ge Jin, and George Moridis. Rock Deformation and Strain-Rate Charac- terization during Hydraulic Fracturing Treatments: Insights for Interpretation of Low-Frequency Distributed Acoustic-Sensing Signals. SPE Journal, 25(05):2251–2264, 2020. 10.2118/202482-PA. 138 D. Loggia, N. Rakotomalala, D. Salin, and Y. C. Yortsos. Evidence of new instability thresholds in miscible displacements in porous media. Europhys. Lett., 32:633–638, 1995. G.L. Manjunath, Zhongqi Liu, and Birendra Jha. Multi-stage hydraulic fracture monitoring at the lab scale. Engineering Fracture Mechanics, 289:109448, 2023. https://doi.org/10.1016/ j.engfracmech.2023.109448. Stephan K Matthäi, Hamidreza M Nick, Christopher Pain, and Insa Neuweiler. Simulation of solute transport through fractured rock: a higher-order accurate finite-element finite-volume method permitting large time steps. Transport in porous media, 83(2):289–318, 2010. John McCarthy and John Zachara. Subsurface transport of contaminants. Environmental science & technology, 23(5):496–502, 1989. S. Meguerdijian and B. Jha. Quantification of fault leakage dynamics based on leakage magnitude and dip angle. Int. J. Numer. Anal. Methods Geomech., 2021. 10.1002/nag.3267. James W Mercer and Robert M Cohen. A review of immiscible fluids in the subsurface: Properties, models, characterization and remediation. Journal of contaminant hydrology, 6(2):107–163, 1990. Bert Metz, Ogunlade Davidson, HC De Coninck, Manuela Loos, and Leo Meyer. IPCC special report on carbon dioxide capture and storage. Cambridge: Cambridge University Press, 2005. Christian Miehe and Steffen Mauthe. Phase field modeling of fracture in multi-physics problems. part iii. crack driving forces in hydro-poro-elasticity and hydraulic fracturing of fluid-saturated porous media. Computer Methods in Applied Mechanics and Engineering, 304:619–655, 2016. Christian Miehe, Martina Hofacker, and Fabian Welschinger. A phase field model for rate- independent crack propagation: Robust algorithmic implementation based on operator splits. Computer Methods in Applied Mechanics and Engineering, 199(45-48):2765–2778, 2010. A Mikelić, Mary F Wheeler, and Thomas Wick. Phase-field modeling of a fluid-driven fracture in a poroelastic medium. Computational Geosciences, 19(6):1171–1195, 2015. Nicolas Moës, John E. Dolbow, and Ted Belytschko. A finite element method for crack growth without remeshing. International Journal for Numerical Methods in Engineering, 46:131–150, 1999. Ali Moinfar, Kamy Sepehrnoori, Russell T Johns, and Abdoljalil Varavei. Coupled geomechanics and flow simulation for an embedded discrete fracture model. In SPE Reservoir Simulation Symposium. Society of Petroleum Engineers, 2013. Lawrence C Murdoch and Leonid N Germanovich. Analysis of a deformable fracture in permeable material. International Journal for Numerical and Analytical Methods in Geomechanics, 30(6): 529–561, 2006. 139 Evgeniy Myshakin, Hema Siriwardane, Carter Hulcher, Ernest Lindner, Neal Sams, Seth King, and Mark McKoy. Numerical simulations of vertical growth of hydraulic fractures and brine migration in geological formations above the Marcellus shale. Journal of Natural Gas Science and Engineering, 27:531–544, 2015. SundararajanNatarajan, RatnaKAnnabattula, andEmilioMartínez-Pañeda. Phasefieldmodelling of crack propagation in functionally graded materials. Composites Part B: Engineering, 169:239– 248, 2019. Lin Ni, Xue Zhang, Liangchao Zou, and Jinsong Huang. Phase-field modeling of hydraulic fracture network propagation in poroelastic rocks. Computational Geosciences, 24(5):1767–1782, 2020. Christos Nicolaides, Birendra Jha, Luis Cueto-Felgueroso, and Ruben Juanes. Impact of viscous fingering and permeability heterogeneity on fluid mixing in porous media. Water Resources Research, 51(4):2634–2647, 2015. Jack H Norbeck, Mark W McClure, Jonathan W Lo, and Roland N Horne. An embedded fracture modeling framework for simulation of hydraulic fracturing and shear stimulation. Computational Geosciences, 20(1):1–18, 2016. C. Ollivier-Gooch and M. Van Altena. A high-order-accurate unstructured mesh finite-volume scheme for the advection–diffusion equation. J. Comput. Phys., 181:729–752, 2002. R Olsson and N Barton. An improved model for hydromechanical coupling during shearing of rock joints. International Journal of Rock Mechanics and Mining Sciences, 38(3):317–329, 2001. W. K. S. Pao, R. W. Lewis, and I. Masters. A fully coupled hydro-thermo-poro-mechanical model for black oil reservoir simulation. Int. J. Numer. Anal. Methods Geomech., 25(12):1229–1256, 2001. Donald W Peaceman. Interpretation of well-block pressures in numerical reservoir simulation (in- cludes associated paper 6988). Society of Petroleum Engineers Journal, 18(03):183–194, 1978. G. P. Peters and D. W. Smith. Solute transport through a deforming porous medium. Int. J. Numer. Anal. Methods Geomech., 26:683–717, 2002. Justin Pogacnik, Derek Elsworth, Michael O’Sullivan, and John O’Sullivan. A damage mechanics approach to the simulation of hydraulic fracturing/shearing around a geothermal injection well. Computers and Geotechnics, 71:338–351, 2016. S. B. Pope. Turbulent Flows. Cambridge University Press, Cambridge, MA, 2000. M. M. Poulsen and B. H. Kueper. A field experiment to study the behavior of tetrachloroethylene in unsaturated porous media. Environ. Sci. Tech., 26:889–895, 1992. 140 GuotongRenandRamiMYounis. Anintegratednumericalmodelforcoupledporo-hydro-mechanics and fracture propagation using embedded meshes. Computer Methods in Applied Mechanics and Engineering, 376:113606, 2021. P Richter, T Parker, C Woerpel, Y Wu, R Rufino, and M Farhadiroushan. Hydraulic fracture monitoring and optimization in unconventional completions using a high-resolution engineered fibre-optic Distributed Acoustic Sensor. First Break, 37(4):63–68, 2019. R. Roded, X. Paredes, and R. Holtzman. Reactive transport under stress: Permeability evolution in deformable porous media. Earth Planet. Sci. Lett., 493:198–207, 2018. Guan Rong, Jie Yang, Long Cheng, and Chuangbing Zhou. Laboratory investigation of nonlinear flow characteristics in rough fractures during shear process. Journal of Hydrology, 541:1385–1394, 2016. https://doi.org/10.1016/j.jhydrol.2016.08.043. N. P. Roussel and S. Agrawal. Introduction to poroelastic response monitoring – quantifying hy- draulic fracture geometry and SRV permeability from offset-well pressure data. In 50th US Rock Mechanics/Geomechanics Symposium, page 2645414. Unconventional Resources Technology Con- ference, 2017a. 10.15530/urtec-2017-2645414. Nicolas P Roussel and Samarth Agrawal. Introduction to poroelastic response monitoring- quantifying hydraulic fracture geometry and srv permeability from offset-well pressure data. In SPE/AAPG/SEG Unconventional Resources Technology Conference. OnePetro, 2017b. D.F. Rugh and T.J. Burbey. Using saline tracers to evaluate preferential recharge in fractured rocks, floyd county, virginia, usa. Hydrogeol J., 16:251–262, 2008. 10.1007/s10040-007-0236-3. Jonny Rutqvist and Ove Stephansson. The role of hydromechanical coupling in fractured rock engineering. Hydrogeology Journal, 11(1):7–40, 2003. Jonny Rutqvist, Colin Leung, Andrew Hoch, Yuan Wang, and Zhen Wang. Linked multicontinuum and crack tensor approach for modeling of coupled geomechanics, fluid flow and transport in fractured rock. J. Rock Mech. Geotech. Eng., 5(1):18–31, 2013. M. Sahimi. Flow and Transport in Porous Media and Fractured Rock: From Classical Methods to Modern Approaches. Wiley–VCH, Germany, 2nd edition, 2011. Anusarn Sangnimnuan, Jiawei Li, and Kan Wu. Development of efficiently coupled fluid- flow/geomechanics model to predict stress evolution in unconventional reservoirs with complex- fracture geometry. SPE Journal, 23(03):640–660, 2018. David Santillán, Ruben Juanes, and Luis Cueto-Felgueroso. Phase field model of hydraulic fractur- ing in poroelastic media: Fracture propagation, arrest, and branching under fluid injection and extraction. Journal of Geophysical Research: Solid Earth, 123(3):2127–2155, 2018. 141 Nicolas Schwenck, Bernd Flemisch, Rainer Helmig, and Barbara Wohlmuth. Dimensionally reduced flow models in fractured porous media: crossings and boundaries. Computational Geosciences, 19, 12 2015. 10.1007/s10596-015-9536-1. A. Settari and D. A. Walters. Advances in coupled geomechanical and reservoir modeling with applications to reservoir compaction. Soc. Pet. Eng. J., 6:334–342, 2001. Mahmood Shakiba. Modeling and simulation of fluid flow in naturally and hydraulically fractured reservoirs using embedded discrete fracture model (EDFM). PhD thesis, University of Texas Austin, 2014. Rilin Shen, Haim Waisman, and Licheng Guo. Fracture of viscoelastic solids modeled with a modified phase field method. Computer Methods in Applied Mechanics and Engineering, 346: 862–890, 2019. Erik Svenson, Todd Schweisinger, and Lawrence C. Murdoch. Field evaluation of the hydrome- chanical behavior of flat-lying fractures during slug tests. Journal of Hydrology, 359:30–45, 2008. 10.1016/j.jhydrol.2008.06.004. M. R. Sweeney and J. D. Hyman. Stress effects on flow and transport in three-dimensional fracture networks. Journal of Geophysical Research: Solid Earth, 125(8):e2020JB019754, 2020. https:// doi.org/10.1029/2020JB019754. C.-T. Tan and G. M. Homsy. Stability of miscible displacements in porous media: Rectilinear flow. Phys. Fluids, 29:3549, 1986. C. T. Tan and G. M. Homsy. Simulation of nonlinear viscous fingering in miscible displacement. Phys. Fluids, 31:1330–1338, 1988. C.-T. Tan and G. M. Homsy. Viscous fingering with permeability heterogeneity. Phys. Fluids, 4(6): 1089–1093, 1992. Matei Tene, Sebastian BM Bosma, Mohammed Saad Al Kobaisi, and Hadi Hajibeygi. Projection- based embedded discrete fracture model (pedfm). Advances in Water Resources, 105:205–216, 2017. M. Tran, F. Aminzadeh, and B. Jha. Effect of coupled flow and geomechanics on transport of a fluid slug in a stress-sensitive reservoir. In 52nd US Rock Mechanics/Geomechanics Symposium, pages ARMA–18/143. American Rock Mechanics Association, 2018. Minh Tran and Birendra Jha. Coupling between transport and geomechanics affects spreading and mixing during viscous fingering in deformable aquifers. Advances in Water Resources, 136:103485, 2020. Minh Tran and Birendra Jha. Effect of poroelastic coupling and fracture dynamics on solute 142 transport and geomechanical stability. Water Resources Research, 57(10):e2021WR029584, 2021. https://doi.org/10.1029/2021WR029584. J.R.Waggoner, J.L.Castillo, andL.L.Lake. Simulationofeorprocessesinstochasticallygenerated permeable media. Soc. Petrol. Engrs Form. Eval., 7:9619–9628, 1992. 10.2118/21237-PA. H. F. Wang. Theory of Linear Poroelasticity with Applications to Geomechanics and Hydrogeology. Princeton University Press, 2000. Jing Wang, Haishan Luo, Huiqing Liu, Fei Cao, Zhitao Li, and Kamy Sepehrnoori. An integrative model to simulate gas transport and production coupled with gas adsorption, non-Darcy flow, surface diffusion, and stress dependence in organic-shale reservoirs. SPE Journal, 22(01):244–264, 2017. NR Warner, TH Darrah, RB Jackson, Romain Millot, Wolfram Kloppmann, and Avner Vengosh. Newtracersidentifyhydraulicfracturingfluidsandaccidentalreleasesfromoilandgasoperations. Environmental science & technology, 48(21):12552–12560, 2014. Noriaki Watanabe, Nobuo Hirano, and Noriyoshi Tsuchiya. Determination of aperture structure and fluid flow in a rock fracture by high-resolution numerical modeling on the basis of a flow-through experiment under confining pressure. Water Resources Research, 44(6), 2008. C. Welty and L. W. Gelhar. Stochastic analysis of the effects of fluid density and viscosity variability on macrodispersion in heterogeneous porous media. Water Resour. Res., 27:2061, 1991. J. A. White and R. I. Borja. Stabilized low-order finite elements for coupled solid-deformation/fluid- diffusion and their application to fault zone transients. Comput. Geosci., 197:4353–4366, 2007. Thomas Wick, Gurpreet Singh, and Mary F Wheeler. Fluid-filled fracture propagation with a phase-field approach and coupling to a reservoir simulator. SPE Journal, 21(03):981–999, 2016. Paul A Witherspoon, NG Cook, and JE Gale. Geologic storage of radioactive waste: Field studies in sweden. Science, 211(4485):894–900, 1981. Jian-Ying Wu, Vinh Phu Nguyen, Chi Thanh Nguyen, Danas Sutula, Sina Sinaie, and Stéphane PA Bordas. Phase-field modeling of fracture. Advances in Applied Mechanics, 53:1–183, 2020. LiangXia, JulienYvonnet, andSiavashGhabezloo. Phasefieldmodelingofhydraulicfracturingwith interfacial damage in highly heterogeneous fluid-saturated porous media. Engineering Fracture Mechanics, 186:158–180, 2017. Yifei Xu. Implementation and application of the embedded discrete fracture model (EDFM) for reservoir simulation in fractured reservoirs. PhD thesis, University of Texas Austin, 2015. Xia Yan, Zhaoqin Huang, Jun Yao, Zhao Zhang, Piyang Liu, Yang Li, and Dongyan Fan. Numerical 143 simulation of hydro-mechanical coupling in fractured vuggy porous media using the equivalent continuum model and embedded discrete fracture model. Advances in Water Resources, 126: 137–154, 2019. 10.1016/j.advwatres.2019.02.013. Zhi Ye and Ahmad Ghassemi. Injection-induced shear slip and permeability enhancement in granite fractures. Journal of Geophysical Research: Solid Earth, 123(10):9009–9032, 2018. Liang-Ping Yi, Haim Waisman, Zhao-Zhong Yang, and Xiao-Gang Li. A consistent phase field model for hydraulic fracture propagation in poroelastic media. Computer Methods in Applied Mechanics and Engineering, 372:113396, 2020. Y. C. Yortsos and F. J. Hickernell. Linear stability of immiscible displacement in porous media. SIAM J Appl. Math, 49:730–749, 1989. H. J. Zhang, D.-S. Jeng, B.R. Seymour, D. A. Barry, and L. Li. Solute transport in partially- saturated deformable porous media: Application to a landfill clay liner. Adv. Water Resour., 40: 1–10, 2012. X. Zhao and B. Jha. Role of well operations and multiphase geomechanics in controlling fault stability during CO 2 storage and enhanced oil recovery. J. Geophys. Res. Solid Earth, 124, 2019. 10.1029/2019JB017298. X. Zhao and B. Jha. A new coupled multiphase flow–finite strain deformation–fault slip framework for induced seismicity. J. Comput. Phys., 433, 2021. 10.1016/j.jcp.2021.110178. Zhihong Zhao, Lanru Jing, Ivars Neretnieks, and Luis Moreno. Numerical modeling of stress effects on solute transport in fractured rocks. Computers and Geotechnics, 38(2):113–126, 2011. W. B. Zimmerman and G. M. Homsy. Three-dimensional viscous fingering: A numerical study. Phys. Fluids, 4:1901–1914, 1992. 144

Abstract (if available)

Abstract The main focus of my research is to develop a holistic computational modeling framework that couples reservoir geomechanics, fluid flow, solvent transport, and fracture dynamics. The first objective is to quantitatively characterize the effect of geomechanical coupling on key macroscopic transport phenomena such as spreading, mixing, and viscous fingering. The second objective is to investigate how flow-transport coupling and fracture dynamics modulate the stress state and geomechanical stability of stress-sensitive subsurface formations. The third objective is to probe the impact of poroelasticity on the temporal and spatial evolution of embedded fracture networks in terms of induced permeability, dynamic propagation, and fracture interactions.

The key conclusions are summarized here. I examined how viscous fingering-induced fluctuations in the flow field can increase the strength of coupling between poroelastic stresses and the mixing/spreading characteristics of the tracer slug. Failure to account for the stress-transport coupling will result in an inaccurate prediction of the tracer breakthrough time, tracer recovery, and the degree of mixing. I characterized the effect of geomechanical coupling on key transport metrics and quantified the impact each coupling parameter had on tracer breakthrough and the degree of mixing. I explored how the characteristics of the tracer retrieval curve were modulated by both fracture morphology and the strength of geomechanical coupling. I discovered that the transport-geomechanics synergy exists and can be modeled via numerical coupling to the pressure field. I also found that the geomechanical coupling strength governed the evolution of fracture permeability with different quantifiable impacts. I investigated how flow-transport coupling and fracture dynamics influence the stress state and fracture stability in the domain. Last but not least, I had shown that the poroelastic pressure signal resulting from deformation could serve as a diagnostic tool for identifying fracture morphology and reservoir heterogeneity.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Reactivation of multiple faults in oilfields with injection and production

$A coupled poromechanics-damage mechanics approach to model fracturing in multiphase flow$

PDF

A coupled poromechanics-damage mechanics approach to model fracturing in multiphase flow

$Optimization of CO2 storage efficiency under geomechanical risks using coupled flow-geomechanics-fracturing model$

PDF

Optimization of CO2 storage efficiency under geomechanical risks using coupled flow-geomechanics-fracturing model

$An extended finite element method based modeling of hydraulic fracturing$

PDF

An extended finite element method based modeling of hydraulic fracturing

PDF

A coupled multiphase flow-geomechanics framework for modeling material, kinematic and contact nonlinearities in faulted reservoirs

$Coarse-scale simulation of heterogeneous reservoirs and multi-fractured horizontal wells$

PDF

Coarse-scale simulation of heterogeneous reservoirs and multi-fractured horizontal wells

PDF

On the transport dynamics of miscible fluids in porous media under different sources of disorder

PDF

Stress and deformation analysis on fluid-exposed reservoir rocks

PDF

Effective flow and transport properties of deforming porous media and materials: theoretical modeling and comparison with experimental data

PDF

Continuum modeling of reservoir permeability enhancement and rock degradation during pressurized injection

PDF

Modeling and simulation of complex recovery processes

PDF

The role of counter-current flow in the modeling and simulation of multi-phase flow in porous media

PDF

Asphaltene deposition during co₂ injection in enhanced oil recovery applications

PDF

Efficient simulation of flow and transport in complex images of porous materials and media using curvelet transformation

PDF

Flow and thermal transport at porous interfaces

PDF

Modeling and simulation of multicomponent mass transfer in tight dual-porosity systems (unconventional)

PDF

Mechanical and flow interactions facilitate cooperative transport and collective locomotion in animal groups

PDF

Chemical and mechanical deformation of porous media and materials during adsorption and fluid flow

PDF

Investigation of gas transport and sorption in shales

PDF

A study of diffusive mass transfer in tight dual-porosity systems (unconventional)

Asset Metadata

Creator Tran, Minh Tuan (author)

Core Title Synergistic coupling between geomechanics, flow, and transport in fractured porous media: applications in hydraulic fracturing and fluid mixing

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Petroleum Engineering

Degree Conferral Date 2023-12

Publication Date 09/11/2023

Defense Date 10/11/2021

Publisher Los Angeles, California (original), University of Southern California (original), University of Southern California. Libraries (digital)

Tag flow transport in porous media,fluid mixing,fracture dynamics,geomechanics,hydraulic fracturing,numerical simulation,OAI-PMH Harvest,poroelasticity

Format theses (aat)

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Jha, Birendra (committee chair), Ben-zion, Yehuda (committee member), Elsworth, Derek (committee member), Ershaghi, Iraj (committee member), Jessen, Kristian (committee member)

Creator Email trademark152@gmail.com,tranmt@usc.edu

Permanent Link (DOI) https://doi.org/10.25549/usctheses-oUC113305014

Unique identifier UC113305014

Identifier etd-TranMinhTu-12335.pdf (filename)

Legacy Identifier etd-TranMinhTu-12335

Document Type Dissertation

Format theses (aat)

Rights Tran, Minh Tuan

Internet Media Type application/pdf

Type texts

Source 20230911-usctheses-batch-1093 (batch), University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the author, as the original true and official version of the work, but does not grant the reader permission to use the work if the desired use is covered by copyright. It is the author, as rights holder, who must provide use permission if such use is covered by copyright.

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA

Repository Email cisadmin@lib.usc.edu