The effect of ﬂuid pulse driven fractures (FPDF) propagating in poroelastic media on fault slip in the presence of natural fractures is a complicated in- terplay between fracture propagation, fracture-fracture interaction, fracture- fault interaction, friction model governing fault slip and wave propagation associated with pulsing injection. Furthermore, the problem is stochastic due to the uncertainty associated with the existing fracture-fault topology. 1 Introduction Figure 1: States where hydrofracturing in ultra-low permeability rocks is currently ongoing. Source: [1] Hydraulic fracturing (HF) in porous media is a complex multi-physics problem [2 –6] which involves nucleation and propagation of cracks by injecting a combination of padding ﬂuid and additives in the form of proppant at high pressure down a wellbore. HF is omnipresent in the United States (Fig. 1) due to its historical use for stimulation of oil and gas reservoirs [3, 7]. In the recent past, techniques like pulse or fatigue fracturing [8 –12] have been deployed in which a controlled dynamic release of ﬂuid energy pressurizes a wellbore to initiate multiple fractures, which are then propagated by static or dynamic pressurization. Pulsing been proposed as a strategy to enhance the productivity and mitigate the environmental impact of HF hydrocarbon projects due to lower quantities of fracturing ﬂuid required [11, 13]. Pulse fracturing in hot dry rock has been documented to have long lasting implications for enhanced geothermal energy (EGS) production [4, 14 –17]. Figure 2: Four conceptual models for EGS stimulation. The black dot represents the wellbore. New fractures are represented with red lines, and preexisting fractures are represented with blue lines. Source: [14] The conceptual models that have been proposed in literature as stimulation mechanisms in EGS [14] include (1) Only new, propagating fractures contribute to permeability en- hancement, (2) Stimulation occurs only through induced slip on preexisting fractures, (3) Continuous new fractures propagate away from the wellbore, but ﬂuid leaks off into natural fractures, which slip and experience enhanced transmissivity, and (4) Continuous pathways for ﬂow involve both new and preexisting fractures. The fourth mechanism has been widely accepted to model hydraulic stimulation of EGS [18–21]. The fracture initiated from a sheared preexisting natural fracture is called a wing crack or a splay fracture. Wing cracks found in rock are tensile fractures but they are different from hydraulic fractures because (1) wing cracks are initiated from the tension ﬁeld induced by shear slip of a preexisting natural fracture while hydraulic fractures are initiated from an injection well by ﬂuid pressure; (2) wing cracks are curved cracks while hydraulic fractures propagate straight and perpendicularly to the least principal stress; and (3) tension forces to open wing cracks are supported both by ﬂuid pressure and shear slip of a preexisting natural fracture while tension forces to open hydraulic fractures are supported only by ﬂuid pressure [21]. Hitherto, HF has been associated with microseismic events due to their small magnitudes. However, there is now growing evidence to suggest that HF induced fault activity in EGS projects across the world has caused larger earthquakes [18, 20, 22]. Fig. 3 shows the distribution of seismicity in the Pohang, South Korea EGS site. Postulates suggest the actual amount of injected ﬂuid was only 0:2%of what is theoretically expected for the largest 2 Figure 3: Seismicity in the Pohang EGS site. Source: [22] recorded earthquake at the site [22, 23]. The causation is established via spatio-temporal correlation; the events occur close ( <10km) to the injection well and shortly (minutes to hours) after injection. Most of the sites are in geographical regions with little to no history of prior natural seismicity, which further supports the causative link between fracturing and the observed seismic events. The HF-induced earthquake triggering mechanisms (see Fig. 4) that have been proposed are: (1) propagation of the ﬂuid pressure diffusion front through the newly created fracture network, followed by pressurization/lubrication of a basement fault that slips seismically [24], (2) poroelastic stress transfer from slipped fractures on to the fault causing an increase in shear and/or tension on the fault followed by seismic slip, and (3) ﬂuid injection-induced aseismic rupture front which propagates faster and to larger distances than the ﬂuid pressure diffusion front, thereby transmitting stresses to unstable portions of the fault beyond the ﬂuid-pressurized region [25]. All three mechanisms can activate a critically stressed fault (initial stress state close to frictional failure), which can culminate into an induced seismic event [1, 23, 26–32]. In quasi-static simulations, which neglect elastodynamics and wave propagation effects, the induced seismicity mechanisms are often quantiﬁed using Coulomb Failure Function (CFF) [33 –35], which includes contributions from the effective stress and fault’s friction coefﬁcient. The effective stress is deﬁned as the difference between the total stress (from tectonics and expansion/contraction of surrounding rock) and the pore pressure (from injection) scaled by the Biot coefﬁcient. The friction coefﬁcient evolution often requires a fault rheological model e.g., the Slip Weakening or the Rate- and State-dependent Friction (RSF) model. RSF is the gold standard for modeling earthquake cycles (inter-seismic loading followed by co-seismic relaxation) on faults [36 –38], whereas slip weakening provides a reasonable approximation for modeling the drop in friction during a single seismic/aseismic event. CFF is computed using the fault traction vector which is given by the product of the 3 Figure 4: Three proposed HF earthquake triggering mechanisms in literature. Source: [23] effective stress tensor and the fault normal vector. Initially, when the fault is locked and under equilibrium with the tectonic and hydrostatic stresses, CFF<0. As injection proceeds, the pore pressure increases and the near-wellbore region expands volumetrically, which applies compression on the region farther from the wellbore. The spatial distribution of the total stress tensor and pressure around the fault change. An increase in CFF, i.e. CFF>0, indicates ﬂuid-induced destabilization of the fault. When CFF exceeds fault’s intrinsic cohesion, the fault may slip. The slip releases the host rock strain energy, accumulated from natural and/or anthropogenic loading, over multiple years. The slip can be seismic or aseismic depending on fault’s frictional properties (e.g., velocity-weakening vs. velocity- strengthening) and the fault loading rate. Similarly, a decrease in fault’s CFF, e.g., due to stresses transferred from a newly created HF, can induce mechanical stabilization of the fault. 2 Literature Review The modeling and simulation of ﬂuid driven fracture is important for a variety of appli- cations in geotechnical engineering [39, 40], biomedical engineering [41, 42], and glaciol- ogy [43]. 2.1 Pulse/fatigue loading: In a typical fracturing job, a padding ﬂuid (mostly water) is ﬁrst pumped into the formation to generate fractures. Few hours later, a slurry with proppant (granular material like beach sand or ceramic particles) is injected into the fracture. The pad volume is usually around 10 to 50 percent of the total injected ﬂuid volume. In a non-pulse fracturing job, the prop- pant concentration monotonically increases to ensure the near-wellbore region maintains high fracture conductivity as the treatment continues. After 2012, the concentration of the proppant has increased in the industry from 0.5 lb/gal to 5 lb/gal [44, 45] to create fractures with larger aperture and permeability. In the pulsed proppant fracturing tech- nique, proppant-laden and proppant-free ﬂuid is intermittently pumped into the fracture to promote heterogeneous proppant placement inside the fracture [11, 12]. Such intermittent pumping is expected to create wave-like features in both fracture propagation and proppant transport with a hope of enhanced productivity. 4 2.2 Elastodynamics and ﬂuid inertia: The HF-induced seismicity problem is expected to be marked with multiple loading and unloading episodes such as dynamic HF propagation, ﬂuid diffusion, and fault-fracture interaction. High-rate transient phenomena such as fracture/fatigue pulsing is also ex- pected. These processes will likely generate dynamic stress wave fronts [46, 47] triggering stick-slip behavior of the fault [48]. Modeling these processes require a poro-elastodynamics framework with solid and ﬂuid inertia and the RSF fault rheology. Quasi-static simulations are limited to aseismic slip velocities ( <1m/s) and cannot be applied to faults that follow the RSF rheology, e.g. due to non-monotonic friction vs. slip behavior or due to multiple loading-unloading episodes [49, 50]. The quasi-dynamics approach [51], which approxi- mates solid inertia via a seismic radiation damping term containing a tunable viscosity parameter, has been used as a middle ground between quasi-static and full elastodynamics for modeling induced seismicity [52, 53]. However, the approach cannot model dynamic stress transfer and wave propagation effects [54], which are vital for modeling HF-fault interactions in the problem considered here, and suffer from poor accuracy in the simulated seismic magnitudes. Some of the prior studies are also limited by their assumption that the host rock is impermeable [52] or ﬂuid ﬂow is decoupled from rock mechanics [55, 56], both of which are not justiﬁed for the reservoir rocks we consider here. Although Darcy’s law is the de facto standard for the ﬂuid momentum equation for ﬂow in porous media in the limit of slow ﬂow [57 –60], there are multiple reasons why the slow ﬂow assumption does not work for HF under pulse or fatigue loading. Firstly, dimensional analyses [61, 62] and experimental studies [63, 64] reveal that hydraulic fractures propagate in essentially three regimes, out of which the ﬁrst and third regimes may violate the slow ﬂow assumption. The ﬁrst regime is the initiation regime, and the initial fracture velocity can be as high as 100m/s. Since the fracture velocity is dependent on the pumping rate in the third regime, the inertial effects on ﬂuid ﬂow cannot be ignored for pulse/fatigue loading. Secondly, recent studies [65 –68] suggest that including the solid acceleration causes oscillations in the ﬂuid pressure in the high permeability fracture zone. These oscillations are potentially numerical artifacts which can only be removed by incorporating ﬂuid inertia in the mathematical formulation. Although theoretical frameworks have been proposed [69, 70], current dynamic analyses of fractured porous media are often based on the inconsistent extension of Biot’s quasi- static theory that neglects ﬂuid inertia. For this reason, numerical methods and tools to solve dynamic hydraulic fracture problems do not yet exist, although some attempts have recently been made employing the extended ﬁnite element method [10, 71], cohesive zone models [72] and discontinuous pressure model [65]. 2.3 Coupling of proppant transport to dynamic HF: Proppants play a critical role in creating and maintaining the shape and size of hydraulic fractures by controlling the rheology of the ﬂuid inside HF. As a result, we hypothesize proppant transport exerts important control on HF-induced stresses. Yet most modeling studies, which account for HF propagation and ﬂow inside HF, ignore the dynamics of proppant ﬂow due to multiphysics modeling challenges. They approximate the effect of proppants on ﬂuid mobility via empirical functions [73]. Proppant modeling has been rec- ognized as one of the pressing issues in HF modeling [74]. Proppant is usually modeled by 5 representing the slurry as a two-component, interpenetrating continuum in which the ﬂow equations are solved for the mixture [44, 75 –77]. However, none of these studies investigate HF-induced seismicity and sometimes suffer from idealized assumed HF geometry [76]. The four equations that are solved in proppant transport coupled to HF are deformation, ﬂuid ﬂow, fracture propagation and proppant transport. The existing solution algorithms in literature do not iterate between the deformation/ﬂow/propagation and the proppant concentration, thereby yielding a loosely coupled scheme which lead to inaccuracies in the ﬁnal solution [78 –83]. The ﬁxed stress split iterative scheme [34, 84 –109] is leveraged in this realm. 2.4 Quantiﬁcation of induced seismicity: Elastodynamics [54, 110] and quasi-dynamic [52, 53] approaches allow calculation of the seismic magnitude from dynamic fault rupture simulations. In case of quasi-static simu- lations (inﬁnite wave speed assumption), with or without aseismic fault slip [34, 35, 111], it is common to use the CFF approach to quantify the likelihood of an induced seismic event. The CFF approach is based on the Mohr-Coulomb frictional failure theory, which deﬁnes CFF = j sj+f0 n, where j sjis the shear traction magnitude, 0 n=n+ pis the effective normal traction (positive under tension), nis the total normal traction, pis the pore pressure, is the Biot coefﬁcient, and fis the friction coefﬁcient. All quantities are deﬁned on the fault and 0 nand sare computed from the fault traction vector which is the product of the effective stress tensor and the fault normal vector. In quasi-dynamics, the radiation damping term, G_d=(2Vs), is added to the right hand side of the CFF equation, whereGis the shear modulus, Vsis the shear wave speed, and _dis the fault slip velocity along the rake direction. In case of post-mortem studies of induced seismicity [112 –115], the change in CFF, com- puted over a loading period, is often compared against the recorded seismic magnitude and shear stress released, to conﬁrm or deny the anthropogenic causality. It has been a popular approach to monitor and forecast induced seismicity because solving a quasi-static poroe- lastic problem in absence of dynamic slip is mathematically and numerically simpler than solving a poro-elastodynamics problem with RSF slip and possibly transient ﬂow effects. However, the CFF approach neglects dynamic stress transfer processes that we hypothesize play an important role in modulating the interaction between HF and faults in the induced seismicity problem considered in this proposal. Recent work on modeling injection-induced fault slip in conﬁrms that wave propagation is different between a quasi-static poroelastic system and a poro-elastodynamic system [110]. Recent studies [48,116] has proposed that stable fault creep in hydraulically fractured shales loads distant unstable parts of the fault capable of hosting an earthquake, even hundreds of metres above the hydraulic fracturing zone. They suggest that injection of fracturing ﬂuids can lead to slow slip on a fault, and over time, this creates sufﬁcient strain on another section of the fault to cause it to slip seismically. These studies postulate that following questions arise during elastodynamic extension of induced seismicity models: How far do these wave fronts traverse before they die out? Do these wave fronts get ampliﬁed or mitigated by natural fractures and faults? How do the reﬂected waves affect the inﬂuence of the incident waves? Do the seismic waves emanating from hypocenters of induced slip on a fault interact with hydraulic fractures? 6