Rupture Energetics in Crustal Rock From Laboratory‐Scale Seismic Tomography

The energy released during earthquake rupture is partly radiated as seismic waves and mostly dissipated by frictional heating on the fault interface and by off‐fault fracturing of surrounding host rock. Quantification of these individual components is crucial to understand the physics of rupture. We use a quasi‐static rock fracture experiment combined with a novel seismic tomography method to quantify the contribution of off‐fault fracturing to the energy budget of a rupture and find that this contribution is around 3% of the total energy budget and 10% of the fracture energy Gc. The off‐fault dissipated energy changes the physical properties of the rock at the early stages of rupture, illustrated by the 50% drop in elastic moduli of the rock near the fault, and thus is expected to greatly influence later stages of rupture and slip. These constraints are a unique benchmark for calibration of dynamic rupture models.


Introduction
Strain energy released during earthquakes is partly radiated as seismic waves that cause ground shaking and mostly dissipated by frictional heating on the fault interface and by fracturing of the rocks surrounding the fault. The latter energy sink, and a component of frictional heating, constitutes the fracture energy (G c , sometimes referred to as rupture energy) that dictates the dynamics of rupture propagation (Rice, 1980). Quantification of these individual components is crucial to understand the physics of rupture, to better understand the feedback between rupture and slip processes, and to improve ground motion predictions.
Fracture energy is the work associated with the breakdown of the rock strength toward its residual frictional strength. G c is a collective term for several dissipative processes in the breakdown zone around the rupture tip, both on and off fault (Kanamori & Rivera, 2006). These dissipative processes may include shear heating, plastic yielding, on-and off-fault creation of fractures (surface energy), and grain comminution. A measure of G c can be inferred from earthquake data (Tinti et al., 2005), and from laboratory mechanical (Nielsen et al., 2016;Wong, 1982Wong, , 1986 and acoustic data (Lockner et al., 1991), but such data do not provide a breakdown of the individual components of G c . Cumulated surface energy measured on fault rocks (Chester et al., 2005;Faulkner et al., 2011;Rockwell et al., 2009;Savage & Brodsky, 2011;Wilson et al., 2005) provide onand off-fault components of G c , but these estimates are measured on rocks that have recorded numerous earthquakes and deformation episodes and therefore do not represent a single earthquake, nor can they constrain energy dissipation into heat. Cumulated surface energy obtained from microstructural studies on off-fault damage in laboratory samples (Moore & Lockner, 1995;Zang et al., 2000) are only static snapshots of the dynamic breakdown process. To establish the off-fault energy dissipation component G off , one possibility is to use the change in off-fault elastic properties caused by off-fault deformation. Such changes must be measured in situ during rupture, ideally under realistic crustal conditions (i.e., at elevated pressure and temperature). This may be done by the acquisition of active seismic surveys during laboratory rupture experiments (Lockner et al., 1977). Yet, the size and geometry of the off-fault damage zone and the actual local wave speeds therein remain unconstrained because of the lack of spatial resolution of conventional laboratory ultrasonic measurements.
Here, we combine stress, strain, and acoustic emission (AE) and ultrasonic velocity measurements obtained in situ during a laboratory rock fracture experiment. From this, we determine the time-resolved 3-D seismic velocity structure of a growing fault zone that provides the size and geometry of the off-fault elastic properties. Taken together, our measurements allow us to estimate the partitioning of G c into off-fault (G off ) and on-fault energy dissipation.

Method
We performed a triaxial rupture experiment on a 40 mm diameter, 100 mm length sample of dry Lanhélin granite (Brittany, France). The sample was placed into a rubber jacket equipped with 16 piezoelectric P wave transducers and two pairs of axial-radial strain gauges ( Figure 1). Acoustic signals were recorded by a digital oscilloscope at a 50 MHz sampling frequency, after being amplified to 40 dB. Active ultrasonic velocity surveys were performed every 5 min by sending a 1 MHz pulse at a voltage of 250 V to one transducer, while the other transducers recorded the resulting waveforms. During one survey, all 16 transducers were used as a source, and the results of six pulses for each source transducer were stacked to enhance the signal-to-noise ratio. In between the acoustic velocity surveys, the waveforms of AEs were recorded provided that a signal amplitude of 250 mV was surpassed on at least two transducers. All waveforms, both of active acoustic velocity surveys and AEs, consisted of 4,096 datapoints (82 μs length).
The jacketed sample was placed into a triaxial deformation rig and pressurized to 100 MPa confining pressure. Axial stress was measured by a load cell; axial shortening was measured by a pair of linear variable differential transducers. The axial deformation was then applied by a piston that moved with a strain rate of 10 −5 s −1 for the elastic portion of the stress-strain curve and 10 −6 s −1 for the remainder of the experiment. The axial shortening rate was controlled in such a way as to hamper the dynamic propagation of shear rupture, using a technique similar to that of Lockner et al. (1991): The AE rate was monitored visually and, when the AE rate showed acceleration (about eight hits or more per second, recorded on at least two channels), the direction of movement of the piston was reversed to reduce the load. More than hundred of such load reductions were performed. The overall fracture propagation across the sample occurred over a time interval of around 16 hr.
From the ultrasonic data set, we computed the AE source locations together with the evolution of the seismic velocity structure within the sample by using the 3-D seismic tomography code FaATSO, specifically designed for laboratory rock deformation experiments (Brantut, 2018). Twelve thousand of the highest quality AE events were considered for the tomography. For all these AE events; the first P wave arrivals were picked using an automated picking algorithm. The automated picks were manually checked and, where necessary, improved or cancelled. The acoustic data were split into 38 time intervals of 100-400 AEs that cover a more or less equidistant interval of the stress-strain curve ( Figure 2). The mean posterior model determined by FaATSO yielded the V P structure during each time interval. The spatial resolution of the velocity structure is 5 mm, which is similar to the wavelength of the 1 MHz surveys. The velocity structure of the previous time interval was used as the a priori velocity structure for the next time interval. For the a priori velocity structure of the first time interval, we took a transverse isotropic structure derived from the first active ultrasonic surveys. The inverse problem was constrained by several inversion parameters that ascribe Gaussian variances to the input data, expressed as standard deviations (Table 1). The velocity model was ascribed a covariance as a function of the variance of the velocity and a correlation length. Heterogeneities and velocity contrasts smaller than the correlation length are suppressed and smoothed, eliminating unrealistic sharp velocity contrast from the results. Such a covariance was ascribed to the anisotropy structure as well.
The AE source locations were recomputed by using the 3-D seismic velocity structure from FaATSO. An approximation of the position of the rupture front and the seismically active fault surface for each time interval were obtained from the relocalized AE locations ( Figure S1 and Text S1).

Results
Before the peak stress and the onset faulting, we observe an overall decrease in V P from around 6 down to 5 km/s ( Figure 2a). Then, rupture starts at the bottom of the sample and propagates upward ( Figure S1 and Movie S1). A low-velocity zone develops parallel to the rupture plane and migrates along with the growing fault (delineated by the AE source locations, Figures 2b, 2c, and S1). Velocities in the localized zone are as low as 4.6 km/s-a 25% drop relative to the areas outside of the fault zone where V P remains nearly constant ( Figure 2c). This corresponds to a drop of around 50% in P wave modulus. We interpret this low-velocity zone as the fault damage zone, which is generated ahead and along the propagating rupture tip. In the wake Correlation length 25 mm of the rupture tip the damage zone width decreases slightly by several millimeters. There is a widespread partial recovery of V P throughout the damage zone (the minimum value rising from 4.6 to about 4.7 km/s) at the onset of the frictional sliding stage (Figure 2d). The V P anisotropy at the peak stress is 13% (i.e., vertical V P is 13% higher than horizontal V P , Figure S3) and increases during fault growth up to 20% in two large zones adjacent to the fault. In the wake of the passing fault tip, the change in velocity anisotropy is much less than the change in velocity (Figure 1). The anisotropy decreases again as the axial stress is reduced during the frictional sliding stage.
The robustness of the tomographic inversion results is demonstrated by the mean posterior velocity and 500 individual posterior solutions along a fault-perpendicular transect (Figures 3a-3d). The individual solutions (grey

Rupture Energetics
Now, we can determine the fracture energy G c , the total dissipated energy, and the off-fault dissipated energy G off . G c and the total dissipated energy are computed from the shear stress versus fault slip record (Wong, 1982) up to the slip-weakening distance 0 (Figure 4). We assume that all axial shortening is caused by slip along the fault from localization onward. Axial shortening is corrected for elastic strain by using the intact elastic moduli for the rock. G c is 27 kJ m −2 , similar to previous experimental results (13-29 kJ m −2 ; Lockner et al., 1991;Wong, 1982Wong, , 1986).
Next, we estimate the off-fault dissipated energy G off during the slip-weakening stage. G off is given by the change in stored elastic strain energy (i.e., elastic softening) around the fault interface. These elastic compliance changes are caused by off-fault dissipative processes, of which microcracking is dominant. Strain derived from mechanical data includes slip along the fault and thus cannot be used to obtain the off-fault strain components necessary to obtain changes in elastic compliance. Instead, the spatial and temporal evolution of the seismic velocity structure is used to obtain these. For each time interval, G off is approximated as The ratio between the total fracture energy G c (solid black curve) and the off-fault energy dissipation G off (dashed black curve) shows a large initial off-fault contribution (brown curve). The slip-weakening distance ( 0 ) is the interval between onset of rupture and convergence of the shear stress toward the frictional residual strength as shown in the shear stress versus slip plot below. 0 is the slip at which the entire fault interface has formed, which is evidenced by the distribution of acoustic emission source locations. For the total slip-weakening distance, the relative energy contribution is around 10%.
wherēi are the average stress components between two time intervals, ΔS ijkl is the change in the elastic compliance tensor between two time intervals, and w is the width of the damage zone. S ijkl is estimated from the seismic velocities, where we assumed that microcracks were oriented parallel to the loading direction (see Text S2). Stress rotations within the damage zone were neglected for simplicity, and for each time interval a single value for ΔS ijkl represented the entire damage zone (see Text S2).
A first-order approximation of w is established by analyzing the spatial extent of permanent damage in our data. Nearly elastic behavior (i.e., full recovery of V P ) is observed at some distance from the fault interface after passing of the rupture front (point p 3 , Figures 3e and 3f). This matches qualitatively with the predictions of a stress field around a passing rupture tip (Freund, 1990), whereby the stiffness reduction and the full recovery that follows is caused by the rupture tip stress field that dilates or contracts preexisting flaws. Closer to the fault interface such full recovery is not observed (points p 1 and p 2 , Figure 3f). Within this zone of at most 20 mm in width, permanent microcrack damage is generated. We adopt a more conservative damage zone width w of 10 mm, since the AE source locations are clustered in a more narrow band of 5 to 10 mm around the fault (Figure 2; Zang et al., 2000).
G off is around 3 kJ/m 2 and is a similar order of magnitude to G off established from postmortem microstructures (Moore & Lockner, 1995). G off is approximately 10% of the fracture energy ( Figure 4) and about 3% of the total dissipated energy during the fracture process. Hence, 90% of G c is dissipated by on-fault processes. G off accounts for over 15% of G c for the earliest stage of rupture (up to 0.1 mm of slip, Figure 4).

Discussion
The constraints on rupture energetics were obtained during quasi-static rupture propagation and are representative of the nucleation phase of an earthquake. Once the rupture velocity has accelerated toward the critical wave speed of the rock (equal to the Rayleigh wave speed in in-plane conditions, or to the shear wave speed in antiplane conditions), the stress field around the rupture tip is distorted relative to that of a quasi-static rupture tip (Freund, 1979). Such a distortion is more likely to increase the size of the off-fault regions where fracturing might occur (Poliakov et al., 2002;Rice et al., 2005). In addition, the transient rupture tip stress field imposes high strain rates on the off-fault rock volume, which results in increased fragmentation (Bhat et al., 2012;Grady, 1982). G off is thus expected to increase with increasing rupture velocity, changing the ratio G off over total energy dissipation. Therefore, the quasi-static ratio of 3% estimated from our experimental data is a lower bound for dynamic ruptures but provides a unique calibration benchmark for dynamic rupture models that allow for off-fault damage (Thomas & Bhat, 2018;Xu et al., 2015). Such models predict a maximum drop in V P of around 30% (Thomas & Bhat, 2018;Xu et al., 2015), which is consistent with the maximum drop of 25% observed here. This maximum drop in V P of 25% during experimental rupture is of a similar order to geophysical observations on coseismic V P reduction near recently ruptured faults (Allam & Ben-Zion, 2012;Cochran et al., 2009;Froment et al., 2014). However, these velocity reductions were the product of multiple ruptures with a higher rupture velocity than our experimental rupture, and they also reflect the postseismic state rather than the coseismic state that we document here-thus without the transient reduction of elastic properties. Geophysical observations of seismic velocity reductions caused by single-rupture events are of the order of 20-45% but are restricted to S wave velocity only (Karabulut & Bouchon, 2007;Wu et al., 2009). Rather than directly comparing the absolute values observed here with those measured on faults, our constraints on rupture energetics of laboratory-sized samples can be upscaled to larger faults by relying on scaling relations established by other studies.
First, studies along exhumed faults suggest that damage zone width increases linearly with total fault displacement (Faulkner et al., 2011;Savage & Brodsky, 2011), and total fault displacement is linearly proportional to fault length (Cowie & Scholz, 1992). Earthquake fault slip increases linearly with fault length (Scholz, 1982). From this it follows that damage zone width scales linearly with fault slip (Savage & Brodsky, 2011;Faulkner et al., 2011). G off will increase with increasing damage zone width, assuming that the drop in elastic moduli is independent of (i.e., the off-fault drop in elastic moduli occurs only during the slip-weakening phase). This implies that G off ∝ , and the ratio of G off to total energy dissipation, which is 3% for our experiment, thus remains constant. Field estimates of this ratio are of the order of 1% (Chester et al., 2005;Rockwell et al., 2009), which further supports this conclusion.
Second, seismological estimates and theoretical predictions indicate that fracture energy G c scales with fault slip with an exponent (Brantut & Viesca, 2017;Viesca & Garagash, 2015). By adopting aforementioned linear scaling of G off with , we obtain G off ∕G c ∝ 1− . The exponent ≈ 2 for small earthquakes (slip less than around 10 cm, M w ≲ 4) and < 1 for larger events (Brantut & Viesca, 2017;Viesca & Garagash, 2015). Thus, we expect the ratio G off ∕G c (for which we measure G off /G c ≈ 10%) to decrease initially with increasing slip and earthquake magnitude and subsequently to stabilize or slightly increase with earthquake slip and magnitude.
These scaling relations are valid for total fault displacements up to a kilometer, because at larger displacements the extent of the damage zone is not proportional to displacement anymore (Savage & Brodsky, 2011) but depends on the width of the fault (which equals the seismogenic depth for strike-slip faults; Ampuero & Mao, 2017;Scholz, 2019). Similarly, there may be a crossover from being proportional to fault length at smaller fault lengths, to being proportional to the fault width (Scholz, 1994). In that case, our scaling relations for G c are valid for larger earthquakes as well.
Our experimental results have been obtained in an initially intact material with a large cohesion, whereas many earthquakes occur along preexisting faults, possibly containing clay-rich gouge. In that case, the fault is likely to have a lower peak strength and a shorter slip-weakening distance (Ohnaka, 2003), which decreases the size of the rupture tip processes zone and reduces the amount of damage. The off-fault dissipated energy for ruptures along natural faults is thus expected to be lower than for rupture in intact material. In terms of the ratio G off ∕G c , a rigorous estimate for natural faults should include the nature of the fault zone material and consider the roughness of faults. Fault roughness affects both G c (Ohnaka, 2003) and G off (Johri et al., 2014), and experiments conducted on initially intact materials (generating a spontaneous fault roughness during fault growth) provide a useful benchmark for more advanced numerical models.
The data presented in this study provide new insight in rupture processes, particularly the strong influence of the transient in situ stress state during rupture on the seismic velocities and elastic moduli both inside and outside the fault damage zone (Figure 2). We infer moduli variations even in the original host rock, due to the presence of preexisting flaws. The effect of damage on the seismic velocity structure of faults is hence strongly coupled to the local stresses around them. Therefore, observations from structural analysis in terms of fault rock microcracking are not expected to match necessarily the in situ elastic moduli distribution (and anisotropy) under realistic crustal stress states.
The damage-induced reduction in elastic moduli around the propagating rupture tip (Figure 2) suggests local dilatancy of the rock. When pore fluids are present, dilatancy causes a drop in pore fluid pressure that stabilizes brittle failure (Martin, 1980). The pore fluid pressure drop may decrease the efficiency of fluid-driven slip-weakening processes such as thermal pressurization (Lachenbruch, 1980), but its effect can only be assessed fully when the changes in local pore fluid pressure, permeability, and storage capacity of the rock near the fault interface are known (Brantut, 2019).
A substantial component of the seismic waves radiated from the rupture tip process zone may be caused by off-fault reduction of elastic moduli in addition to radiation from a classic planar rupture (Ben-Zion & Ampuero, 2009). Theoretical first-order magnitude estimates of this off-fault component are based on an arbitrary off-fault drop in prerupture to postrupture stiffness of 50%, assuming isotropic elasticity (Ben-Zion & Ampuero, 2009). Here, we verify that the P wave modulus indeed drops by 50%, and we provide the temporal and spatial evolution of the in situ stiffness matrix including anisotropy. Thus, our results can help to improve estimates of earthquake source properties and predictions of strong ground motion caused by seismic radiation.
The laboratory-scale seismic tomography method used here provides a unique constraint on the elastic properties of propagating faults and on the energetics of rupture. We show a significant drop of elastic parameters