Knowledge of the shear stress on a fault during slip is necessary for a physically-based understanding of earthquakes. Borehole temperature measurements inside the fault zone immediately after an earthquake can record the energy dissipated by this stress. In the first Wenchuan Earthquake Fault Zone Scientific Drilling Project hole (Sichuan province, China) we repeatedly measured temperature profiles from 1.3 to 5.3 yr after the 12 May 2008, Mw 7.9 Wenchuan earthquake. The previously identified candidate for the principal slip surface had only a small local temperature increase of at most 0.02 °C with no obvious decay. The small amplitude of the temperature increase provides an upper bound for the frictional heat–generated coseismic slip, but is unlikely to be a frictionally generated signal. Two larger temperature anomalies are located above and within the fault zone. However, neither anomaly evolves as expected from a frictional transient. We conclude that the frictional heat from the Wenchuan earthquake remains elusive and the total heat generated at this location is much less than 29 MJ/m2. Low friction during slip is consistent with the temperature data.


The magnitude of the shear stress resisting slip along a fault during an earthquake has long been unknown. Recent measurements have shown that a straightforward application of Byerlee’s law with a coefficient of friction of 0.6 for most rocks and 0.2 for clays may not yield the correct shear stress at high speeds (Byerlee, 1978). At typical earthquake slip velocities of 1 m/s, laboratory values of friction plummet and span a range of values from 0.05 to 0.4, depending on lithology and experimental conditions (Di Toro et al., 2011). The theory controlling high-velocity friction is vigorously debated, and nonfrictional processes can alter the local shear stress in natural systems. Field measurements are needed to constrain the magnitude of shear resistance during earthquakes on actual faults.

To address this observational gap, rapid response drilling projects have measured the temperature in fault zones directly after major earthquakes. Single profiles in both shallow and deep boreholes after the 1999 Mw 7.7 Chi-Chi earthquake (Taiwan) recorded temperature anomalies interpreted to be equivalent to a dynamic coefficient of friction of ∼0.1 (Tanaka et al., 2006; Kano et al., 2006) and a similar value was inferred for the 2011 Mw 9.0 Tohoku earthquake (Japan) from the Japan Trench Fast Drilling Project (Fulton et al., 2013). All of these experiments tracked the temperature over <1 yr.

In this paper, we report on borehole temperature measurements made across the fault zone that ruptured during the 12 May 2008 Mw 7.9 Wenchuan (Sichuan province, China) earthquake, continuing from 1.3 to 5.3 yr after the event. We use the temperature data and thermal conductivity measurements to place an upper bound on the coseismic friction on the fault.


The Wenchuan Fault Earthquake Zone Scientific Drilling Project (WFSD) has successfully collected a data set over a 4 yr period with 23 temperature profiles in a single borehole. The measurement hole (WFSD-1; Fig. 1A) was drilled from November 2008 to July 2009 to a maximum depth of 1201 m (Fig. 1B) with a cased interval to a depth of ∼810 m.

The borehole intersects a major fault at 589 m borehole (core) depth, i.e., as measured along the hole, which is equivalent to 578 m vertical depth (all depths reported hereare borehole depths, unless otherwise specified). The 589 m fault is near the boundary between the Neoproterozoic Pengguan complex, which here consists of diorite, porphyrite, pyroclastics and other volcanics, and the Late Triassic Xujiahe Formation, which here consists of interbedded sandstone and shale (Li et al., 2013). The cataclasite, fault breccia, and fault gouge layers in the core continue to 759 m core depth (746 m vertical depth). Multiple faults exist throughout the 589–759 m zone, and the identification of a single principal slip surface corresponding to the 2008 Wenchuan earthquake is difficult. Togo et al. (2011b) reported that at a nearby surface outcrop, the surface break is localized on the Sichuan basin side of the fault zone packet, which corresponds to the 759 m depth in the core. However, internal evidence from the core suggests that the 589 m fault is the strongest candidate for the principal slip zone of the Wenchuan earthquake because of the fresh gouge appearance, microstructures relating to coseismic slip, high magnetic susceptibility values (Li et al., 2013), clay mineral composition in cores (Si et al., 2014), borehole logging data (Li et al., 2014), and drilling mud gas concentrations (Tang et al., 2013). Breaking both the surface rupture on the basin side and the 589 m fault at depth during the Wenchuan earthquake either requires a 65° dipping fault crossing the gouge packet (Li et al., 2013) or, more likely, slip transfer across the packet on to the subparallel en echelon surfaces.

Repeated temperature measurements were made in the well from October 2009 through September 2013 (Fig. 1C). For each of the 23 profiles, temperature was measured by lowering a string of two or three temperature sensors (RBR Ltd. model 1050/2050) at a rate of ∼1 m/min downward through the well to the bottom of the casing (800 m) and then raising the string at the same rate. Measurements were recorded every second and a built-in pressure sensor was intended to record depth. Only the downgoing data from the bottom sensor are interpreted because both the trailing sensors and upgoing runs are affected by the disturbance of the fluid in the well during temperature logging. For three profiles from June and December 2012, a stop-go logging technique was used with the sensor held stationary for at least 90 s to allow equilibration at fixed intervals that varied between 0.2 and 1 m. The data were then fit through the hold time to find the asymptotic steady-state temperature for comparison with data obtained with the sensors in motion (Harris and Chapman, 2007). The profile from 20 September 2011 was unusable due to a faulty pressure sensor. Profiles on 23 October 2009 and 23 June 2012 do not extend throughout the study zone and the 18 October 2010 profile was lowered too quickly for accurate temperature gradient estimates. The remaining 19 profiles were used for the following analysis.

We used the temperature gradient change associated with a change of lithology documented in the core and gamma ray logs at core depth of 394 m to align the logs (Figs. DR1 and DR2 in the GSA Data Repository1) because instrumental problems and variable densities in the muddy water in the borehole made the depth inferred from the pressure transducer inaccurate. The shifted 19 profiles overlay almost exactly, indicating that the data quality is good and the long-term geothermal gradient of ∼0.02 °C/m is stable (Fig. 1C). This inference of stability is independent of the corrections applied here.

The temperature measurements are interpreted in conjunction with measurements of thermal conductivity of the recovered core (see the Data Repository). The thermal conductivity was measured by an optical scanner at 5 m intervals over the core from 350 to 800 m depth with denser sampling near 589 m under dry conditions. Saturated thermal conductivities are estimated from a composite thermal model, as described in Data Repository equations DR1 and DR2.


Just below the previously identified candidate principal slip zone at 589 m core depth, a small 15-m-wide temperature deviation of 0.02 °C persisted throughout the observation period. This perturbation to the geotherm is most clearly resolved by the highest precision experiments that used stop-go logging (Fig. 2B). These high-precision logs also have the best depth control of any data collected. Similar amplitude temperature perturbations exist elsewhere in all profiles taken after the end of drilling. A temperature log taken during drilling reported a 0.15 °C anomaly at this same location (Li et al., 2014; Fig. 2); however, as the drilling was stuck at ∼590 m for over 30 days due to the difficulties of penetrating the fault zone, the early time data are strongly influenced by the drilling fluid temperature (see the Data Repository). For the data considered here, only the association with the strongest geological candidate for the principal slip zone distinguishes the 0.02 °C feature as worthy of further examination.

A larger feature in the temperature profile is centered at ∼690 m depth with a width of ∼130 m, which nearly spans the fault zone (Fig. 2A). The amplitude decays from 0.25 °C to 0.20 °C over the observation interval, but the width does not change. Above the primary fault zone, a feature with a similar width but smaller amplitude is centered at ∼450 m. All of these features exist in both the continuous and the stop-go logs, indicating that sensor reequilibration is not a factor in their identification.

We follow standard interpretive procedure and examine temperature as a function of thermal resistance to form a Bullard plot (Beardsmore and Cull, 2001). Thermal resistance incorporates variations in thermal conductivity and the procedure allows us to estimate the steady-state conductive heat flow. The mean resultant heat flow is 69 mW/m2 with a range of 68–70 mW/m2 over all the profiles (see the Data Repository). We remove the steady-state conductive temperature profile at each depth and define the residual temperature as the anomaly (Fig. 3A). The resultant anomalous temperature rise over the 589 m fault is still small and no anomaly >0.01 °C exists over the second candidate fault zone at 759 m. The wider features noted in the raw data are also anomalies in the residual temperature centered at 450 m and 690 m depth.


In our search for the frictionally generated heat, we first consider the largest thermal anomalies: the 450 and 690 m depth features. Although there is no well-developed fault plane documented at either of 450 or 690 m depth, there are fault breccias within each of the anomalies (Fig. 3 of Li et al., 2013).

If these thermal anomalies are frictionally generated, the residual temperature provides a direct constraint on the thermal energy S generated on the faults, i.e.,
where ρ is the density and cp is the heat capacity. Estimating cp = 800 J/m3 (Beardsmore and Cull, 2001) and ρ = 2500 kg/m3 (Li et al., 2014) results in S of 5 ± 2 MJ/m2 at 450 m and 24 ± 6 MJ/m2 at 690 m; error ranges indicate 1 standard deviation on the estimates over the suite of profiles.

To test the consistency of the data with a frictional model, we calculate the evolution of the temperature field over time assuming that the thermal energy is generated on fault planes at 450 and 690 m depth and diffuses into the surrounding rock (Fig. 3B). The diffusive model predicts that the anomalies should widen and decay resolvably in amplitude over time. Although some reduction is seen in the 450 m depth anomaly, neither observed anomaly widens over time and both are narrower than predicted at the end of the study period.

We conclude that the lack of observed time evolution implies that the 450 m and 690 m anomalies are unlikely to be due to diffusion of frictional heat away from a fault. These anomalies could be the result of frictional heating combined with another source of heat transport. For example, advective flow could generate a narrow thermal pulse limited by the high hydraulic diffusivity structure in the damage abutting the main fault zone (Fulton et al., 2010; Xue et al., 2013). A fluid flow explanation is qualitatively the most attractive possibility, and detailed modeling of such a fluid flow is the focus of future work. We confine ourselves here to concluding that the temperature data provide an upper bound of 29 MJ/m2 of frictionally dissipated energy. The heat energy could have been dissipated on any number of planes within the region; there is no direct constraint on the localization of the energy.

The geologically significant 589 m depth fault has a subtle temperature increase. However, the width of a frictionally generated anomaly is predicted to be much greater than observed. For a homogeneous medium with a thermal diffusivity of 1.5 × 10−6 m2/s and heat generation at the time of the earthquake, the anomaly is expected to be ∼100 m wide at the time of measurement, which is not observed (Fig. 2B). Furthermore, the signal does not decay or widen (Figs. DR2 and DR3).

The temperature increase at 589 m could be attributed to a gradient change associated with local thermal conductivity structure. We evaluate this possibility in Figure 4 by calculating the steady-state thermal profile consistent with the full suite of thermal conductivity measurements in the 589 m zone and a constant heat flow of 69 mW/m2 constrained by the analysis in Figure 3. The procedure is similar to that of Tanaka et al. (2007). Although the symmetry of the observed temperature increase is not reproduced by a realistic thermal conductivity structure, the measured variation in thermal conductivity suggests an even larger range of temperatures than observed. Here we use only the dry laboratory values and therefore Fig. 4 may over-estimate variations of thermal conductivity. In the model of Fig. DR4, we took the opposite approach and designed an inversion in which porosity completely compensates for thermal conductivity variations in the laboratory samples. Figure 4 is primarily useful as a warning of the potential effects of heterogeneous structure, but is unlikely to provide an accurate estimate of the residual temperature as saturation effects are not included. We have already concluded based on the width and lack of time dependence that the temperature increase below 589 m is unlikely to be frictionally generated, and now modeling shown in Figure 4 permits the feature to be the result of heterogeneous thermal conductivity structure.

The lack of a resolvable frictional heat anomaly yields an upper bound for frictional heat dissipated here during the earthquake. For frictional heating during slip, the decay of the peak temperature, T, on a planar fault due to one-dimensional heat conduction away from the fault is
where μ is the average coefficient of friction during slip, σn is the effective normal stress, d is the slip, cp is the specific heat capacity, ρ is the density, κ is the thermal diffusivity, and t is the time elapsed since the event (Carslaw and Jaeger, 1959). The model of Equation 2 is for an infinitesimally thin fault. The times of the data recorded here are sufficiently long after the earthquake that a thermal anomaly would have spread over a region larger than the shear width and the planar solution is indistinguishable from the finite shear zone solution. Based on Equation 2 and a normal stress on the fault equal to the lithostatic overburden less the hydrostatic pore pressure, the 0.02° C upper bound in the 589 m zone implies that the effective coefficient of friction during the earthquake is <0.02 if 7 m of slip happened on this surface (see Fig. DR7 caption for parameters). Although the bound in terms of the coefficient of friction enables easy comparison with previous work, dynamic reduction of normal stress or inaccurate estimate of slip are also acceptable explanations of the data. The data most directly constrain the dissipated heat energy, which over the depth of Figure 2B is <1.2 MJ/m2.

In summary, the total energy dissipated frictionally at the locale penetrated by the WFSD-1 borehole is <29 MJ/m2 and the dissipated energy on the previously identified fault surface is <1.2 MJ/m2. The upper bound for the entire zone (29 MJ/m2) is less than would be anticipated with a coefficient of friction of 0.6, but could be consistent with the range of coefficients of friction seen in dynamic weakening experiments (Di Toro et al., 2011). Other studies have inferred dramatic weakening for the Wenchuan fault material based on laboratory studies of fault zone rocks from surface outcrops subject to high-velocity shear (Togo et al., 2011a; Chen et al., 2013; Yao et al., 2013; Zhang and He, 2013). The high organic and clay content or presence of graphite may be critical factors (Zhang and He, 2013; Kuo et al., 2014). All of the cited studies concluded that the coefficient of friction during slip was <0.2, which would be consistent with the upper bound on the total dissipated energy. The Wenchuan earthquake fault appears to have been weak during slip.

This research was supported by the Chinese National Science and Technology Planning Project (Wenchuan Earthquake Fault Zone Scientific Drilling Project, WFSD), National Science Foundation of China grant 41330211 (to Li) and U.S. National Science Foundation grant EAR-1220642 (Brodsky). We thank W. Zhang, G. Yang, R. Guo, and Y. Huang for data collection assistance. Detailed reviews from T. Shimamoto, S. Nielsen and an anonymous reviewer are enormously appreciated.

1GSA Data Repository item 2015060, supplemental documentation of profile alignment, thermal conductivity measurement and interpretation, frictional heat model and drilling disturbance model, including Table DR1 and Figures DR1-DR7, as well as the original temperature data, is available online at www.geosociety.org/pubs/ft2015.htm, or on request from editing@geosociety.org or Documents Secretary, GSA, P.O. Box 9140, Boulder, CO 80301, USA.