Hydraulic Fracture Propagation in a Poro-Elastic Medium with Time-Dependent Injection Schedule Using the Time-Stepped Linear Superposition Method (TLSM)

: The Time-Stepped Linear Superposition Method (TLSM) has been used previously to model and analyze the propagation of multiple competitive hydraulic fractures with constant internal pressure loads. This paper extends the TLSM methodology, by including a time-dependent injection schedule using pressure data from a typical diagnostic fracture injection test (DFIT). In addition, the e ﬀ ect of poro-elasticity in reservoir rocks is accounted for in the TLSM models presented here. The propagation of multiple hydraulic fractures using TLSM-based codes preserves inﬁnite resolution by side-stepping grid reﬁnement. First, the TLSM methodology is brieﬂy outlined, together with the modiﬁcations required to account for variable time-dependent pressure and poro-elasticity in reservoir rock. Next, real world DFIT data are used in TLSM to model the propagation of multiple dynamic fractures and study the e ﬀ ect of time-dependent pressure and poro-elasticity on the development of hydraulic fracture networks. TLSM-based codes can quantify and visualize the e ﬀ ects of time-dependent pressure, and poro-elasticity can be e ﬀ ectively analyzed, using DFIT data, supported by dynamic visualizations of the changes in spatial stress concentrations during the fracture propagation process. The results from this study may help develop fracture treatment solutions with improved control of the fracture network created while avoiding the occurrence of fracture hits.


Introduction
Attempts to engineer the creation of fracture networks and their hydraulic conductivity are key for successful hydrocarbon field development strategies. To be successful, any hydraulic fracture intervention must consider the initial state of the reservoir (i.e., the local stress conditions, pore pressures, heterogeneities, anisotropy, stratification, faults, other natural fractures, and any pre-existing hydraulic fractures). A vast body of modeling attempts have focused on investigating how the various parameters affect the development of complex fracture networks during hydraulic fracturing [1,2]. Simultaneously, fracture characterization efforts in dedicated field experiments have revealed the existence of complex fracture networks in the vicinity of horizontal wells completed with multi-stage fracture treatment operations [3][4][5][6][7].

Pore Pressure and Localized Principal Stress Changes
What has remained perhaps poorly articulated in prior work is the careful distinction of the direct effect of pore pressure changes on (1) the local principal stress magnitudes, and (2) local principal stress orientations. First, the alleged changes in the orientation of principal stresses due to pore pressure alterations associated with engineering interventions are here grouped into Type A and Type B interventions. Type A intervention effects are short-term (days to weeks) and nearly instantaneous changes in the scalar reservoir pressure, accompanied by elastic strains in the rock adjacent to the main hydraulic fractures resulting in a local deviatoric stress tensor field, such as those induced by DFIT tests and full-blown fracture treatment operations. Type B intervention effects refer to the longer-term pressure changes (months to years) related to production-induced pressure depletion in the reservoir. When parent-child well situations are modeled, the effects of Type A and B interventions will start to interfere.
For example, the pressurization and depressurization of natural fractures due to injection and depletion in their vicinity can modify the magnitude of the in-situ stress. A hydraulic fracture propagation path due to Type A intervention would favor regions of higher local pore pressure, accompanied by locally increased tensile principal stress magnitudes that are superimposed on the in-situ stresses [18]. A nearby legacy or parent well, when being pressure-depleted due to a Type B Energies 2020, 13, 6474 3 of 22 intervention, may cause local changes in the stress magnitude as compared to the prior far-field stress magnitude. In the latter case, a reduction in the principal stress magnitude follows for a well under study. "Far-field stress" changes induced by an adjacent well may be a misnomer, and should be called "near-field stress" changes due to Type A and/or B interventions, as opposed to a regional native stress which traditionally is a familiar source of a "far-field stress".
Re-completions and infill drilling with Type A interventions may alter the stress state in the reservoir by the generation of locally superposed principal stresses, changing the prior far-field stresses. Empirical evidence for this contention comes from the observation that hydraulic fractures during a re-fracturing process propagated normal to that of the initial hydraulic fractures [19]. This change was attributed to a change in the horizontal stress orientations during production. Based on the new results presented in our study, we claim that a 90 degree flip in fracture propagation direction is simply due to a switch in the minor and major principal stress directions, which can be induced by either over-pressuring or under-pressuring a formation with regional far field stresses. This phenomenon has been extensively studied for wellbores with overbalanced or underbalanced pressure loads [20][21][22][23].
Prior work has unequivocally established that induced shear stresses in rocks may cause the initially horizontal principal stress to reorient out of the horizontal plane and thereby influences the direction of subsequent fractures, which has been observed to be escalated in heterogeneous media [24]. However, such shear stresses cannot be generated by pore pressure changes, which is why we suggest that prior claims that the local principal stress reorientation would occur due to poro-elasticity seem untenable, as follows further from evidence from experiments presented in our study.
While pore pressure can change due to Type A and Type B interventions, it may not act as the primary mechanism for stress reversal and fracture reorientation. In practice, local changes in the prior far-field stress can be induced by a nearby, new child well is hydraulically pressured. Unlike the pore pressure change itself (which is a scalar quantity), the dilation of fracture planes causes elastic stresses that superimpose on a pre-existing stress state (which is a tensor quantity). Stress reversals, therefore, are not a primary, but rather a secondary effect of poro-elasticity, as these occur due to certain anisotropic "far-field" and "near-field" stress superposition situations. Examples will be shown in our study.
The ultimate goal of all modeling efforts on the reorientation of hydraulic fractures due to well interference, is to next plan the parent and infill-well placement in the field such that any stress reversals can be controlled to be conducive for fracture network development leading to production optimization [16,25]. If fracture hits between parent and chilled wells are to be avoided, insight from fracture development simulators are helpful [1,2]. The reorientation of hydraulic fractures and the occurrence of local stress reversals are suggested to be induced by concurring factors such as the propagation of hydraulic fractures, the presence of natural fractures that create local stress re-distributions interfering with the pre-existing far-field stresses.

Accounting for Poro-Elasticity in LSM and TLSM Models
Poro-elasticity is accounted for in this study by simple addition of the pore pressure multiplied by the Biot's constant to the directional stress magnitude. The effect of pore pressure is shown to induce a change in the stress magnitude imposed by the hydraulic fractures on the elastic medium. Fundamental equations used in the Linear Superposition Method (LSM) and the Time-Stepped Linear Superposition Method (TLSM) are briefly outlined, and modified to model stress interference with poro-elasticity, first in this section. A geomechanical sign convention was employed to describe the magnitude of stresses for all model results. Henceforth, compressional stress magnitudes are positive, while tensional stress magnitudes are negative. interactions between pressurized fractures [15] and boreholes [26]. TLSM and LSM are proposed as fast, gridless alternatives to commonly used, grid-dependent fracture propagation models, which are computationally more demanding and henceforth more expensive.
LSM was first introduced to model stress accumulations induced by boundary displacements in an elastic medium perforated by pressurized boreholes [26]. Then, LSM was expanded to model the stress interference between static hydraulic fractures and compared to independent photo-elastic results to evaluate physical compatibility [15]. By assuming perforation tunnels as a set of pressurized fractures from the wellbore, the incorporation of LSM solves for the stress distribution and the potential fracture propagation paths from the perforations in a given wellbore [17]. The results from the combination of fractures and boreholes solved by LSM indicate that the maximum compressive stress occurs at the rim of the wellbore, while the maximum tensile stress occurs at the perforation tips. Additionally, the presence of the pressurized fractures in the near vicinity of the wellbore deviates the principal stresses at the wellbore wall. Wellbore instability may occur with a tensile fracture or by shear failure.
TLSM was introduced as an expansion of LSM with a Eulerian time-stepped component to model dynamic fracture propagation to give insights into Hydraulic Fracture Test Site (HFTS) seismic data from the Wolfcamp Formation in the Permian Basin (Texas) [8]. The results indicate that the dispersed appearance seismic clouds can be explained by hydraulic fracture divergence during the fracturing process. Additionally, TLSM was applied to observe and analyze the dynamic progression of hydraulic fracture hits and stress trajectory reorientation in typical single and multi-well completions [16].
The results indicate that the number of fracture hits and degree of redirection are lower for a typical modified zipper fracturing completion. Additionally, fracture hits and redirection are less prominent when perforation spacing is higher, and when well spacing is tighter.

Basic LSM Equations
The summation of displacement fields imposed by pressurized fractures is employed to calculate the stress magnitudes induced in the elastic medium. The coordinate system of a single fracture modified for LSM and TLSM is shown in Figure 1.
Energies 2020, 13, x FOR PEER REVIEW 4 of 22 mechanical interactions between pressurized fractures [15] and boreholes [26]. TLSM and LSM are proposed as fast, gridless alternatives to commonly used, grid-dependent fracture propagation models, which are computationally more demanding and henceforth more expensive. LSM was first introduced to model stress accumulations induced by boundary displacements in an elastic medium perforated by pressurized boreholes [26]. Then, LSM was expanded to model the stress interference between static hydraulic fractures and compared to independent photoelastic results to evaluate physical compatibility [15]. By assuming perforation tunnels as a set of pressurized fractures from the wellbore, the incorporation of LSM solves for the stress distribution and the potential fracture propagation paths from the perforations in a given wellbore [17]. The results from the combination of fractures and boreholes solved by LSM indicate that the maximum compressive stress occurs at the rim of the wellbore, while the maximum tensile stress occurs at the perforation tips. Additionally, the presence of the pressurized fractures in the near vicinity of the wellbore deviates the principal stresses at the wellbore wall. Wellbore instability may occur with a tensile fracture or by shear failure.
TLSM was introduced as an expansion of LSM with a Eulerian time-stepped component to model dynamic fracture propagation to give insights into Hydraulic Fracture Test Site (HFTS) seismic data from the Wolfcamp Formation in the Permian Basin (Texas) [8]. The results indicate that the dispersed appearance seismic clouds can be explained by hydraulic fracture divergence during the fracturing process. Additionally, TLSM was applied to observe and analyze the dynamic progression of hydraulic fracture hits and stress trajectory reorientation in typical single and multiwell completions [16]. The results indicate that the number of fracture hits and degree of redirection are lower for a typical modified zipper fracturing completion. Additionally, fracture hits and redirection are less prominent when perforation spacing is higher, and when well spacing is tighter.

Basic LSM Equations
The summation of displacement fields imposed by pressurized fractures is employed to calculate the stress magnitudes induced in the elastic medium. The coordinate system of a single fracture modified for LSM and TLSM is shown in Figure 1. (1) The elastic displacements, induced by a single hydraulic fracture, in the x and y direction are ( Figure 1): Energies 2020, 13, 6474 where P is the internal pressure of the hydraulic fracture, v is the bulk Poisson's ratio, and E is the bulk Young's modulus. The strain tensor elements are calculated by taking the partial derivative of the displacements in the x and y direction: Directional stresses in the x and y direction are determined from the strain tensor elements (ε xx , ε yy , ε xy ). The constitutive equation was used as the assumed boundary condition in the original Sneddon solution: The hydrostatic pressure in the bulk rock is assumed to remain isotropic and scalar throughout the poro-elastic medium, and relates linearly to the pore pressure P P multiplied by a scaling factor α (Biot's constant): The pore pressure can be assumed to either be constant throughout the reservoir, for the sake of simplicity, or may spatially vary, as shown later in this study.
Combining Equation (7) with Equations (5) and (6) gives the total stress: The shear stress is not affected by the pore pressure.
Stress tensor components (σ xx , σ yy , σ xy ) or (σ xx , σ yy , σ xy ) are used to obtain the magnitude of the principal stresses, σ 1 and σ 2 , for all coordinate within the observation plane by applying the standard expressions: The maximum principal stress trajectory can be calculated with: According to Equation (13), the principal stress trajectories are not affected by any pore pressure offsets, because the term σ xx − σ yy annuls the poro-elastic terms appearing in Equations (8) and (9). For multiple cracks, the elastic displacements, due to all the superposed hydraulic fractures, are summed directionally: After calculating the total displacement field for multiple cracks, Equations (3)-(13) are applied to solve the strain, stress, and stress trajectories. This methodology was used to model the stress field contours resulting from two active hydraulic fractures, accounting for the effect of poro-elasticity, as shown in Figure 2.
For multiple cracks, the elastic displacements, due to all the superposed hydraulic fractures, are summed directionally: After calculating the total displacement field for multiple cracks, Equations (3)-(13) are applied to solve the strain, stress, and stress trajectories. This methodology was used to model the stress field contours resulting from two active hydraulic fractures, accounting for the effect of poroelasticity, as shown in Figure 2.

Static LSM Model Results
The static Linear Superposition Method (LSM) is applied here to study the spatial effect of poro-elasticity. The linear elasticity assumption and plane strain boundary condition result in a 2D displacement model and a 3D stress model. Additionally, all cracks are assumed to be cylindrical. These assumptions ensure that all 2D solutions are identical for each parallel slice in 3D space. Simulation inputs for an LSM model with two active hydraulic fractures are given in Table 1.

Static LSM Model Results
The static Linear Superposition Method (LSM) is applied here to study the spatial effect of poro-elasticity. The linear elasticity assumption and plane strain boundary condition result in a 2D displacement model and a 3D stress model. Additionally, all cracks are assumed to be cylindrical. These assumptions ensure that all 2D solutions are identical for each parallel slice in 3D space. Simulation inputs for an LSM model with two active hydraulic fractures are given in Table 1. To analyze the influence of poro-elasticity on numerical displacement and stress value, the values of displacement and stress magnitudes are compared at a single coordinate (when x = 0 and y = 0) in static LSM models with Pp = 0 ( Figure 2a) and Pp = 4439 ( Figure 2b). Table 2 shows the numerical Energies 2020, 13, 6474 7 of 22 comparison of the principal stresses spatially. Pore pressure is shown to elevate the magnitude of the two principal stresses, at the origin of the coordinate system (x,y) = (0,0). The local pressure increase creates a stress difference as shown in Table 2, exactly equal to pore pressure (3995 psi) multiplied by Biot's constant.  Table 3 compares the magnitude of the principal stresses at another arbitrary point (x,y) = (10,10) to further observe the spatial effect of poro-elasticity. The principal stress magnitudes are entirely different from values previously recorded at (x,y) = (0,0). However, the results in Table 1 show that the pore pressure again increases the magnitude of both principal stresses, and the stress differences in Table 3 are exactly the same (3995 psi) as in Table 2. Next, assume that there exists a spatial gradient in the pore pressure, within the elastic medium of interest, with a maximum scalar gradient along the x-direction, ranging from 1000 to 4439 psi ( Figure 3). Such a spatial gradient would occur in well models (map view) with a nearby legacy well at the left boundary that has depleted the reservoir pressure to 1000 psi (left boundary), while at the right hand side of the lease boundary the pristine reservoir pressure still exists (4439 psi). The local stress field will change when a horizontal child well is drilled and completed in the middle of the lease region. The contours for the spatial pore pressure, in the absence of any fractures, are shown in Figure 3.  When two hydraulic fractures are introduced ( Figure 2) with the condition as in Table 1, except for pore pressure varying spatially as in Figure 3, then the spatial variation in the maximum principal stress magnitude occurs, as in Figure 4. The stress contours are no longer observed to be symmetric about the wellbore at x = 0, due to the effect of variable pore pressure in the model space. In effect, the left side of Figure 4 reflects Figure 2a, and the right-hand side of Figure 4 reflects Figure 2b. Additionally, the principal stress trajectory remains symmetric for the left-and righthand side, as shown in Figure 4b. When two hydraulic fractures are introduced ( Figure 2) with the condition as in Table 1, except for pore pressure varying spatially as in Figure 3, then the spatial variation in the maximum principal stress magnitude occurs, as in Figure 4. The stress contours are no longer observed to be symmetric about  When two hydraulic fractures are introduced (Figure 2) with the condition as in Table 1, except for pore pressure varying spatially as in Figure 3, then the spatial variation in the maximum principal stress magnitude occurs, as in Figure 4. The stress contours are no longer observed to be symmetric about the wellbore at x = 0, due to the effect of variable pore pressure in the model space. In effect, the left side of Figure   In order to analyze the effect of variable pore pressure in the poro-elastic LSM model on the numerical displacement and stress values are compared at laterally opposite coordinate points in the model of Figure 4 ( Table 4).  In order to analyze the effect of variable pore pressure in the poro-elastic LSM model on the numerical displacement and stress values are compared at laterally opposite coordinate points in the model of Figure 4 ( Table 4).  Tables 2 and 3, the escalation of pore pressure leads to an increase in principal stress magnitudes. The principal stress magnitude difference is exactly the same at both coordinates (−20, 20) and (20,20). The principal stress magnitude difference spatially is observed to be uniform and constant throughout the elastic medium. Additionally, the principal stress magnitude difference between two coordinates is equal to that of the pore pressure difference of the same coordinates multiplied by the Biot's constant: The stress magnitudes induced into the poro-elastic medium by the hydraulic fractures are increased when the pore pressure increases. In the examples given, the stress increases in a linear fashion in every spatial direction towards higher (positive) x-coordinates. In the y-direction, the pore pressure gradient vanishes, and differential stress magnitudes are reduced everywhere along a certain y-coordinate by the same amount, equal to the local pore pressure multiplied by the Biot's constant.

Sub-Conclusions
Several sub-conclusions can be inferred from the poro-elastic models presented in this section:

1.
Principal stress trajectories are not affected by the regional pore pressure variations (Figure 4b).
Energies 2020, 13, 6474 9 of 22 2. However, pore pressure changes have a direct relationship with principal stress magnitudes-as pore pressure increases, the principal stress magnitudes will increase (Figure 4a).

3.
Changes in the principal stress magnitude, induced by the pore pressure, were constant for all spatial coordinates in the elastic medium. 4.
Induced changes in the principal stress magnitudes are equal to the pore pressure multiplied by the Biot's constant.

Advancing TLSM with Time-Dependent Pressure
In this section, the Time-Stepped Linear Superposition Method (TLSM) is integrated with field Diagnostic Fracture Injection Test data. This integration of data into the TLSM model allows for a more realistic interpretation of the hydraulic fracturing process, accounting for the time-dependent changes in the pressure load of individual fractures.

Time-Dependent Pressure Equations and Assumptions
The present TLSM model adopts several assumptions in order to simplify the fracture propagation problem: (1) all hydraulic fractures are assumed to have the same storage capacity and propagate uniformly to contain all injection fluid, (2) linear elasticity, and (3) plane strain. The stress concentration induced by each hydraulic fracture at any point in the elastic medium can be described by the summation of the displacements in each direction. The directional strain tensor solutions are determined using a partial derivative of the spatial displacement field. Directional stress magnitudes are then calculated from these strain tensors. The maximum principal stress trajectory calculated at the fracture tip determines the direction of future fracture growth.
A simplified failure criterion is used such that if the tensile principal stress at the tip of the fracture exceeds that of the rock strength, then the fracture propagates. Conversely, when stresses at the fracture tip due to the pressure load are no longer critical, then the fracture will no longer propagate. Following mechanical failure, the propagation path of the hydraulic fracture is calculated with a Eulerian displacement scheme. The successive predicted future location z j of each fracture tip is [8,16]: where z is the fracture tip position in 2D space, t is the time in seconds, V is the fracture velocity in ft/s, and ∆t is the time step in seconds. Using Equation (17), the fracture growth velocity (V) accounts for any changes in pump rates in the DFIT data. Additionally, the timestep parameter (∆t), is used to account for events happening in the DFIT data such as the start of injection, the duration of injection, and shut-in time.

Dynamic Fracture Propagation
TLSM models of fracture propagation use the changes in the global displacement field at each time-step to adjust the incremental superstition of successive displacement fields. Each time-step accounts for the concurrent fluid volume injected, and the elastic deformation resulting from such enacted force. Static equilibrium is assumed at each time-step, which means the sum of all forces in each direction is constrained to zero at every predetermined time-step. Additionally, no linear or angular accelerations are permitted to occur. This implies that at every successive time-step, the hydraulic fracture in equilibrium can propagate, but its linear and angular velocities remain constant. Similarly to previously used, common fracture propagation frameworks, there is no acceleration involved in the systems modeled. Thus, the inertia term in the equation of motion is neglected.
The justification for neglecting the inertia terms is that any stress perturbations involving inertia would be very small or absent in conventional dynamic simulations of fracture propagation. In other words, such approaches could be called quasi-static in the strict sense that dynamic terms are assumed to remain negligible and thus are neglected altogether. Such simplification may not be warranted, for example, if a fracture propagation case has an abrupt initiation or termination, such that dynamic stress perturbations associated with the occurrence of seismic waves may become significant. Seismicity may cause pressure fronts to run ahead of the stress-tip-induced stress perturbations, and would potentially lower the threshold for failure at the fracture tip, due to dynamically superposed stress field changes resulting from inertia effects [27].

Diagnostic Fracture Injection Tests (DFIT)
In unconventional exploration, diagnostic fracture injection tests (DFIT) have become a primary method to determine key reservoir parameters-for later use in hydraulic fracture treatment operations-such as the fluid efficiency, minimum in situ stress, pore pressure, leak-off coefficient, and fracture gradient. DFITs are typically completed in short time-frames spanning a mere couple of days, and commonly less than two weeks, which is significantly shorter as compared to the traditional pressure build-up tests, which could last several weeks due to the requirement of stable-flowing conditions, before the start of the shut-in period [28]. DFITs involve the injection of a relatively small fluid volume, typically at a first stage in the toe of a newly drilled horizontal well, with a constant injection rate to obtain information about the hydraulic fracture treatment parameters as well as the reservoir. Due to the risk of formation damage, the fluid injection rate is commonly limited to a constant rate of 5-6 bbl/min over a very brief period of 3-5 min [29].
The well is shut-in following the DFIT injection, and the pressure fall-off is recorded over several days-either at the surface or at the bottom-hole level. When the reservoir pressure is equal to or greater than the hydrostatic pressure in the wellbore, the recording of surface pressure response data suffices [30]. During the shut-in period, injection fluid gradually leaks off into the reservoir, causing the hydraulic mini-fracture created by the DFIT to close. Using the data obtained during the pressure de-escalation period, after shut-in, one can determine the fracture closure stress and pore pressure [31].
Traditionally, DFIT analysis methods use the G-function model and the tangent line method to evaluate the fracture closure stress. The G-function represents the shut-in time [32]. Recent developments in DFIT analysis apply forward models and inverse calculations to determine the stress-dependent, unpropped fracture conductivity [33].
Data from DFITs can be used to enhance the realistic representation of fracture propagation models. For example, DFIT-induced fractures in a pressure depleted environment will have a larger footprint near the depleted region [34]. Additionally, the growth of hydraulic fractures during injection and shut-in was observed to be episodic, where the fractures grow intermittently through phases instead of in a steady continuous rate [35]. The presence of any natural fractures has also been shown to influence the responses of hydraulic fractures during DFITs. The pumping schedule and the hydraulic fracture propagation path may impact the DFIT estimation of the minimum horizontal stress [36].

Field DFIT Data and Simulation Parameters
Field data from a Diagnostic Fracture Injection Test are shown in Figure 5. The breakdown pressure, instantaneous shut-in pressure (ISIP), fracture propagation pressure, fracture closure pressure, and injection cycles are labeled. During this injection cycle, the injection rate remains constant while the reservoir pressure is monitored and recorded. The TLSM model applied here assumes the existence of failure nuclei in the elastic medium and that breakdown pressure already occurred. Hence, the second injection cycle is used for hydraulic fracture propagation. The associated data with the start time of the second injection cycle normalized to 0 are shown in Figure 6. constant while the reservoir pressure is monitored and recorded. The TLSM model applied here assumes the existence of failure nuclei in the elastic medium and that breakdown pressure already occurred. Hence, the second injection cycle is used for hydraulic fracture propagation. The associated data with the start time of the second injection cycle normalized to 0 are shown in Figure  6.  The starting pressure recorded in the analyzed data (of Figures 5 and 6) is 35.5 psia; hence, the pressure gauge is assumed to be positioned at the surface. Neglecting pipe friction, the bottom hole tubing pressure (BHTP) is calculated by: where Ps is the recorded surface pressure of 35.5 psia and Pp is the calculated constant pore pressure of 4439 psia (assumed pore pressure gradient of 0.6 psi/ft). The bottom hole pressure vs. time curve is calculated using Equation (18) and shown in Figure 7. The fracture propagation pressure inferred from the bottom hole pressure vs. time curve in Figure 5 is 9370.83 psi and was integrated into the TLSM model. This fracture propagation pressure is constrained to be constant along with the fluid injection rate during the fracture propagation process. Further DFIT pressure and time parameters used in TLSM are given in Table 5. A typical constant injection rate of 6 bpm/min for a DFIT test was assumed to simulate the fracture propagation process due to fluid injection [29]. This allows for the calculation of total injected volume, and the number of fractures to contain those fluid volumes with assumed fracture geometry. Further fracture geometric and volumetric parameters used in TLSM are given in Table 6.  The starting pressure recorded in the analyzed data (of Figures 5 and 6) is 35.5 psia; hence, the pressure gauge is assumed to be positioned at the surface. Neglecting pipe friction, the bottom hole tubing pressure (BHTP) is calculated by: where P s is the recorded surface pressure of 35.5 psia and P p is the calculated constant pore pressure of 4439 psia (assumed pore pressure gradient of 0.6 psi/ft). The bottom hole pressure vs. time curve is calculated using Equation (18) and shown in Figure 7. The fracture propagation pressure inferred from the bottom hole pressure vs. time curve in Figure 5 is 9370.83 psi and was integrated into the TLSM model. This fracture propagation pressure is constrained to be constant along with the fluid injection rate during the fracture propagation process. Further DFIT pressure and time parameters used in TLSM are given in Table 5. A typical constant injection rate of 6 bpm/min for a DFIT test was assumed to simulate the fracture propagation process due to fluid injection [29]. This allows for the calculation of total injected volume, and the number of fractures to contain those fluid volumes with assumed fracture geometry. Further fracture geometric and volumetric parameters used in TLSM are given in Table 6.
where Ps is the recorded surface pressure of 35.5 psia and Pp is the calculated constant pore pressure of 4439 psia (assumed pore pressure gradient of 0.6 psi/ft). The bottom hole pressure vs. time curve is calculated using Equation (18) and shown in Figure 7. The fracture propagation pressure inferred from the bottom hole pressure vs. time curve in Figure 5 is 9370.83 psi and was integrated into the TLSM model. This fracture propagation pressure is constrained to be constant along with the fluid injection rate during the fracture propagation process. Further DFIT pressure and time parameters used in TLSM are given in Table 5. A typical constant injection rate of 6 bpm/min for a DFIT test was assumed to simulate the fracture propagation process due to fluid injection [29]. This allows for the calculation of total injected volume, and the number of fractures to contain those fluid volumes with assumed fracture geometry. Further fracture geometric and volumetric parameters used in TLSM are given in Table 6.    cycle ends, and the hydraulic fractures cease to propagate. When hydraulic fracture propagation ceases at the end of the injection cycle, fracture pressure de-escalation begins due to leak-off. Figure 8 shows the fracture pressure vs. time curve used in the TLSM model, which was set to 0 psi at time 0 to represent an undisturbed and unfractured formation. The fracture pressure was then increased to 9371 psi and remained constant for 4.58 min. At the 4.58 min mark in Figure 8, the injection cycle ends, and the hydraulic fractures cease to propagate. When hydraulic fracture propagation ceases at the end of the injection cycle, fracture pressure de-escalation begins due to leak-off.  To investigate the spatial effect of induced stress by hydraulic fracturing, Figure 9 shows the maximum principal stress recorded at (x,y) = (0,0) vs. time curve with the maximum principal stress at a distance of 25 ft away from the perforations, using with the same fracture pressure inputs as in Figure 8. The maximum principal stress at (x,y) = (0,0) is observed to increase rapidly at the start of injection and peaked at approximately 11.22 s after the injection cycle begins. However, the maximum principal stress begins to decrease exponentially even before the end of injection. This decline in maximum principal stress is continued until the pore pressure is reached. The signature plot trajectory is particularly similar to that of the DFIT data, where the pressure increases dramatically as the fracture propagates, exponentially drops, and plateaus as the fracture closes. To investigate the spatial effect of induced stress by hydraulic fracturing, Figure 9 shows the maximum principal stress recorded at (x,y) = (0,0) vs. time curve with the maximum principal stress at a distance of 25 ft away from the perforations, using with the same fracture pressure inputs as in Figure 8. The maximum principal stress at (x,y) = (0,0) is observed to increase rapidly at the start of injection and peaked at approximately 11.22 s after the injection cycle begins. However, the maximum principal stress begins to decrease exponentially even before the end of injection. This decline in maximum principal stress is continued until the pore pressure is reached. The signature plot trajectory is particularly similar to that of the DFIT data, where the pressure increases dramatically as the fracture propagates, exponentially drops, and plateaus as the fracture closes.

Sub-Conclusions
Several sub-conclusions can be formed based on integrated DFIT data with the TLSM model in this section: 1. The maximum principal stress at (x,y) = (0,0) during injection increases dramatically as the fracture propagates, exponentially drops, and steadies on a plateau value as the fracture closes, while injection pressure remained constant (Figure 9).

Sub-Conclusions
Several sub-conclusions can be formed based on integrated DFIT data with the TLSM model in this section: 1.
The maximum principal stress at (x,y) = (0,0) during injection increases dramatically as the fracture propagates, exponentially drops, and steadies on a plateau value as the fracture closes, while injection pressure remained constant (Figure 9). 2.
Hence, the variation in maximum principal stress is primarily caused by the changes in fracture geometry as the fracture propagates, rather than by changes in injection pressure.

TLSM Results with Imposed Pore Pressure Gradient
In this section, the dynamic propagation of two hydraulic fractures is modeled with a regional pore pressure gradient with and without a near-field stress superposition due to a pre-existing natural fracture. While the pore pressure gradient alone does not lead to any discernable reorientation of propagating hydraulic fractures (Section 5.1), when a near-field stress perturbation exists (in this case a pressurized natural fracture, Section 5.2), such reorientation is very pervasive. Figure 10a,b show, respectively, that the scalar pore pressure magnitude and stress trajectory isoclines do not exist in the initial poro-elastic medium before the onset of the hydraulic fracture treatment. The initial pore pressure varies spatially, but no stress tensor field exists (Figure 10a). The principal stress trajectory isoclines show no spatial variation and are zero everywhere, since hydraulic fractures have not yet been activated to impose a stress tensor field.  Figure 11 shows the TLSM results for the first hydraulic fracture nuclei imposed on the 2D poro-elastic medium used in TLSM. The presence of the pore pressure gradient does impact the maximum principal stress magnitudes but does not affect the principal stress trajectories isoclines. The extent of the altered maximum principal stress is limited to that of the area around the hydraulic fracture, while the principal stress trajectories now appear everywhere in the poro-elastic medium due to the distortion of the scalar pressure field by the stress tensor field induced by the fractures.  Figure 11 shows the TLSM results for the first hydraulic fracture nuclei imposed on the 2D poro-elastic medium used in TLSM. The presence of the pore pressure gradient does impact the maximum principal stress magnitudes but does not affect the principal stress trajectories isoclines. The extent of the altered maximum principal stress is limited to that of the area around the hydraulic fracture, while the principal stress trajectories now appear everywhere in the poro-elastic medium due to the distortion of the scalar pressure field by the stress tensor field induced by the fractures.

Hydraulic Fracture Propagation-Pre-Existing Pore Pressure Gradient
poro-elastic medium used in TLSM. The presence of the pore pressure gradient does impact the maximum principal stress magnitudes but does not affect the principal stress trajectories isoclines. The extent of the altered maximum principal stress is limited to that of the area around the hydraulic fracture, while the principal stress trajectories now appear everywhere in the poro-elastic medium due to the distortion of the scalar pressure field by the stress tensor field induced by the fractures.
(a) (b) Figure 11. TLSM model results for the initial hydraulic fracture nuclei showing (a) maximum principal stress. Color code is shifted with respect to Figure 10a, due to stress peaks at the fracture nuclei. (b) Principal stress trajectory isoclines at t = 0.093 min. TLSM results were generated in 0.466 s including post-processing. Figure 12 shows the next successive hydraulic fracture propagation step after Figure 11. The maximum principal stress and the maximum principal stress trajectory isoclines are observed to be further changed as the fractures propagate. The maximum principal stress contours increase uniformly from the left to the right boundary of the study region, with the pressure depletion attributed to a Type B intervention (see Section 2). Alterations in the magnitude of the stress occur in the vicinity of the imposed hydraulic fractures, and the principal stress trajectories are also altered systematically throughout the medium.  Figure 12 shows the next successive hydraulic fracture propagation step after Figure 11. The maximum principal stress and the maximum principal stress trajectory isoclines are observed to be further changed as the fractures propagate. The maximum principal stress contours increase uniformly from the left to the right boundary of the study region, with the pressure depletion attributed to a Type B intervention (see Section 2). Alterations in the magnitude of the stress occur in the vicinity of the imposed hydraulic fractures, and the principal stress trajectories are also altered systematically throughout the medium.  Figure 13 shows the TLSM results for the final hydraulic fractures at the end of the injection period, immediately prior to the de-escalation of fracture pressure due to shut-in. The hydraulic fractures are observed to diverge uniformly through region of variable pore pressure. The stress induced by the propagated hydraulic fracture at this stage is dominant over the existing pore pressure in the reservoir.  Figure 13 shows the TLSM results for the final hydraulic fractures at the end of the injection period, immediately prior to the de-escalation of fracture pressure due to shut-in. The hydraulic fractures are observed to diverge uniformly through region of variable pore pressure. The stress induced by the propagated hydraulic fracture at this stage is dominant over the existing pore pressure in the reservoir. generated in 1.151 s including post-processing. Figure 13 shows the TLSM results for the final hydraulic fractures at the end of the injection period, immediately prior to the de-escalation of fracture pressure due to shut-in. The hydraulic fractures are observed to diverge uniformly through region of variable pore pressure. The stress induced by the propagated hydraulic fracture at this stage is dominant over the existing pore pressure in the reservoir.   Figure 13). This is due to the neglection of fracture closure in the current TLSM model. In more complex scenarios, the closure of hydraulic fractures would reduce fracture width and fracture length, and thus would change the induced stress at every point in the model. Similar to the results in Figure 13, the propagated hydraulic fracture remains the dominant stress inducer relative to the existing pore pressure. As the pressure inside the fracture decreases due to leak-off, so does the maximum principal stress.  The principal stress trajectory isoclines during leak-off remained unchanged as compared to the result at the end of injection ( Figure 13). This is due to the neglection of fracture closure in the current TLSM model. In more complex scenarios, the closure of hydraulic fractures would reduce fracture width and fracture length, and thus would change the induced stress at every point in the model. Similar to the results in Figure 13, the propagated hydraulic fracture remains the dominant stress inducer relative to the existing pore pressure. As the pressure inside the fracture decreases due to leak-off, so does the maximum principal stress.  Figure S1a shows the dynamic TLSM result for two hydraulic fractures in a region with pore pressure gradient showing maximum principal stress, as shown in Figures 11a-14a. Figure S1b shows the dynamic TLSM result for two hydraulic fractures in a region with pore pressure gradient showing principal stress isoclines, as shown in Figures 11b-14b.

Hydraulic Fracture Propagation-Uniform Pressure Gradient and a Pre-Existing Near-Field Stress
In Section 5.1, the pore-pressure gradient did not appear to modify hydraulic fracture propagation direction. In this section, a pre-existing pressurized natural fracture is assumed, in addition to the pore-pressure gradient. The natural fracture was imposed at a position of x = −25 feet and y = 100 ft, with a half-length of 100 ft, and oriented −30° from the positive x-axis. The average pressure inside the fault was set to 1000 psi.
The far-field stress imposed by the natural fracture changed the propagation direction of the two hydraulic fractures. The final hydraulic fracture propagation patterns at the end of injection period are shown in Figure 15, and after leak-off shown in Figure 16. The final fracture pattern no longer resembles that of Figures 13 and 14. Additionally, the hydraulic fractures propagated to   Figure S1a,b. Figure S1a shows the dynamic TLSM result for two hydraulic fractures in a region with pore pressure gradient showing maximum principal stress, as shown in Figure 11a, Figure 12a, Figure 13a, Figure 14a. Figure S1b shows the dynamic TLSM result for two hydraulic fractures in a region with pore pressure gradient showing principal stress isoclines, as shown in Figure 11b, Figure 12b, Figure 13b, Figure 14b.

Hydraulic Fracture Propagation-Uniform Pressure Gradient and a Pre-Existing Near-Field Stress
In Section 5.1, the pore-pressure gradient did not appear to modify hydraulic fracture propagation direction. In this section, a pre-existing pressurized natural fracture is assumed, in addition to the pore-pressure gradient. The natural fracture was imposed at a position of x = −25 feet and y = 100 ft, with a half-length of 100 ft, and oriented −30 • from the positive x-axis. The average pressure inside the fault was set to 1000 psi.
The far-field stress imposed by the natural fracture changed the propagation direction of the two hydraulic fractures. The final hydraulic fracture propagation patterns at the end of injection period are shown in Figure 15, and after leak-off shown in Figure 16. The final fracture pattern no longer resembles that of Figures 13 and 14. Additionally, the hydraulic fractures propagated to develop mutual fracture hits, shown by their collision, as well as with the natural fracture. Stress reorientation induced by a near-field stress superposition due to a nearby natural fracture may profoundly affect the hydraulic fracture propagation paths. The pressure depletion due to Type B interventions may not be a primary cause for stress reorientation, but Type A interventions may run into local stress perturbations that alter the fracture path. Previous works attributed fracture reorientation and stress reversals to pore pressure changes by reservoir pressure depletion due to production [12,13]. However, pressure depletion in wells would not cause deviatoric stress changes, thus will not change the preferred fracture paths in the reservoir.

Sub-Conclusions
The model results generated by TLSM for the competing effects of poro-elasticity and far-field stress gives rise to the following sub-conclusions: 1. The spatial variation of pore pressure does not affect the fracture propagation direction ( Figure  13). 2. Changes in stress state are concentrated in areas closest to the hydraulic fractures. 3. The superposition of a near-field stress, in this case a pressurized natural fracture, causes stress reversals, as well as the reorientation of what previously would be symmetric hydraulic fracture patterns.

Sub-Conclusions
The model results generated by TLSM for the competing effects of poro-elasticity and far-field stress gives rise to the following sub-conclusions: 1. The spatial variation of pore pressure does not affect the fracture propagation direction ( Figure  13). 2. Changes in stress state are concentrated in areas closest to the hydraulic fractures. 3. The superposition of a near-field stress, in this case a pressurized natural fracture, causes stress reversals, as well as the reorientation of what previously would be symmetric hydraulic fracture patterns.  Figure S2a,b. Figure S2a shows the dynamic TLSM result for two hydraulic fractures in proximity to a pressurized natural fracture showing maximum principal stress, as shown in Figure 15a, Figure 16a. Figure S2b shows the dynamic TLSM result for two hydraulic fractures in proximity to a pressurized natural fracture showing principal stress trajectory isoclines, as shown in Figure 15b, Figure 16b.

Sub-Conclusions
The model results generated by TLSM for the competing effects of poro-elasticity and far-field stress gives rise to the following sub-conclusions: 1.
The spatial variation of pore pressure does not affect the fracture propagation direction ( Figure 13).

2.
Changes in stress state are concentrated in areas closest to the hydraulic fractures. 3.
The superposition of a near-field stress, in this case a pressurized natural fracture, causes stress reversals, as well as the reorientation of what previously would be symmetric hydraulic fracture patterns.

Pore Pressure and Principal Stress Reorientation
Prior suggestions that principal stress re-orientations would occur due to pore pressure depletion [12,37] (i.e., Type B interventions) seem untenable, as follows from the model results in this study. Stress reversals are not a primary effect of gradients in pore pressure, but occur due to certain far-field and near-field stress superposition situations.
In practice, local changes in near-field stress can be induced by Type A interventions when a nearby, new child well is hydraulically pressured [38], which then causes a deviatoric stress field to arise in the vicinity of the well. Alternatively, a pre-existing natural fracture may also cause a change in the so-called near-field stress. The reorientation of hydraulic fractures observed in high and low pore pressure areas may not be an artifact of spatial variation of pore pressure, but is instead due to spatial variations of the shadowed stresses.

Strengths and Weaknesses
Similar to other hydraulic fracture propagation modeling methods, there are strengths and weaknesses associated with the use of TLSM. Among the strengths of TLSM are that it offers grid-less, closed-form, semi-analytical solutions for non-planar fracture propagation. By avoiding adaptive grid-refinement, fast simulation of multiple competitive hydraulic fracture propagation and infinite resolution are achieved. Among the weaknesses of both LSM and TLSM is that the current codes are constrained to a linear elasticity assumption.
The present study adapted the code to account for linear poro-elasticity and also accommodates field-based pressure inputs (DFIT). The growth of hydraulic fractures is constrained to account for all the fluid volumes injected at each timestep. Hence, the hydraulic fractures must propagate in order to contain the volumes injected. Since all hydraulic fractures are assumed to have the same storage capacity and a positive propagation velocity, fracture closure and the reduction in hydraulic fracture length were neglected.

Future Work
This study focused on the modeling of fracture propagation with TLSM, assuming internally pressurized fractures with the same storage capacity to accommodate all injection fluid. Due to fracture closure during leak-off, reduction in the conductive length of the hydraulic fracture occurs [38]. Instead of the current constant volume injection scheme used in the current code, a volume-dependent fracture geometry scheme may be implemented. In future studies, the TLSM code can be modified to account for injection volume-dependent fracture propagation of each fracture. The fluid leak-off process can be modeled by coupling the TLSM fracture propagation simulator with our semi-analytical simulator for fluid flow in fractured porous media based on Complex Analysis Methods (CAM). Key