A phase-field model for quasi-dynamic nucleation

Despite its critical role in the study of earthquake processes, numerical simulation of the entire stages of fault rupture remains a formidable task. The main challenges in simulating a fault rupture process include complex evolution of fault geometry, frictional contact, and off-fault damage over a wide range of spatial and temporal scales. Here, we develop a phase-field model for quasi-dynamic fault nucleation, growth, and propagation, which features two standout advantages: (i) it does not require any sophisticated algorithms to represent fault geometry and its evolution; and (ii) it allows for modeling fault nucleation, propagation, and off-fault damage processes with a single formulation.

 

Built on a recently developed phase-field framework for shear fractures with frictional contact, the proposed formulation incorporates rate- and state-dependent friction, radiation damping, and their impacts on fault mechanics and off-fault damage. We show that the numerical results of the phase-field model are consistent with those obtained from well-verified approaches that model the fault as a surface of discontinuity, without suffering from the mesh convergence issue in the existing continuous approaches to fault rupture (e.g. the stress glut method). Further, through numerical examples of fault propagation in various settings, we demonstrate that the phase-field approach may open new opportunities for investigating complex earthquake processes that have remained overly challenging for the existing numerical methods. Keywords: Phase-field modeling, Fault rupture, Rate-and-state friction, Radiation damping, Off-fault damage, Earthquake 1. Introduction Earthquakes are classically thought of as manifestations of frictional instabilities associated with the nucleation and subsequent rapid growth of slip on geologic fault surfaces. The earthquake processes involve a wide range of complexities arising from the fault geometry, frictional contact, and off-fault damage (i.e. microcracks and temporal variations in the bulk wave speeds).

 

. Since these complexities are not tractable using analytical methods, except for a limited set of idealized problems, numerical methods play an indispensable role in investigating source physics and its implications on seismological, geological, and geodetic observations. For decades, the primary means for numerical modeling of fault rupture has been discontinuous approaches that treat the fault motion as a displacement jump (relative displacement) across a zero-thickness interface. Notable examples include the boundary integral method [11–14], and the traction at split node method which has been applied to finite difference [15–17] and finite element computations [18, 19]. The popularity of these discontinuous approaches may be credited to the widespread applications of fracture mechanics in earthquake studies, in which the fault is considered a shear fracture with zero thickness. This discontinuous viewpoint has further been justified by the fact that the fault slip is primarily localized in a region with a negligible thickness relative to both the fault length and the wavelengths of interest in the seismic wavefield [20–23]. For a discontinuous description of fault rupture, however, one should explicitly handle the geometric complexities involved in the seismogenesis—a particularly unsettling challenge when the fault path and/or the damage in the surrounding bulk would evolve. For example, because the boundary integral method cannot model the material nonlinearity in surrounding rock masses, it is intrinsically unable to capture the growth of off-fault damage. The traction at split node approach can simulate damage in rock masses by incorporating bulk inelasticity [9, 19, 24], but it requires the fault interface to be aligned with the element boundaries. Therefore, the traction at split node method needs a sophisticated remeshing algorithm to model fault propagation in arbitrary directions. This limitation may be overcome by the use of embedded discontinuity methods such as the extended/generalized finite element method (XFEM/GFEM), which can accommodate the fault interface in the interior of discretized elements [25, 26].

 

Yet such embedded discontinuity methods may demand significant efforts for implementation as they require enrichment of shape functions and sophisticated algorithms for numerical integration. An alternative approach to the numerical modeling of fault rupture is to treat the fault surface as a continuous entity with finite thickness. Commonly used continuous approaches include the stress glut method [16, 27, 28] and thick fault zone methods [28–31]. In contrast to the discontinuous approaches, these continuous approaches model the fault as an inelastic zone of finite thickness across a single or multiple layers of elements. The fault motion is then represented by an inelastic shear strain, which can be simply calculated according to a specific constitutive law for friction. As a result, neither a remeshing algorithm nor explicit calculation of displacement jump is required for these continuous approaches, allowing them to be easily implemented in standard numerical methods. However, the results of these continuous approaches are quite different from those produced by the discontinuous approaches. Dalguer and Day [28] have shown that the results of the thick fault zone methods are not even qualitatively consistent with reference solutions obtained by standard discontinuous methods. The results of the stress glut method show qualitative agreement with the reference solution, but they do not converge to the reference solution upon mesh refinement [28, 31]. This non-convergence problem is attributed to the fact that these formulations lack a length scale that is necessary to retain mathematical 2 well-posedness during softening behavior. Such a convergence issue is absent in the very recent model proposed by Gabriel et al. [32], which describes dynamic rupture processes in a low velocity fault zone using a regularized damage formulation.

 

Yet, because the model is derived from damage rheology, there is no a priori constraint on the results to converge to the laboratory-derived rate-and-state friction limits [33]. Also, while this approach is particularly promising, it entails several parameters that are not standard in earthquake studies and so would not be easy to calibrate. Therefore, it is highly desirable to develop a continuous approach that consistently converges to the discontinuous fault description while relying on standard mechanics theory for shear fracture with frictional contact [34]. Over the last decade, in the computational mechanics community, the phase-field method has emerged as a novel approach to continuous modeling of fracture. The method diffusely approximates the sharp geometry of a crack surface with a spatially distributed variable—the phase field—and describes the evolution of the crack surface by a partial differential equation formulated from fracture mechanics theory. In this way, the phase-field method handles complex crack geometries that are neither aligned with mesh boundaries nor do they require enrichment functions, lending itself to a simple implementation in standard numerical frameworks such as the finite element method. Furthermore, the phase-field method is free of the problems of the stress glut and thick fault methods. Specifically, as the phase-field formulation is rooted in fracture mechanics theory, the phase-field solutions are consistent with those of discontinuous approaches based on the same theory. Also, because the phase-field approximation introduces a length scale for regularizing a sharp discontinuity, it is mathematically well-posed and hence its numerical solutions converge upon mesh refinement. However, as the phase-field method was originally developed for tensile (mode I) fracture, the vast majority of phase-field models in geomechanics have focused on tensile fracture, e.g. [35–40]. Recently, the capabilities of phase-field modeling have been extended to frictional shear fractures in geologic materials.

 

Fei and Choo [41] have developed the first phase-field formulation for discontinuities with frictional contact. Later, Fei and Choo [42] have extended the formulation to propagating shear fractures, in a way that is demonstrably consistent with the fracture mechanics theory proposed by Palmer and Rice [34]. Also, Bryant and Sun [43] have incorporated rate- and state-dependent friction, which is considered essential to earthquake modeling, into the phase-field formulation for frictional discontinuities. Yet none of the existing phase-field formulations is sufficiently credible for modeling fault rupture processes in earthquakes. The current phase-field models for frictional shear fracture are all restricted to the quasi-static condition, neglecting dynamic effects or their quasi-dynamic approximation. The (quasi- )dynamic nature, however, is essential to the fault rupture modeling as it is responsible for several seismic features including fast rupture propagation and fault slip rates in the m/s range which correlates with strong ground motion. Also, a quasi-static model becomes unstable as the rupture accelerates and static equilibrium is no longer possible. It is thus important to explicitly retain the inertia term (a fully dynamic approach) or its approximation through radiation damping (a quasi-dynamic approach), to ensure the numerical stability of the calculation during the coseismic period, as well as to obtain results that are seismologically relevant. In this work, we develop the first phase-field approach to fault rupture and propagation. First, we incorporate rate- and state-dependent friction into the most recently proposed phase-field framework for geologic discontinuities [44], whereby the displacement-jump-based kinematics of faults are transformed into 3 a strain-based version and then inserted into the phase-field formulation for frictional interfaces. Second, we adopt, in this initial study, a quasi-dynamic formulation that approximates the inertial effects by a radiation damping term in the calculation of the shear stress on the fault plane [45]. Compared with the fully-dynamic modeling, the quasi-dynamic formulation significantly reduces the computational cost, while retaining the numerical stability and qualitative patterns of the simulation results. Therefore, the quasi-dynamic formulation has been popular in numerical investigations of fault rupture processes, e.g. [46–54].

 

Drawing on the same idea, we incorporate quasi-dynamic effects by augmenting a radiation damping term to the phasefield formulation for frictional shear fracture under quasi-static conditions [42]. The resulting phase-field formulation is consistent with the earthquake energy budget and enables energy partitioning between strain energy, frictional work, and radiated energy approximated through the radiation damping. We demonstrate the potential capabilities of the proposed formulation as well as its performance through a number of numerical examples with varying complexity. At this point, we clarify two important limitations of the scope of this initial study. First, as explained above, the formulation is limited to quasi-dynamic conditions. Second, we shall limit the fault kinematics to a two-dimensional antiplane condition. Extension of the work to the fully dynamic limit and/or more complex kinematics will be presented in future work.