Nonlinear Radiation Damping: A New Method for Dissipating Energy in Dynamic Earthquake Rupture Simulations

Dynamic, also called “spontaneous,” earthquake rupture computational simulations are an appealing approach to understanding how earthquakes operate. This type of simulation is able to incorporate many features affecting earthquake rupture behavior. These include interactions among fault geometry, rock properties, stress conditions, and fault friction, and the resulting wave propagation (e.g., Harris et al., 2018; Ramos et al., 2022). The simulations are used to better understand earthquakes which have already occurred, and also to foresee potential earthquake behavior in the future, for generic or specific tectonic settings (e.g., Lozos and Harris, 2020; Harris et al., 2021; Taufiqurrahman et al., 2022). Dynamic rupture simulations often assume that rocks behave as a linear elastic material. This assumption means that it is not necessary to know the initial stress tensor throughout the model volume to solve the dynamic elasticity equations; the solution depends only on the change in stress from an unspecified initial state. It also means that one is free to specify any desired initial tractions on the fault, and that the model need not include gravitational body forces, because it can be assumed that the unspecified initial stress state is in static equilibrium and gives the desired tractions when resolved onto the fault surface. All of this makes it straightforward for a linear elastic model to incorporate features such as a nonplanar 3D fault surface and a 3D seismic velocity model with heterogeneous rock density.

The assumption of elastic material properties often produces peak slip rates and ground motions higher than have been inferred from seismological observations. This has led some studies to implement viscoplastic off‐fault yielding (e.g., Ma and Andrews, 2010; Dunham et al., 2011; Gabriel et al., 2013; Roten et al., 2017) in their simulations. Viscoplasticity can be an effective way to dissipate energy and thereby reduce slip rates and ground motions, but with very few exceptions (e.g., Andrews et al., 2007; Duan and Day, 2010), the appropriate viscoplastic parameters for modeling rocks in the Earth are unknown. It is also necessary to include the full stress tensor throughout the model volume, another unknown. A 3D viscoplastic model must include gravitational body forces to achieve static equilibrium in the initial state, and the initial fault tractions can no longer be specified arbitrarily but instead are determined by resolving the initial stress tensor onto the fault surface. (Specifically, the initial stress tensor must satisfy the static equilibrium equation $σij,j=ρϕ,j$⁠, in which $ρ$ is density and $ϕ$ is the gravitational potential, with a zero‐traction boundary condition at the Earth’s surface. At the fault it must satisfy $Ti=σijnj$⁠, in which $Ti$ is the initial traction on the fault and $nj$ is the normal to the fault surface. The initial fault traction must be consistent with the friction law in that the initial shear traction is not so high that the fault ruptures immediately, nor so low that the fault can never rupture. Likewise, the initial stress must be consistent with the assumed viscoplastic parameters in that the stress is not so large that the viscoplastic material yields immediately.) It can be difficult to obtain an initial stress tensor that satisfies all the requirements, especially if the model includes a nonplanar 3D fault surface or a 3D velocity and density model. In addition, a viscoplastic model can be expensive to compute; for the computer code used in this study, a viscoplastic model takes three times as long to execute as a similar elastic model.

Our goal is to present a new method for dissipating energy in a dynamic rupture simulation, thereby reducing peak slip rates and ground motions. The method works by adding a new term to the friction law in an elastic simulation, which represents nonlinear radiation damping. Because the method only modifies the friction law, and retains linear elastic rock properties, it is simple and efficient to implement. We demonstrate that our method can emulate many of the properties of a viscoplastic model, while keeping the relative simplicity of a linear elastic model.

For an earthquake fault, the emission of seismic waves is accompanied by inertial body forces that tend to oppose rapid slip. This effect is called “radiation damping.” Our method is inspired by the way that radiation damping is implemented in quasi‐dynamic earthquake simulations (Rice, 1993; Tullis et al., 2012; Richards‐Dinger and Dieterich, 2012). Because quasi‐dynamic simulations lack seismic waves, the effects of radiation damping are often emulated by adding a term to the friction law. The additional frictional traction is $GV(t)/2VS$⁠, in which G is the shear modulus, V(t) is the slip rate at time t, and $VS$ is the shear‐wave velocity. This is a linear function of V. The additional traction represents a back action on the fault due to inertia, which is otherwise not included in a quasi‐dynamic simulation.

Dynamic rupture simulations ordinarily do not contain explicit radiation damping terms, because such simulations already include the effects of inertia and wave propagation. But inelastic off‐fault yielding creates additional damping beyond what occurs in a purely elastic simulation. The physical process by which this occurs is that fault slip deforms the viscoplastic material near the fault, which causes the material to yield, and the yielding, in turn, reduces the shear stress acting on the fault and thereby reduces further slip.

Here, $Cr$ is a coefficient (with dimensions of stress) that controls the overall strength of the radiation damping, V(t) is the fault slip rate at time t, $V0$ is the slip rate at which radiation damping starts to become significant, and n > 1. The idea is that $τdamp$ is small when $V<V0$⁠, and $τdamp$ increases approximately linearly in V, when $V>V0$⁠. The exponent n controls the sharpness of the transition between these two regimes. Conceptually, $τdamp$ represents a reduction in the shear stress acting on the fault due to inelastic yielding near the fault, which is otherwise not accounted for in an elastic simulation.

This is the only change that needs to be made when applying radiation damping, showing that our method is simple to implement and computationally efficient. Although we use linear slip‐weakening friction in this study, the same method could be applied to rate and state friction or any other desired friction law.

To test and demonstrate our method, we use modified versions of Southern California Earthquake Center‐U.S. Geological Survey (SCEC‐USGS) dynamic rupture benchmarks TPV26v2 and TPV27v2 (Harris et al., 2018). These are a matched pair of benchmarks, designed to illustrate the effects of off‐fault inelastic yielding. TPV26v2 uses linear elastic material properties, whereas TPV27v2 is identical except that it uses Drucker–Prager viscoplastic material properties.

Our simulations assume a vertical planar strike‐slip fault in a half‐space, as illustrated in Figure 1. The supplemental materials contain complete details of our models, including the modifications that we made to the SCEC‐USGS benchmarks, the model parameters, and the methods used to calculate radiated and dissipated energy. We present and compare the results of four simulations, as shown in Table 1. Simulation E1 has elastic material properties and no radiation damping. Simulation V1 has viscoplastic material properties within 2000 m of the fault, and elastic material properties elsewhere. Simulations R1 and R2 have elastic material properties and differing values of the radiation damping parameters $Cr$⁠, $V0$⁠, and n. The radiation damping parameters shown in Table 1 apply at a depth of 10 km, the hypocentral depth. Because our model includes depth‐dependent stresses (which is necessary to perform runs with viscoplastic material properties), the radiation damping parameters used in the simulations also have depth dependence. In other models, with depth‐independent stresses, it would be appropriate for the radiation damping parameters to be depth‐independent as well.

We demonstrate that the simulations with nonlinear radiation damping, R1 and R2, can emulate some of the properties of the viscoplastic simulation V1. Parameters are selected so that R2 has stronger nonlinear radiation damping than R1. We tried using linear radiation damping (by setting n = 1 in equation 1), but it did not produce the same effect. Of course, one cannot expect radiation damping to produce results that agree with viscoplastic results in every respect.

For each simulation, Figure 2 shows the slip rate history at four locations on the fault surface, which are at along‐strike distances of 5, 10, 15, and 20 km from the hypocenter. The four locations are illustrated in Figure 1. For the elastic simulation E1, peak slip rate is approximately proportional to distance from the hypocenter. Because a linear model like E1 has no “preferred” slip rate, there may be a tendency for slip rates to grow large in proportion to the rupture size (Andrews, 1976). Introducing plasticity, as in simulation V1, not only lowers the peak slip rates but also causes them to increase more slowly, consistent with the peak slip rate leveling off at some “preferred” value rather than growing without limit. Simulation R1 has peak slip rates very similar to V1, demonstrating that nonlinear radiation damping can affect fault slip in much the same way as viscoplastic yielding. Simulation R2 has stronger radiation damping, and correspondingly lower peak slip rates. As the slip rate peaks get lower, they also get wider; this occurs because, at each of the four locations shown, all four simulations have similar amounts of final slip. (The final slip is similar because all the simulations have the same elastic modulii, initial stress, and friction coefficients, and therefore must slip by a similar amount in order for the shear traction to fall below the dynamic friction.)

Table 1 shows the radiated seismic energy for each simulation. It is the wave energy transmitted into the far‐field (Haskell, 1964). The viscoplastic simulation V1 radiates about 18% less seismic energy than the elastic simulation E1. Nonlinear radiation damping can achieve a similar reduction in seismic radiation; simulation R1 radiates about 24% less energy than E1. Simulation R2, with greater radiation damping, reduces the radiated energy by about 60%.

The situation is different when one looks at inelastic energy dissipation. For the viscoplastic simulation it is the amount of energy converted to heat by plastic yielding, whereas for the radiation damping simulations, it is the amount of energy converted to heat due to the radiation damping term in the modified friction law. The energy dissipated by the nonlinear radiation damping term in simulation R1 is only 12% of the energy dissipated by plastic yielding in run V1. This indicates that when attempting to emulate the effects of inelastic yielding one should not attempt to match the amount of energy dissipated; instead, one should look at the amount of energy radiated or the slip rates.

Fracture energy is the amount of work done at the fault in excess of the work done against the final sliding stress (Andrews, 2004). A simulation with inelastic yielding can have fracture energy that appears to increase with the rupture size if one counts the off‐fault inelastic energy dissipation as part of the fracture energy (Andrews, 1976; Andrews, 2004). Nonlinear radiation damping can produce a similar effect, creating an apparent increase in critical slip distance and fracture energy as the rupture propagates. Figure 3a shows shear stress plotted as a function of slip, for the same four locations on the fault illustrated in Figure 1, for simulation R1. The apparent critical slip distance is the amount of fault slip at which shear stress reaches its final sliding value, and the apparent fracture energy is the area under the curve and above the final sliding stress, and both are seen to increase with increasing rupture distance from the hypocenter. The actual slip‐weakening critical distance $Dc$ appearing in the friction law has the same value throughout the fault surface. Figure 3b shows the same plots for simulation R2. It is difficult to pick out an apparent critical slip distance for simulation R2, but the apparent fracture energy increases with rupture distance from the hypocenter.

Figure 4 illustrates how nonlinear radiation damping can reduce simulated ground motions. The figure shows particle velocity as a function of time, for two near‐fault stations on the Earth’s surface, at fault‐perpendicular distances of 1 and 3 km, respectively. In each plot, the elastic simulation E1 has the highest peak ground velocity (PGV), and simulations V1 (viscoplastic), R1, and R2 (both elastic with nonlinear radiation damping) have progressively lower PGV. There is considerable variation among the plots shown in Figure 4, but for simulations V1 and R1 the reductions in PGV appear generally consistent with the results of Roten et al. (2017), who studied the effects of plasticity on near‐fault PGV in dynamic rupture simulations. Simulation R2, which has a larger amount of nonlinear radiation damping, has a much greater reduction in PGV.

We have presented a new method for reducing slip rates and ground motions in dynamic rupture models with linear elastic rock properties. The method has only three parameters. It is simple to implement and efficient, requiring only that a nonlinear radiation damping term be added to the friction law. Despite its simplicity, the method can emulate some of the behaviors seen in viscoplastic models, including reduced peak slip rates, reduced radiated seismic energy, an apparent increase in critical slip distance with increasing distance from the hypocenter, and reduced PGV close to the fault.

Some experimentation may be needed to select parameter values, but a starting point can be obtained by examining fault slip rates and shear stress. $V0$ selects the slip rate at which radiation damping begins to take effect, so it is less than the highest peak slip rates. $Cr$ can be comparable to the stress drop (defined to be initial shear stress minus dynamic sliding shear stress) or strength excess (defined to be yield stress minus dynamic sliding shear stress) on the fault. If $Cr$ is much smaller, then the radiation damping will have little effect, whereas if $Cr$ is much larger, then the radiation damping will dominate the friction. For models in which the stresses vary with depth, $Cr$ (and possibly $V0$ as well) can vary with depth in the same manner as the stresses. The value n = 16 seems reasonable, although smaller values can be used to make radiation damping take effect more gradually. Because the slip rate rises sufficiently far above $V0$ to get past the “corner” in the radiation damping function, the precise value of n should not be too important.

The parameters used for our radiation damping simulations are consistent with the aforementioned guidance. Referring to panel V1 in Figure 2, we see that in our viscoplastic simulation, the peak slip rate begins to level off around 2 m/s and ultimately reaches about 3 m/s, suggesting that values of $V0$ in the range of 2 m/s would be reasonable. For our models, at a depth of 10 km, the stress drop is about 6 MPa, and the strength excess is about 10 MPa, suggesting that $Cr$ in this range would be reasonable.

For this demonstration, we assigned the parameters of simulation R1 to closely match the viscoplastic simulation V1. We used trial and error to adjust $V0$ and $Cr$ to produce peak slip rates like those in V1. Our trial‐and‐error procedure did not vary n. This resulted in simulation R1 having somewhat lower radiated seismic energy than simulation V1; it would have been possible to adjust $V0$ and $Cr$ to achieve a closer match in radiated seismic energy with a less‐close match in peak slip rates. For simulation R2, we arbitrarily reduced $V0$ and increased $Cr$ to show the effects of having a larger amount of radiation damping.

Nonlinear radiation damping works by dissipating energy at high slip rates. Our nonlinear radiation damping term can be thought of as emulating the reduction in shear traction on the fault caused by inelastic yielding in the fault zone, which is not otherwise present in a linear elastic model. In effect, we relocate energy dissipation from inelastic material near the fault to a nonlinear radiation damping term directly on the fault. Because nonlinear radiation damping operates on the fault, it does not emulate the effects of inelastic yielding that occurs in the far‐field; such yielding might contribute to the attenuation of seismic waves.

Another way to dissipate energy is to add viscoplasticity to the model, but doing so increases the model’s complexity and computational cost. Viscoplastic parameters are poorly known and only indirectly related to the behavior of the fault rupture. The absolute state of stress in the crust is also poorly known. Furthermore, the initial stress tensor must be in static equilibrium, and when the initial stress tensor is resolved onto the fault surface it must agree with the initial fault tractions. All of this makes it difficult to establish parameters and initial conditions, particularly in models with 3D nonplanar fault geometry or a 3D heterogeneous velocity model. In contrast, adding nonlinear radiation damping to an elastic model achieves many of the same results, without changing anything except for adding one term to the friction law.

There have been some previous attempts to achieve a similar goal. Aochi and Ide (2004) proposed the idea of implementing an increase in the slip‐weakening critical distance with increasing rupture distance from the hypocenter. Our new method achieves this result organically, as an output of the simulation, instead of artificially enforcing it before the simulation begins.

Another approach was suggested by Andrews (2005), who proposed imposing a maximum slip rate on the fault surface. This approach produces unusual slip‐time histories, in which the slip rate remains constant at the imposed maximum for some period of time. During that time, the operation of the friction law is inverted; normally the code computes the slip rate induced by an applied shear stress, but during that time the code must instead compute the shear stress that produces a desired slip rate. Our new method produces typical‐looking slip time histories, and the friction law operates normally.

We mention that our new method cannot reduce the strong ground motions created when a propagating rupture hits a hard barrier at the end of a fault. Viscoplastic yielding can reduce these ground motions by reducing the high stresses associated with a hard barrier. In an elastic simulation, a “soft” boundary that gently slows the rupture as it nears the end of the fault can achieve a similar effect, as we have done in our simulations. We have not investigated how to emulate the effects of inelastic yielding at geometrical complexities internal to a fault, such as branches or sharp bends.

We have implemented and tested nonlinear radiation damping for use with slip‐weakening friction. The method can easily be applied to models with other friction laws, such as rate‐and‐state friction, but we did not test the use of nonlinear radiation damping with other friction laws.

The numerical simulations shown in this article were performed using the FaultMod finite‐element code (Barall, 2009). They were run on the U.S. Geological Survey (USGS) Denali supercomputer (Falgout et al., 2022). Descriptions of dynamic rupture benchmarks TPV26v2 and TPV27v2 (Harris et al., 2018) are available at https://strike.scec.org/cvws/benchmark_descriptions.html (last accessed March 2023). The supplemental material to this article includes details of our test procedures.

The authors declare that they have no conflicts of interest.

Thanks to Shuo Ma, David Lockner, Nick Beeler, and Ralph Archuleta for rewarding conversations in the early stages of this work; to Brad Aagaard and Shuo Ma for thoughtful U.S. Geological Survey (USGS) internal reviews; and to Frantisek Gallovic and an anonymous reviewer for thoughtful journal reviews. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.