Frictional heating in a thick fault zone with empirical slip-weakening friction: implications for slip parameter estimation from temperature observations in deep fault drilling

The strain energy released during an earthquake is consumed by processes related to seismic radiation or dissipation. Deep fault drilling and subsequent temperature measurements in a thick fault zone immediately after an event have provided important insights into this dissipation process. By employing an analytical solution to the heat conduction problem, which involves the sudden injection of an infinitesimally thin heat source into an infinite medium, previous drilling projects have estimated the strength of the heat source and the level of shear stress from observed temperature anomalies. However, it is unclear under what conditions this analytical source solution can be regarded as a good approximation for the thick fault problem, a situation which has led to uncertainty of the approximation error in these previous studies. In this study, I first derived an analytical solution for the thick fault problem that accounted for experimentally derived slip-weakening friction. I then validated the derived solution both analytically and numerically. Using the derived thick solution, I next demonstrated that the thick, planar, and source solutions can be considered equivalent under the typical conditions of the previous drilling projects. Therefore, the slip parameters estimated by using the source solution obtained by these studies are appropriate. These results suggest that coseismic information with spatio-temporal extent, such as shear stress and friction coefficient, are lost due to heat diffusion when the temperature observations are conducted; thus, they cannot be inferred directly from observed temperature anomalies. These results also suggest that for most drilling projects, including future ones, the observed temperature distribution can be well explained by using the source solution instead of the thick solution as long as coseismic slip is not markedly delocalized and the spatial extent of the temperature measurements is not significantly larger than the diffusion length.


Introduction
Strain energy stored in the Earth is released when earthquakes occur. The two principal sinks for the released strain energy are radiation (seismic waves) and dissipation (work done by friction). Whereas the former can be studied by examining earthquake source parameters (e.g., Abercrombie and Rice 2005), the latter is difficult to study using a seismological approach because the level of absolute stress must be known. Laboratory experiments have shown that most of the dissipation energy is likely to be associated with frictional heat (e.g., Lockner and Okubo 1983). This suggests that precise estimates of frictional heating in fault zones should lead to better understanding of dissipation process and, hence, earthquake energetics.
In general, there are two representative approaches to the examination of frictional heating in fault zones. In the first approach, anomalies in the temperature indices of fault rocks exposed at the Earth's surface have been investigated to detect heat signatures recorded in natural faults (Rowe and Griffith 2015). Although this approach is applicable to various natural faults, the lack of a modern seismological network at the time the corresponding event took place makes it impossible to compare the obtained information with the event's seismic source parameters. The second approach is to directly measure residual heat in post-seismic fault zones immediately after an event occurs. Some drilling projects have succeeded in detecting residual heat in a targeted active fault, such as the Chelungpu fault in Taiwan (Kano et al. 2006;Tanaka et al. 2006), which slipped during the 1999 Chi-Chi earthquake, and the plate-boundary fault in the Japan Trench (Fulton et al. 2013), which slipped during the 2011 Tohoku-Oki earthquake. However, when either advective fluid flow or spatial variation in thermal properties are efficient, it may be difficult to directly detect residual heat by temperature measurements (Li et al. 2015). Nevertheless, post-seismic temperature observations are of great importance because they provide the only practical way to quantitatively examine the dissipation process of one particular event in conjunction with seismological observations.
Observed anomalies in post-seismic fault temperatures can be used to constrain dissipation-related slip parameters, such as the strength of the heat source, by solving the heat conduction problem in and around a fault zone. Because numerical solutions always contain numerical errors and because investigations of a large number of parameter sets or long-term evolution may require large computational resources, analytical solutions, if available, are preferable, because they can cause numerical errors to vanish and reduce computational costs. Therefore, to estimate the strength of the heat source and the absolute level of shear stress, previous drilling projects have adopted an analytical solution for a one-dimensional system in which a planar heat source is suddenly injected into an infinite medium (Kano et al. 2006;Tanaka et al. 2006;Fulton et al. 2013;Li et al. 2015). Although these studies have assumed that it is reasonable to apply this analytical source solution to the heat conduction problem with a thick fault zone under their observation conditions, they provide no quantitative discussion of this point. In fact, seismic slip rates are usually expected to be no less than 10 -2 to 10 -1 m s -1 , which is sufficiently fast to cause a dramatic decrease in the shear stress, as has been demonstrated by numerous laboratory friction experiments (Di Toro et al. 2011). Therefore, an analytical solution to the thick fault problem that takes account of experimentally obtained slip-weakening friction is required to provide a quantitative criterion for using the source solution as a reasonable approximation. Although Carslaw andJaeger (1959) and, subsequently, Lachenbruch (1986) and Fulton and Harris (2012) have reported analytical solutions to the thick fault problem in a one-dimensional system with constant shear stress, no analytical solution accounting for typical slip-weakening friction has yet been reported. Therefore, it has not been possible to examine whether the source solution adopted by various previous drilling projects ( Kano et al. 2006;Tanaka et al. 2006;Fulton et al. 2013;Li et al. 2015) is reasonably applicable under the observation conditions of those projects.
In light of these backgrounds, the main objective of this study is to investigate the validity of the slip parameters reported in previous drilling projects using source solutions. I first derive analytical solutions for the one-dimensional thermal conduction problem with a thick fault zone by considering typical slip-weakening friction at seismic slip rates. I then verify the derived solutions by comparing them with previous analytical solutions, as well as with numerical solutions. Next, on the basis of the derived solutions, I investigate whether the source solution adopted by previous drilling projects as an approximation is appropriate and examine the validity of the reported slip parameters. Finally, on the basis of the obtained results, I discuss some implications of estimating coseismic slip parameters in deep drilling projects.

Methods
This section is devoted to describe the problem setup that will be addressed throughout this paper. I start from the general heat conduction equation in a three-dimensional, homogeneous, and infinite solid whose thermal properties have no spatial or temperature dependences: where is the Laplacian, T (x, t) is temperature, x is the position vector, t is time, and α is thermal diffusivity. For simplicity, I consider only the case in which heat conduction is the principal process governing the temperature change and the effects of other possible heat transfer mechanisms, such as fluid advection, are neglected. If the heat flow is linear and parallel to the x direction, Eq. (1) becomes Page 3 of 13 Kaneki Progress in Earth and Planetary Science (2021) 8:67 Under the assumption that the center of the fault zone is x = 0 , the system is symmetric with respect to the x = 0 plane and it is only necessary to solve the heat conduction problem in the semi-infinite region x ≥ 0 . If it is assumed that a fault starts to slip and to supply heat to the fault zone at t = 0 , then it is necessary to solve with the initial condition where Q(x, t) is the rate of heat production per unit volume, K is thermal conductivity, and T 0 (x) is the initial temperature distribution. The system discussed in this study is shown schematically in Fig. 1. I consider the case in which Q is related only to frictional heat. If deformation is homogeneous and confined to the fault zone, then Q can be expressed as where Q s (t) is the rate of heat generation by slip per unit volume, H(x) is the Heaviside function, and l is the halfwidth of the fault zone. If it is assumed that all of the work done by friction is converted into heat, then where τ (t) and v(t) are the shear stress and slip rate, respectively. Because work done by friction can also be consumed by other processes such as gouge fracturing and chemical reactions, this assumption provides an (3) where τ ss and τ p are the steady-state and peak shear stresses, respectively, d = v s t is the nominal slip displacement, and d c is the characteristic length (Fig. 2). Applying Eqs. (6)-(8) to (5), one obtains where c 1 = v s τ ss /2l , c 2 = v s τ p − τ ss /2l , and c 3 = v s /d c . The temperature distribution before slip occurs is assumed to be spatially homogeneous; thus, If temperature and heat flux are postulated to be spatially continuous, it follows that where q(x, t) = −K (∂T /∂x) is the heat flux per unit area. Assuming no effect of frictional heating when x → ∞ , one obtains another boundary condition, Therefore, I focus here on solving the following partial differential equations for T (x, t): with the initial and boundary conditions described by Eqs. (10) and (11).

Analytical solutions to the problem
Denote the Laplace transform of some function F (X, t) by F (X, s) . Then, it follows from Eqs. (10) and (12) that and from Eq. (11) that (13) can be solved for T in terms of x as follows: where A n (n = 1, 2, 3, 4) is a constant that is constrained by the boundary conditions (Eq. 14) as follows: Applying Eqs. (16) and (17) to (15), one obtains .
where erfc(x) is the complementary error function, i 2 erfc(x) is the second integral of erfc(x) , erfcx(x) is the scaled complementary error function, and i is the imaginary unit. By combining Eqs. (19) and (20), I finally arrive at analytical solutions for the problem of one-dimensional heat conduction with a thick fault zone ( Fig. 1) that account for empirical slip-weakening friction, as described in Eq. (8). Note that these solutions can also be applied to the semi-infinite region x ≤ 0 by replacing x with −x because the system is symmetric with respect to the x = 0 plane. The derived analytical solutions are useful for estimating slip parameters in fault zones that record frictional heating due to ancient seismic events. However, since the main target of this study is temperature measurements by deep fault drilling, I will not discuss this issue in depth.

Comparison with previous analytical solutions for the thick fault problem
To validate the derived solutions (19) and (20), I compare them with previously developed analytical  (21), describe the coseismic spatio-temporal temperature distribution and are equivalent to Eqs. (9) and (10) in Carslaw and Jaeger (1959, p. 80). Therefore, the solutions (21)

Comparison with analytical solutions for the planar fault problem
It is reasonable to expect that solutions (19) and (20) for the thick fault problem must approach the solution for the planar fault problem if the fault is sufficiently thin. The solution for the planar fault problem, which imposes a time-dependent heat flow ]v s τ (t)/2 on the x = 0 plane in a semi-infinite solid, can be expressed as follows (Additional file 1: Text S2): where ierfc(x) is the first integral of erfc(x) , c 4 = v s τ ss /2 , and c 5 = v s τ p − τ ss /2 . As in the case of the thick fault problem, because of system symmetry, this solution can also be applied to the semi-infinite region x ≤ 0 by replacing x with −x.
Here I consider the solutions for the thick fault problem under the conditions that l/ 2 √ αt ≪ 1 when t s ≥ t ≥ 0 and l/ 2 √ αt < l/ 2 √ α(t − t s ) ≪ 1 when t > t s . By performing the Taylor expansion of the functions for values of these parameters around zero and ignoring terms with degree higher than two, the thick solutions (19) and (20) can be approximated as follows (Additional file 1: Text S3): This approximate Eq. (24) is identical to the planar solution (23); therefore, provided that the fault thickness is sufficiently smaller than the diffusion length, the onedimensional heat conduction problem for a thick fault zone is identical to that for a planar fault zone when empirical slip-weakening friction (Eq. 8) is considered.

Comparison with numerical solutions for the thick fault problem
To further investigate the validity of the derived solutions (19) and (20), I compare them with numerical solutions in this subsection. As a typical example, I focus on the solutions for the Chelungpu fault, where thermal diffusivity α = 5.0 × 10 −7 m 2 s −1 , specific heat capacity c = 1300 J kg −1 K −1 , bulk rock density ρ = 2200 kg m −3 , initial temperature T 0 = 46.5 • C , the half-width of the fault zone l = 1 × 10 −2 m , slip duration t s = 6 s , slip displacement is 8.3 m, and deformation depth is 1111 m (Ma et al. 2003(Ma et al. , 2006Hirono et al. 2006;Kano et al. 2006;Hirono and Hamada 2010). Gravitational acceleration g and fluid density ρ f are set to 9.8 m s −2 and 1000 kg m −3 , respectively. Thermal conductivity K ≡ αρc is calculated to be 1.43 W m −1 K −1 . Effective pressure σ e under hydrostatic pore-fluid pressure condition is At a depth of 1111 m, σ e ≈ 13 MPa . If σ e is assumed to be constant during slip, it follows from Eq. (8) that where f ss and f p are the steady-state and peak friction coefficients, respectively. Here, the values of f ss , f p , and d c are set to 0.2, 1.0, and 2.5 m, respectively, based upon the experimental results of Tanikawa and Shimamoto (2009). The parameters used and their values are also summarized in Table 1. Figure 3a shows the analytical solutions evaluated at t = 0 , 6 , and 60 s with Eqs. (19) and (20). To calculate erfc(x) and erfcx(x) , I used the open-source Python functions scipy.special.erfc and scipy.special. erfcx (SciPy v1.6.2). To obtain the numerical solutions, a system with width 1.8 × 10 −1 m was discretized using N = 297 elements (corresponding to an element size of x = 6.06 × 10 −4 m ), based on the central difference scheme and postulating Dirichlet boundary conditions at both edges of the system. Equation (12) was then solved numerically by the explicit Runge-Kutta method of order 5(4) using the open-source Python function scipy.integrate.solve_ivp (SciPy v1.6.2) with absolute and relative tolerances of 10 -10 . Figure 3b shows the obtained numerical solutions evaluated at t = 0 , 6 , and 60 s , which are in good agreement with the analytical solutions ( Fig. 3a) with numerical errors | T | of up to ∼ 3 • C (Fig. 3c). It is reasonable to expect that the numerical error | T | must decrease as the element size x decreases if the derived solutions (19) and (20) are correct. Figure 4 shows the | T | values for various x values evaluated at (x, t) = (0 m, 6 s) on a log-log plot. It is clear that the numerical solutions show almost second-order accuracy in space, which is generally regarded as a sufficient rate of convergence. These findings also strongly support the validity of the derived solutions (19) and (20).

The approximate equation used by previous fault drilling projects
As  (19) and (20) were used to calculate the analytical solutions. The numerical solutions were obtained by the explicit Runge-Kutta method of order 5(4) (see main text for details). All of the parameters used to obtain the solutions are listed in Table 1 Element size,   Convergence test for numerically solving Eq. (12). The numerical error | T | was evaluated at (x, t) = (0 m, 6 s) and plotted against element size x . The numerical solutions were obtained by the explicit Runge-Kutta method of order 5(4) (see main text for details). All parameters, except for x , are the same as those used the numerical solution shown in Fig. 3b and listed in Table 1 Page 9 of 13 Kaneki Progress in Earth andPlanetary Science (2021) 8:67 et al. 2013;Li et al. 2015). All of these previous studies employed the same analytical solution to a one-dimensional system in which a planar heat source is suddenly injected into an infinite medium at (x, t) = (0, 0) (Eq. 4 in Carslaw and Jaeger 1959, p. 259): where Sρc is the heat liberated at the x = 0 plane per unit area. For the case in which shear stress shows exponential decay with slip (Fig. 2), it follows from Eq. (8) that Consider only the case when τ p = τ ss = τ leads to S = ατ v s t s /K . This source solution has been considered to be a reasonable approximation to explain the temperature distribution observed in post-seismic temperature measurements because the diffusion length is much greater than the width of the fault zone (Materials and Methods section in Fulton et al. 2013). However, as mentioned in the Introduction, no quantitative criterion has been reported for the use of the source solution (27), (28) as an approximate solution for the thick fault problem. This subsection focuses on this point and examines whether the conditions of the previous drilling projects ( Kano et al. 2006;Tanaka et al. 2006;Fulton et al. 2013;Li et al. 2015) were appropriate for the use of this approximation. In the case of the Chelungpu fault, I focus on the study by Kano et al. (2006) because in their targeted fault zone, information is available on all of the physical properties that must be known to address this issue. First, I examine the conditions under which the planar solution (23) can be a good approximation of the thick solutions (19) and (20). Thermometers are not installed in boreholes until after an event; thus, t > t s , and it follows from Eq. (24) that the thick solution T thick (x, t) approaches the planar solution T plane (x, t) when Eqs. (19), (20), and (23) shows that there are five nondimensional parameters that govern the ratio between T thick (x, t) − T 0 and T plane (x, t) − T 0 : namely, and c 2 /c 1 = c 5 /c 4 . Figure 5 shows the difference between the thick and planar solutions in relation to the nondimensional parameter l/2 √ α(t − t s ) for various combinations of x/2 √ α(t − t s ) , t s /t , c 3 t , and c 5 /c 4 . Previous drilling projects examined temperature data in the spatial range Kano et al. 2006;Fulton et al. 2013;Li et al. 2015). I thus set the value of x/2 √ α(t − t s ) to either 0 or 5 as a typical minimum or maximum value, respectively. The values of the other nondimensional parameters can be roughly constrained by using the typical conditions for a seismic event and post-seismic temperature measurements: t s ≈ 10 0 − 10 2 s , t ≈ 10 7 − 10 8 s , c 3 = V s /d c ≈ 10 −2 − 10 3 s −1 , and f p /f ss ≈ 1 − 11 , which yield t s /t ≈ 10 −8 − 10 −5 , c 3 t ≈ 10 5 − 10 11 , and c 5 /c 4 ≈ 0 − 10 . Note that when c 5 /c 4 = 0 , the shear stress remains constant during slip, as was assumed in the previous analytical solutions (21) and (22) (Carslaw and Jaeger 1959;Lachenbruch 1986;Fulton and Harris 2012). For both x/2 √ α(t − t s ) = 0 and 5 , the ratio between T thick (x, t) − T 0 and T plane (x, t) − T 0 approaches unity as the l/2 √ α(t − t s ) value decreases. No significant differences were observed among the eight curves drawn with different values of the nondimensional parameters other than l/2 √ α(t − t s ) in the parameter ranges considered here. Among the previous drilling projects, the l/2 √ α(t − t s ) value was calculated to be ≤ 5.1 × 10 −4 for the Chelungpu fault (Kano et al. 2006) and ≤ 3.2 × 10 −4 for the Wenchuan fault (Li et al. 2015) (Table 2). In the case of the plate-boundary fault in the Japan Trench (Fulton et al. 2013), the precise thickness of the primary slip zone during the 2011 Tohoku-Oki earthquake has not been reported. However, the thickness of the slip zone responsible for an individual event is generally in the range 2l ≤ 1 × 10 −1 m , as reviewed by Sibson (2003).  Ma et al. (2003Ma et al. ( , 2006, Hirono and Hamada (2010)
Next, I investigate the conditions that allow us to use the source solution (27), (28) as a good approximation of the planar solution (23). The source solution T source (x, t) with the empirical exponential slip-weakening friction can be expressed from Eqs. (27) and (28) as The results of the dimensional analysis suggest that there are four nondimensional parameters that govern the ratio between (T source (x, t) − T 0 ) and T plane (x, t) − T 0 : namely, x/2 √ αt , t s /t , c 3 t , and c 5 /c 4 . When t s /t ≪ 1 , it can be shown that the planar solution T plane (x, t) Comparison of the thick solution T thick and the planar solution T plane . The ratio between T thick (x, t) − T 0 and T plane (x, t) − T 0 and its deviation from unity are plotted against the nondimensional parameter l/2 √ α(t − t s ) when a, b x/2 √ α(t − t s ) = 0 , and c, d x/2 √ α(t − t s ) = 5 , respectively. The solutions were calculated using Eqs. (19), (20), and (23). The other nondimensional parameters were set to t s /t = 10 −8 , 10 −5 , c 3 t = 10 5 , 10 11 , and c 5 /c 4 = (0, 10) Page 11 of 13 Kaneki Progress in Earth and Planetary Science (2021) 8:67 approaches the source solution T source (x, t) (Additional file 1: Text S4). Figure 6 shows the difference between (T source (x, t) − T 0 ) and T plane (x, t) − T 0 plotted against the nondimensional parameter t s /t for various combinations of x/2 √ αt , c 3 t , and c 5 /c 4 . As in Fig. 5, the values of these nondimensional parameters were set to x/2 √ αt = (0, 5) , c 3 t = 10 5 , 10 11 , and c 5 /c 4 = (0, 10) . For both values of x/2 √ αt , the ratio between T plane (x, t) − T 0 and (T source (x, t) − T 0 ) approaches unity as the t s /t value decreases. Except for the case of (c 5 /c 4 , c 3 t) = 10, 10 5 , no significant difference was observed between the curves obtained by using different values of the dimensionless parameters other than t s /t . According to previous temperature measurements in boreholes, t s /t is ≤ 3.2 × 10 −8 for the Chelungpu fault  Table  2). In this parameter range, These results clearly indicate that the planar solution (23) can be well approximated by using the source solution (29) under the conditions of the previous temperature observations.
From the results presented in the last two paragraphs, it is clear that the temperature anomalies observed by the previous drilling projects ( Kano et al. 2006;Fulton et al. 2013;Li et al. 2015) can be well approximated Comparison of the planar solution T plane and the solution T source . The ratio between T plane (x, t) − T 0 and (T source (x, t) − T 0 ) and its deviation from unity are plotted against the nondimensional parameter t s /t when a, b x/2 √ αt = 0 , and c, d x/2 √ αt = 5 , respectively. The solutions were calculated with Eqs. (23) and (29). The other nondimensional parameters were set to c 3 t = 10 5 , 10 11 and c 5 /c 4 = (0, 10) using the source solution (29), whether the shear stress is constant or shows exponential decay with slip. This finding suggests that the slip parameters reported in these studies are appropriate.

Discussion
The results of this study clearly show that the thick, planar, and source solutions can be considered equivalent provided that the half-width of the fault zone l and the slip duration t s are negligibly small compared with twice the diffusion length 2 √ α(t − t s ) and the elapsed time after the event t , respectively. Although Fulton et al. (2013) and Li et al. (2015) have already addressed this point qualitatively, in this study, I demonstrated that these three solutions are equal, under these limiting conditions, by using the Taylor expansion to transform the solutions. Under these limiting conditions, coseismic information about the slip parameters with spatio-temporal extent is lost due to heat diffusion and thus cannot be directly inferred from the observed temperature anomalies. The most well-constrained parameter is the strength of the heat liberated in fault zone S , as already indicated by Li et al. (2015). The values and errors of the other parameters, such as shear stress and friction coefficient, therefore, need to be estimated based upon those of S , with some additional assumptions.
It is worth discussing under what circumstances the source solution (29) can no longer be regarded as a good approximation for the thick solutions (19) and (20). As a requirement for a good approximation, I here investigate the case where the approximation error is less than 1%. As shown in Fig. 6, the approximation error becomes less than 1% when t s /t = 4.0 × 10 −2 and 8.7 × 10 −3 for x/2 √ αt = 0 and 5 , respectively. If t s = 10 3 s is assumed to be the upper limit of slip duration, a value which can be considered unrealistically long, the approximation error exceeds 1% when t < 14 days . The previous drilling projects considered here were carried out at least one year after the targeted event, and the installation of thermometers within such a short time period after the event would be almost impossible. Therefore, as long as temperature data obtained from 0 ≤ x/2 √ αt < 5 are used, the source and planar solutions can be considered equivalent for all drilling projects.
As shown in Fig. 5, the approximation error between the planar and thick solutions is 1% when l/2 √ α(t − t s ) = 1.7 × 10 −1 and 2.5 × 10 −2 at x/2 √ α(t − t s ) = 0 and 5 , respectively. The thermal diffusivity of various types of rocks is typically α > 1 × 10 −7 m 2 s −1 (e.g., Vosteen and Schellschmidt 2003). Given that the previous drilling projects have been carried out at least one year after the target event, the diffusion length of a typical drilling project is expected to satisfy √ α(t − t s ) > 1.8 m . As a consequence, the fault-zone thickness 2l must be less than 0.18 m in order to keep the approximation error more than 1%. However, the thickness of a brittle fault zone is commonly less than 0.1 m (Sibson 2003), and thus the thick and planar solutions can be regarded as equivalent in most cases at 0 ≤ x/2 √ α(t − t s ) < 5 . Given the above considerations, the source solution (29) can be considered equivalent to the thick solutions (19) and (20), as long as the coseismic slip does not become significantly delocalized and the spatial range of the observation area is not very large compared with the diffusion length.

Conclusion
In this paper, I derived analytical solutions for a onedimensional heat conduction problem in and around a thick fault zone that accounts for empirical slip-weakening friction. I then verified the proposed solutions by comparing them with previously reported solutions, the solution for a planar fault problem, and numerical solutions. Using the derived solutions, I demonstrated that the temperature anomalies observed in previous fault drilling projects (Kano et al. 2006;Fulton et al. 2013;Li et al. 2015) can be well approximated by source solution (29) under the measurement conditions of those drilling projects. This result suggests that the slip parameters reported in these studies are appropriate. My results also indicate that for most drilling projects, including future ones, temperature anomalies within a thick fault zone can be well explained by using the source solution (29), instead of the thick solutions (19) and (20), as long as the coseismic slip is not significantly delocalized and the spatial extent of the temperature observations is not considerably larger compared to the diffusion length.
Additional file 1: Detailed derivations of some equations in the main text.