Stick-slip like behavior in shear fracture propagation including the effect of fluid flow

Shear-based fracture propagation in fluid-saturated porous materials is investigated using a displacement–pressure formulation that includes acceleration and inertial effects of the fluid. Pressure-dependent plasticity with a nonassociated flow rule is adopted to realistically represent the stresses in the porous bulk material. The domain is discretized using unequal order T-splines and cast into a finite element method using Bézier extraction. An implicit scheme is used for the temporal integration. The solid acceleration-driven fluid flow reacts to stress waves, but it results in pressure oscillations. Adding fluid acceleration terms dampens these oscillations and increases the fluid pressure near the fracture tips. By simulating a typical shear fracture case, it is shown that stick-slip like, or stepwise, fracture propagation occurs for a high permeability, also upon mesh refinement. The acceleration driven fluid flow results in a build-up of pressure near the fracture tip. Once this pressure region encompasses the fracture tip, propagation arrests until the pressure has diffused away from the crack tip, after which propagation is resumed and the build-up of pressure begins anew. This results in a stick-slip like behavior, with large arrests in the fracture propagation. Stepwise propagation related to the initial conditions has also been observed, but disap-pears once the fracture length exceeds the size of the region influenced by the initial conditions.

Simulations on shear-based fracture propagation are typically combined with off-fault plasticity with a nonassociated flow rule, either assuming drained or undrained behavior for the porous bulk material. [10][11][12][13] These simulations show that inclusion of the fluid has a large impact on the plastic strains. However, using a "drained" assumption neglects the pressure changes due to compression and the inertial effects of the fluid required to distribute local pressure increases. While the undrained assumption gives a better representation in the low-permeability limit, it neglects fluid diffusion. Indeed, allowing fluid diffusion has been shown to alter the plastic strains significantly, 14 and simulations without plasticity have shown the influence of the fluid diffusion on the fracture propagation. 15 Simulations without plasticity and accelerationdriven flow have shown that the interstitial fluid pressure can potentially stop the fracture from propagating. 16 The acceleration of the fluid is often neglected in the expression for the fluid velocity in Darcy's law. [17][18][19] This is a reasonable assumption for low-frequency responses, but it is less clear whether this assumption holds close to a propagating fracture considering the accelerations that occur around the crack tip. Simulations have shown that the presence of these accelerations invalidates the assumption of a low-frequency response and that neglecting the fluid acceleration terms results in different solutions. 20 Additionally, the fluid acceleration term in Darcy's law damps the pressure and stress waves, [21][22][23] and therefore, results in different values for the stresses near the fracture tip, which has been shown using a simplified formulation. 24 Also when simulating porous materials with a relatively high porosity and permeability, neglecting this term causes differences in the solution. 25 Recently, several authors have argued that the propagation of fractures in porous materials is not smooth, but stepwise. [26][27][28][29] This stepwise propagation is observed in experiments with pressurized fractures, [30][31][32][33] resulting in pressure oscillations at the inlet. Stepwise fracture propagation has also been argued to be one of the sources of intermittent earthquakes. 34 Results from simulations, however, rarely show a stepwise propagation. It was observed to depend on the element size, 35,36 but disappeared upon mesh refinement, [37][38][39] and therefore, was concluded not to have a physical origin.
Simulations exhibiting stepwise fracture propagation were, however, published in Refs. 40, 41, but the numerical scheme employed may have promoted this phenomenon. Stepwise propagation of pressurized fractures that are independent of the element size has also been shown using a hybrid peridynamics/finite element approach, 42 resulting in pressure oscillations at the inlet similar to those observed experimentally.
Most of this research was aimed at pressurized fractures, while stepwise fracture propagation in shear has also been observed experimentally. Stepwise propagation was observed by Ref. 43 when using coarse meshes, but this disappeared upon mesh refinement, even though it was claimed to have a physical origin. 29 Therefore, we believe that so far there are no conclusive simulations that reveal the existence of a stepwise shear fracture propagation, and under which conditions and for which physical parameter values, they occur.
Our aim is to investigate stepwise shear fracture propagation observed in fluid-saturated porous materials. For this purpose, we will first derive the governing equations, neglecting as few terms as possible, while adhering to the standard displacement-pressure formulation. Nonassociated plasticity is included in this scheme to accurately represent reality, and a Cosserat continuum is used to resolve mesh sensitivity and convergence issues that can originate from the use of nonassociated plasticity. The result is a monolithic scheme that includes nonassociated plasticity, shear-based fracture propagation, and acceleration-driven fluid flows.
After deriving the governing equations and describing the method used to include all acceleration terms in Section 2, we will show the effects of including these acceleration terms in Section 4, informing the decisions made for which terms are included in the remaining analysis. Through a mesh refinement study and simulations with varying permeability, we will analyze several types of stepwise fracture propagation in Sections 5. Finally, a second case is simulated in Section 6, showing that the occurrence of stepwise propagation of shear fractures is case dependent.

GOVERNING EQUATIONS
The porous material is modeled as a Cosserat continuum, describing the kinematics of the (two-dimensional) solid by the displacements and and a Cosserat microrotation . 44,45 A Cosserat continuum introduces an internal length scale . It removes the ill-posedness and the ensuing the mesh dependence caused by the nonassociated plasticity, [46][47][48][49] and maintains the quadratic convergence of Newton-Raphson schemes used to solve nonlinear boundary value problems. 50

Including fluid inertia terms
The porous medium is described using the total momentum balance of the solid and fluid, the mass balance for the fluid, and the momentum balance for the fluid combined with Darcy's law. 18,19 These equations can be written in a form that incorporates a Cosserat continuum as: using the Biot modulus , the Biot coefficient , the intrinsic permeability , the fluid viscosity , the porosity , and the interstitial fluid pressure (with compression taken as positive ], using the matrix to only select the displacements, and to eliminate the Cosserat micro-rotation. The matrices used in Equations (1)-(3) are for a two-dimensional domain given by: with and the density of the solid and the fluid, respectively, and Θ = 2 2 ∕(1 + ) the rotational inertia. 48 The terms related to the convective momentum transport in Equations (1) and (3), ( ⋅ ) , require separate degrees of freedom to be used for the fluid flux in order to obtain the gradients of the fluid flux. However, the effect of these terms is usually very small, 18 and these terms will therefore be neglected. This allows the interior of the porous material to be described using only the solid displacement and the interstitial fluid pressure (and history variables).
To eliminate the fluid acceleration term in Equation (1), the definition of the Darcy fluid flux is substituted for the fluid acceleration term, resulting in: with volume averaged density = (1 − ) + . The effect of including separate inertia terms, compared to the more commonly taken approach of assuming the fluid inertia terṁto be negligible, 18,42,[51][52][53] will be investigated in Section 4.
The acceleration of the fluid at time + Δ is determined using a -scheme: Using this temporal discretization, the momentum balance of the fluid combined with Darcy's law, Equation (3), is rewritten in a partial time discretized form as: Using the values at the old time step, anḋ, as history variables results in an explicit expression for the fluid flux, only dependent on the displacement and pressure: This expression is substituted in the mass conservation for the fluid, Equation (2), resulting in the partial time-discretized mass balance at time + Δ : Similarly for Equation (4), discretizing the time derivative of the fluid flux using Equation (5), and substituting the expression of the fluid flux at time + Δ , Equation (8), results in: with the constant A defined as: Simplifications are often made by assuming the porous medium to be either drained or undrained. In the drained limit, the pressure throughout the domain is equal to a reference pressure, and taken to be constant. This assumes that the interstitial fluid instantly attains a uniform pressure, and that the effects of the inertia related to the fluid fluxes required for this are negligible. In contrast, the undrained limit assumes no fluid flow, representing an impermeable porous material. Then, the only source of pressure changes is the compression of the solid and the interstitial fluid.
When the intrinsic permeability is low, the constant and the fluid flux both approach zero, and as a result, Equation (10) contains only a single inertia term. This is consistent with undrained formulations, in which the volume averaged density is used. All terms related to fluid fluxes also cancel from the mass balance, Equation (9), and = 0 corresponding to the undrained formulation.
In contrast, when approaching the drained limit by using a large permeability, becomes one. For this limit, the solid inertia term only depends on the density of the solid, with the fluid using a separate inertia term solely based on the pressure gradient. The effect of using a large permeability in the mass balance results in the acceleration-driven fluid flow becoming more significant. This can be observed from Equation (8). When a high permeability is used such that ( )∕( Δ ) >> 1, the expression for the fluid flux becomes: This shows that the fluid flux is limited for higher permeabilities, and becomes independent of the intrinsic permeability and viscosity in this limit. This is in contrast with the drained assumption, where it is assumed that in the highpermeability limit, the fluid flux becomes large enough to limit any interstitial pressure changes. Therefore, differences between the drained assumption and the high-permeability cases are to be expected, with the cases with actively simulated fluid flow not approaching the drained limit. Instead, these high-permeability cases will approach the limit indicated by the above equation.

Plasticity model
The stresses in the solid material are assumed to be linearly related to the elastic strain: with the linear-elastic stiffness matrix for a Cosserat continuum 49 and the total strain +Δ = +Δ . A nonassociated Drucker-Prager plasticity model is used to model the plastic deformations, with the yield function and plastic potential defined as: with the solid pressure (tension positive), 2 the second invariant of the deviatoric stresses in a Cosserat continuum 45,54 and , , and related to the angle of internal friction , dilatancy angle , and cohesion through: More implementation details related to the integration of the plastic strain and the used tangential stiffness matrix are given in Ref. 50.

Discontinuity models
The discontinuity is modeled using interface elements, which are only present for the fracture. Ahead of this fracture, the material is considered fully intact and is governed by the equations for the poroelastic medium. It is assumed the fracture propagates in mode II, and is solely governed by the shear stress. Once the stresses in the continuum ahead of the fracture exceed the fracture criterion | | > − , 55 the discontinuity is propagated and additional interface elements are inserted. This propagation criterion uses the peak coefficient to model a cohesionless material, with the required shear stress for fracture propagation depending on the normal stress.
At the discontinuity, a discontinuous pressure model 39,56,57 has been used, which allows for a discontinuous pressure distribution across the fracture, with an independent pressure degree of freedom inside the fracture. For closed fractures, without fluid flow inside the fracture in the tangential direction, the mass balance inside the fracture is given by: and the fluid flux imposed on the porous material is given by: using the pressure inside the discontinuity and the interstitial fluid pressures at the top and bottom walls of the fracture, + and − , respectively. The normal vectors at the top and bottom of the fracture, + and − , are taken as positive when pointing outward. The interface permeability governs the amount of fluid that can flow across the fracture. It is an indication of the resistance to the fluid flow originating from plastic deformations and damage close to the discontinuity caused by the fracturing process. This interface permeability therefore allows to include the effects of damage and plasticity on the fluid flow, even though the region of these effects is smaller than the element size. To keep the formulation as simple as possible, the interface permeability will be assumed as constant.
In a local coordinate system, the tractions at the discontinuity are given by: To prevent negative fracture opening heights, we set the following relations: with high values for the dummy stiffnesses and , and and the jumps in the normal displacement and the microrotation, respectively. The tangential component of the interface traction is determined through an exponential traction separation law: with and the peak and residual coefficients of friction. is the weakening distance, given by = 2 ∕( − ) where  is the fracture energy. The used exponential relation was preferred over the more commonly used linear degradation model 10,55,58 in order to retain a smooth function, which improves the convergence of the nonlinear solver. The maximum jump in the tangential displacement dx is treated as a history variable, and therefore, it is equal to the current tangential displacement jump , or constant when the current displacement is lower than the previously obtain tangential displacement.

Discretization
The time discretization of Equations (9) and (10) has been carried out using a Newmark scheme for the displacements and a -scheme for the interstitial fluid pressure: This introduces history variables for the solid velocity and acceleration, and the fluid pressure change at time . These history variables, and the old fluid flux and fluid flux changė, are updated at the end of each time step. In contrast, the history variable for the plastic strain, , and the interface displacement, dx , are updated at the end of each time step and after each fracture propagation step. This ensures these path-dependent variables accurately capture the loading and unloading that can occur during fracture propagation. The fracture propagation is considered to be part of the time step. Therefore, once a converged solution is achieved, the stresses ahead of the fracture tip are checked and if all these stresses exceed = − , the fracture is propagated a single element (checked at both ends of the fracture). Afterwards, more iterations are carried out to again establish an equilibrium solution for the current time step and the new fracture length. Next, the propagation criterion is checked again. This is repeated until no further propagation occurs and a converged solution is reached, resulting in an implicit scheme for the governing equations and the fracture length.
T-splines have been used for the spatial discretization. 59,60 This allows a higher-order interelement continuity, resulting in continuous stresses and fluid fluxes between the elements. 61,62 Furthermore, by using T-splines, the displacement gradient, and thus the strains and stresses, can be directly evaluated at the boundaries of elements. This facilitates the fracture propagation criterion as it is not necessary to extrapolate stresses from the interior of the elements. The T-splines are used in a similar manner as traditional Lagrangian interpolants by using Bézier extraction, decomposing the splinebased functions in a linear combination of Bernstein polynomials. 63 This allows the displacement, the interstitial fluid pressure, and the pressure in the discontinuity to be expressed as a sum over all elements: with: To fulfill the inf-sup condition, 64 and thus to prevent spurious oscillations, the interpolants for the solid displacement need to be an order higher than the interpolants for the fluid pressure. Therefore, quartic T-splines were used for the solid displacement, and , whereas cubic T-splines were used for the interstitial and discontinuity fluid pressures, and (only defined at the discontinuity). The inf-sup condition does not impose any requirements on the order of the interpolants for the Cosserat rotation , and it was chosen to use cubic T-splines to reduce the resulting number of degrees of freedom. For the elements and geometry of the quartic and cubic meshes to match, the interelement continuity of the meshes needs to be the same. This is achieved by repeating the mesh lines of the quartic T-spline mesh, resulting in both meshes having a 2 interelement continuity, thus allowing for continuous and smooth stresses and fluid fluxes.
The discontinuity is represented using spline-based interface elements. [65][66][67][68] By using T-splines, these interface elements can be inserted only for the fractured elements, in contrast to B-splines that require the interface elements to be inserted along a line passing through the complete domain. This removes the need for a formulation for nonfractured interface elements. 69 The displacements and Cosserat microrotation are discontinuous across the interface elements, 70 with the tractions at the interface determined by Equation (18). The interstitial fluid pressure is also discontinuous across the discontinuity, with an additional degree of freedom added at the discontinuity to represent the pressure inside the discontinuity. The fluid flux between the porous material and the interior of the discontinuity is governed by Equations (15)- (17).
Using the discretization from Equations (21)-(26), the momentum balance from Equation (10) is discretized as: with the external force vector defined in a standard manner as: where is the traction imposed on the boundary of the domain, and is the damping matrix used for the absorbing boundary conditions: 50,71 using the pressure and shear stress wave speeds and . The internal force vector is given by: using = . Finally, the forces related to the discontinuity are obtained from: with the rotation matrix at the discontinuity. Similarly, mass conservation, Equation (9) is discretized as: with the external fluxes defined as: and the internal flux vector as: The fluxes due to the fluid flow normal to the fracture are given by: and the mass balance for the interior of the fracture, Equation (15), is discretized as: Equations (28), (32), and (36) are solved in a monolithic manner using a Newton-Raphson scheme. The system matrices used for this scheme are given in Appendix A. A quadratic convergence rate was obtained using these system matrices.

DESCRIPTION OF THE CASE STUDIED
The effect of including pressure-and acceleration-driven fluid flows is now demonstrated by simulating a typical case that features a propagating shear fracture. 10 The spatial discretization uses cubic T-splines for the interstitial pressure, the discontinuity pressure, and the Cosserat microrotation, while quartic T-splines were used for the displacements. Three different meshes were used. The first mesh consists of a layer of 320 × 10 Bézier-extracted elements near the interface, surrounded by layers of 160 × 5 and 80 × 15 elements (hereafter referred to as 320 × 50). The second mesh adds an extra refinement layer of 640 × 10 elements at the interface, for a total of 60 elements in the vertical direction (referred to as 640 × 60). Finally, the third mesh halves all the elements, resulting in a mesh of 1280 × 120 elements.
The conditions at = 0 are a zero interstitial fluid pressure and a zero pressure at the discontinuity, with the solid being in equilibrium and without any plastic deformations at this time. This is imposed by simulating the first 0 = 1 s of the simulation with time steps of Δ = 0.01 s with no plastic strains or fracture propagation being allowed, and a constant and zero pressure imposed in the domain. During this initial period, the discontinuity is not allowed to propagate and no plastic strains are allowed to occur (both the propagation criteria and the yield criteria are ignored). Furthermore, the interstitial pressure is constrained to = 0 MPa in the complete domain. After this initial period, starting from = 0 , the simulation is carried out using a time step Δ = 0.04 ms without the constraints that have been applied initially. The time step is such that even with the finest mesh, the fracture will only propagate over two elements (due to the symmetry) every 5-10 time steps. The temporal discretization for the solid displacement used = 0.4 and = 0.75 as parameters, and for the temporal integration of the fluid = 1.0.

EFFECT OF THE FLUID INERTIA TERMS
The effect of terms that account for acceleration-driven fluid flow in the mass balance and separate inertia terms in the momentum balance is analyzed by adding these terms one-by-one. The first simulation uses only the pressure gradientdriven flow, the term in Equation (6), and uses a single inertia term, thus neglecting the acceleration of the fluid relative to the porous material,̇= 0 in Equation (4). Next, the solid acceleration-driven flow is added,̈in Equation (6), and after that the acceleration of the fluid relative to the solid is taken into account for the mass balance,̇in Equation (6). Finally, two separate inertia terms are used to include the relative acceleration in the momentum balance by includinġfrom Equation (4). All simulations have been carried out using the 640 × 60 mesh, and an intrinsic permeability = 10 −8 m 2 . The interstitial fluid pressure that results from these simulations is shown in Figure 2. Without acceleration-driven flow, only slight changes in pressure occur due to the compression of the porous material, as shown in Figure 2(A). These high-and low-pressure regions are diffused over a large area because of the high permeability.  Figure 2(B). For the sides of the fracture that are under compression (top right and bottom left), this causes lowpressure regions at the location of the stress waves that originate from the fracture tips, encapsulated by two high-pressure regions that result from the accelerations forcing the fluid away from the stress waves. However, pressure oscillations can be seen throughout the domain due to this coupling between the fluid fluxes and solid acceleration.
Adding the fluid acceleration term in Darcy's law removes these interstitial pressure oscillations, as shown in Figure 2(C). This is a result of the fluid acceleration term adding additional damping to the formulation, thereby decreasing the magnitude of the stress and pressure waves away from the discontinuity. Adding the fluid acceleration term also increases the pressure at the fracture tips, resulting in higher pressures on the compression sides, and lower F I G U R E 3 Shear stress at time − 0 = 52 ms when different acceleration terms are included in the formulation pressures on the extension sides compared to the cases without this term. Finally, using separate inertia terms in the mass balance does not result in any significant changes in the interstitial fluid pressure ( Figure 2D).
The shear stress is shown in Figure 3. Similar to the fluid pressure, slight stress oscillations occur when only the solid acceleration term is included, and disappear if the fluid acceleration term is added. Using separate inertia terms results in slight changes in shear stress, mainly near the center of the domain. However, the effect of the acceleration and separate inertia terms on the shear stress (and other stress components, not shown here) is small compared to the effects on the fluid pressure.
From these results, it can be concluded that the acceleration-driven flow has a significant effect on the fluid pressure, and needs to be included in the formulation. The inclusion of separate inertia terms, however, results in lengthy terms F I G U R E 4 Effect of mesh refinement on the fracture propagation in the internal force vector, while having near to no effect on the resulting pressures and stresses. In the remainder of this paper, only the acceleration-driven fluid flow terms will therefore be included, while neglecting the effects of using separate inertia terms.

Mesh refinement
The three different meshes were used to simulate the undrained, the drained, = 10 −8 m 2 , and = 10 −10 m 2 cases to investigate the effect of the interface element size on the fracture propagation. For the drained (not shown), the undrained ( Figure 4B), and the = 10 −10 m 2 cases, the interface element size does not significantly alter the fracture propagation.
With the coarse 320 × 60 mesh, small steps in the fracture length are seen, with the size of these steps corresponding to F I G U R E 5 Effect of mesh refinement on the fracture propagation using = 10 −8 m 2 F I G U R E 6 Fracture propagation length for different values of the intrinsic permeability two interface element lengths. These steps originate from the element-wise fracture propagation resulting in a minimum propagation length corresponding to the element size. If the shear stress is not sufficient to fracture a complete element, as evaluated in all the integration points of the interface element, no fracture propagation occurs. Once the stress rises and exceeds the fracture toughness, the complete element fractures at once, resulting in a jump in the fracture length. Due to the symmetry of the simulated cases, this results in either two or no elements fracturing, and thus creating jumps in the fracture length equal to two element lengths. These jumps decrease upon mesh refinement, but do still occur even for the finest mesh as shown in Figure 5. Therefore, this stepwise fracture propagation is solely caused by the spatial and temporal discretization of the problem, and does not originate from a physical phenomenon. It is important to note that if the time step size is increased such that the fracture propagates at least once during each time step, these steps in the fracture propagation will disappear since the propagation is assumed to occur as part of the time step. This is in contrast to fracture propagation schemes that assumes the fracture to propagate in between time steps, 40,41 where these steps will remain independent of the used temporal discretization.
The results for the = 10 −8 m 2 case are shown in Figure 4(A). In contrast to the other cases, there is a clear difference between the three meshes. Next to the steps originating from the element size, the coarsest mesh also shows a temporary arrest in the fracture propagation. Upon mesh refinement, these arrests become more common. The finest mesh also shows that the duration of the arrests increases. The difference between this manner of stepwise fracture propagation and the element-sized steps observed for all cases is clearly seen in Figure 5. The small steps occur for all meshes, but the pronounced stick-slip like behavior results in large plateaus when plotting the fracture length as a function of time. While this behavior occurs for all meshes, only the fine meshes resolve the fracture propagation in sufficient detail to distinguish the smaller arrests close to the start of the simulation from the element-wise fracture propagation. Furthermore, while the propagation velocity does not differ between the coarsest and finest mesh, the additional stick-slip behavior results in a significantly slower fracture propagation. F I G U R E 7 Plastic strain (( 2 + 2 ) 1∕2 ) at time − 0 = 120 ms around the right half of the discontinuity

Effect of the permeability
Results for the complete range of permeabilities that have been analyzed are presented in Figure 6 for the 1280 × 120 mesh. As also shown in the previous section, stick-slip like behavior is observed for the simulations with a high permeability, whereas this does not occur for lower permeabilities, and neither for the drained and undrained cases. Due to the lack of stick-slip behavior, the = 10 −10 m 2 case approaches the fracture propagation length of the undrained case, while for high values of the permeability the fracture propagates slower than for the drained case. Relative little difference is seen between = 10 −8 m 2 and = 10 −7 m 2 , suggesting that for these values, the behavior of a highly permeable porous material is approached. This indicates that while the undrained limit is a good approximation for cases with a low permeability, the drained limit neglects important physical phenomena and, accordingly, does not correctly approximate the behavior for high values of the permeability.
The plastic strain near the right fracture at time − 0 = 120 ms is shown in Figure 7. The undrained and drained cases show a similar behavior as was observed in Ref. 11, with the plasticity being limited to the extensional side of the fracture, and plastic deformations occurring over a larger area for the drained case compared to the undrained case. While the = 10 −10 m 2 case approximates the results from the undrained case, the plastic strains obtained for the higher permeability cases differ significantly from the drained case. Large plastic deformation peaks are observed at the locations of the pauses, while more plastic deformation is seen between these pauses. In contrast to the other cases, the higher permeability cases also display a slight amount of plastic deformation on the compression sides of the discontinuity. This further confirms that the high permeability cases do not approach the behavior of the drained limit case.

Explanation for stick-slip like fracture propagation
To explain the cause of the stick-slip behavior, we consider the simulation for = 10 −8 m 2 at four different times, see Figure 8. The interstitial fluid pressure is given in Figure 9 and the shear stress, causing propagation of the fracture, is shown in Figure 10. When the fracture propagation restarts, there is a large low-pressure region on the compression side next to a small high-pressure region, while the reverse takes place at the extension side ( Figure 10A). This small high-pressure region increases in size and attains a maximum during crack propagation, while the low-pressure region detaches and diffuses away from the discontinuity ( Figure 10B). Simultaneously, a new low-pressure region emerges behind the growing high-pressure region. Due to the inclusion of the acceleration-driven fluid flow, fluid is transported away from this low-pressure region, further decreasing its pressure, toward the high-pressure region. This allows a build-up of the pressure near the crack tip. Once this region starts to overtake the crack tip ( Figure 10C), crack propagation arrests. However, due to the acceleration terms in Darcy's law, the high-pressure region continues to overtake the crack tip ( Figure 10D). Only when the high-pressure region above the fracture, and the low-pressure region below it have overtaken the crack tip, fracture propagation is resumed and a new high-pressure region starts to build up behind the crack tip. The effect of alternatingly building up high-and low-pressure regions influences the shear stresses near the interface, as shown in Figure 10. While the shear stress is concentrated around the crack tip at the beginning of fracture propagation ( Figure 10A), the high-and low-pressure regions diffuse this stress over a growing region (Figures 10B and 10C). Once the fracture has arrested, this shear stress is distributed over a significantly larger region compared to the beginning of the fracture propagation. The effect of the pressure overtaking the fracture is also clearly seen in Figure 10(D). Once the F I G U R E 9 Interstitial fluid pressure around the right crack tip for = 10 −8 m 2 shear stress concentrations emerge again, as a result of the pressure regions being away from the discontinuity, crack propagation restarts.
Considering the results for a lower permeability ( Figure 11) or when the fluid acceleration terms are not included (Figure 12), we observe that these alternating pressure regions do not occur. For the low-permeability case, the acceleration-driven fluid flow terms are small compared to the pressure changes induced by the compression of the porous material. The shear stresses are also concentrated near the fracture tip for both these cases, in contrast to the case with a high permeability, where these stresses diffused over a larger region. This shows that the stick-slip like behavior is purely caused by the acceleration-driven fluid flow, and that excluding acceleration-driven fluid flow terms significantly alters the crack propagation.

SUPERCRITICAL SHEAR FRACTURE PROPAGATION
A second case is now considered with a combination of external forces to cause the propagation velocity of the fracture to exceed the shear wave speed. Now = 19 MPa, = 10 MPa, = = 2.4 MPa, while all everything else corresponds to the previous case study, Section 3. The crack velocity, obtained using the 1280 × 120 mesh, is shown in Figure 13(A). Different from the previous case, no stick-slip like behavior occurs. Furthermore, the simulations for different values of the permeability all approach the undrained case, while the crack in the drained case propagates slightly slower. This can be explained by considering the interstitial fluid pressure in Figure 14. The simulation for = 10 −8 m 2 still results in a high-pressure region next to a low-pressure region on the compressive side of the fracture. However, due to the high propagation velocity, this highpressure region cannot overtake the fracture tip, and therefore, the interstitial fluid pressure does not initiate the stick-slip behavior observed before. Comparing the results for = 10 −8 m 2 , = 10 −10 m 2 and the undrained case also shows why the propagation velocity is similar in those cases. No interstitial fluid pressure is present in front of the fracture, whereas the combination of acceleration-driven fluid flow (for the high-permeability case) and compression of the porous material results in similar pressures just behind the fracture tip. While no stick-slip behavior originating from the interstitial fluid pressure is observed, small arrests in the crack propagation still occur for all cases, see Figure 13(B). They are independent of the value of the permeability, and only occur up to − 0 ≈ 40 ms (corresponding to a fracture from = 125−375 m). Thereafter, arrests are no longer observed, only elementwise fracture propagation. The cause of these small arrests lies in the shear stress ( Figure 15). Before the onset of fracture propagation, the domain is in equilibrium. However, once propagation occurs, it happens faster than the shear stress can react, resulting in the fracture propagating through the initially present shear stress and causing oscillations. These oscillations are responsible for the short arrests. Eventually, the fracture exits the region that was influenced by the initial position of the fracture and these stress oscillations disappear.

CONCLUDING REMARKS
A pressure-displacement formulation has been presented for fluid-saturated porous media that includes accelerationdriven fluid flow and separate inertia terms. It has been shown that these acceleration-driven flow terms can have a significant impact. Including the solid acceleration results in the interstitial fluid pressure reacting to the pressure and the shear stress waves, but also causes oscillations in the interstitial fluid pressure. Adding the fluid acceleration adds extra damping, removing these oscillations and limiting changes in fluid pressure away from the source of the stress waves. The fluid acceleration term is also responsible for causing small high-and low-pressure regions near the fracture tip. Finally, simulations using separate inertia terms for the fluid and solid show near to no effect of using these separate inertia terms.
Simulations of a shear fracture in a porous material with varying permeability show several types of stepwise fracture propagation. The first type of stepwise fracture propagation is related to the elementwise fracture propagation, and was shown to gradually vanish upon mesh refinement. These steps originate from the small time step used during these simulations, causing either zero or two elements to fracture. This necessitates small time steps. This type of stepwise propagation shows a consistent pattern of very short arrests followed by a jump in fracture length equal to a single element (or two-element lengths in case of symmetric problems). This type of stepwise propagation is a numerical artifact, and does not have a physical origin.
The second type of stepwise propagation observed was a stick-slip like behavior that originates from the accelerationdriven fluid flow. This results in alternating high-and low-pressure regions near the fracture tip, arresting the fracture propagation once one of these regions starts to overtake the fracture tip. This stick-slip like behavior could be observed only when a sufficiently fine mesh was employed, and only occurred for simulations using a high value of the permeability, so that the acceleration terms could dictate the changes in interstitial fluid pressure.
Finally, simulations have been carried out using a combination of external forces that resulted in a crack propagation velocity exceeding the shear wave speed. For these simulations, the interstitial fluid pressure build-up is unable to overtake the fracture, and therefore, the second type of stepwise propagation does not occur. However, a third type of stepwise propagation was observed, originating from the initial conditions not adapting fast enough to the propagating fracture. The fast fracture propagation resulted in stress oscillations when propagating through the initial stress, with these stress oscillations causing small arrests in the fracture propagation. These arrests disappear once the fracture exits the region influenced by the initial stress field.

A C K N O W L E D G M E N T
Financial support through H2020 European Research Council Advanced Grant 664734 "PoroFrac" is gratefully acknowledged.

C O N F L I C T O F I N T E R E S T
The authors declare that there is no conflict of interest regarding the publication of this paper.

D ATA AVA I L A B I L I T Y S TAT E M E N T
Data sharing is not applicable to this article as no datasets were generated or analyzed during the current study.