Silica precipitation potentially controls earthquake recurrence in seismogenic zones

Silica precipitation is assumed to play a significant role in post-earthquake recovery of the mechanical and hydrological properties of seismogenic zones. However, the relationship between the widespread quartz veins around seismogenic zones and earthquake recurrence is poorly understood. Here we propose a novel model of quartz vein formation associated with fluid advection from host rocks and silica precipitation in a crack, in order to quantify the timescale of crack sealing. When applied to sets of extensional quartz veins around the Nobeoka Thrust of SW Japan, an ancient seismogenic splay fault, our model indicates that a fluid pressure drop of 10–25 MPa facilitates the formation of typical extensional quartz veins over a period of 6.6 × 100–5.6 × 101 years, and that 89%–100% of porosity is recovered within ~3 × 102 years. The former and latter sealing timescales correspond to the extensional stress period (~3 × 101 years) and the recurrence interval of megaearthquakes in the Nankai Trough (~3 × 102 years), respectively. We therefore suggest that silica precipitation in the accretionary wedge controls the recurrence interval of large earthquakes in subduction zones.


Calculation of sealing time using the model of crack in reservoir via advection and kinetic reactions
The sealing time of an extension crack by quartz precipitation on pre-existing quartz surfaces is estimated based on the kinetic equation for precipitation as follows 1 : where tr is the residence time of fluid in the crack. The SiO2 concentration in the input fluid (CSiO2) in mg/kg(H2O) is the quartz solubility under the initial P-T conditions of the host rock. The value of CSiO2,Qtz,eq is the quartz solubility at fluid pressure in a crack (Pcrack) after a fluid pressure drop (ΔP) as follows: where Phost is fluid pressure in the host rock before a fluid pressure drop. The precipitation rate constant of quartz (k-) was determined empirically as follows 2 : where T is temperature in K. The model of an extension crack assumes a disk shape. The area of one crack wall (Acrack) in m 2 and the volume of the crack (Vcrack) in m 3 are as follows: where lv and wv are the length and width of a quartz vein as measured in field studies (e.g., Figs. 2a,b), which are assumed to be the diameter and height of a disk-shaped crack in our model, respectively. The area of quartz precipitation and the growth of quartz on crack walls (AQtz) in m 2 is computed as follows: where Qtz is the proportion of quartz in the host rock. The mass of fluid in the crack (Mcrack) in kg is determined as follows: where Vsp is the specific volume of water in m 3 /kg at 250 °C and at the pressure achieved after the fluid pressure drop in m 3 /kg. At the beginning of the calculation, the apparent flow rate (Qap) in m 3 /s is set to estimate the apparent residence time of fluid in the crack (tr,ap) in s as follows: Replacement of the residence time (tr) in equation (2) with tr,ap results in a change in SiO2 concentration due to quartz precipitation in the crack (ΔCSiO2) in mg/kg(H2O) under the conditions of the apparent flow rate. The total volume of quartz precipitate in 1 kg of fluid (VQtz) in m 3 /kg is where  Qtz is the density of quartz (2.65 g/cm 3 ). Thus, the sealing time of the whole crack volume (ts) in s and the total volume of fluid required to seal the whole crack (Vf,seal) in m 3 , or equivalently, for the volume of silica precipitated to equal that of the crack, are estimated as follows: The volume of a spherical reservoir around the crack (Vres) in m 3 required to seal the crack volume (Vf,seal) in m 3 can be written as the following pair of functions: where  is the porosity of the host rock and R is the radius of the reservoir in m. As a result, R is determined using both equations (13) and (14). Thus, the flow rate into a crack (QD) in m 3 /s can be determined from Darcy's law as follows: where  is the permeability of the host rock (sandstone),  is the viscosity of the fluid, and ΔP in Pa is the pressure difference between the host rock and the crack, or equivalently, a fluid pressure drop in equation (3). Here, the apparent flow rate (Qap) used to determine R in equations (2)(14) is not the same as the flow rate into a crack (QD) in Darcy's law in equation (15). Thus, we calculate the optimal flow rate (Qk) in m 3 /s in the loop calculation. First, Qap is set (small value) to calculate R and QD by using equations (2)(15). Qap is compared with QD, and if Qap is lower than 99.9% of QD, Qap is increased in set step (e.g., Qap,i = 0.1 m 3 /s, step = 0.01 m 3 /s, Qap,i+1 = 0.11 m 3 /s). During loop calculation, when Qap = (or slightly higher than) QD  0.999, we set Qk = Qap = QD, which results the optimal volume of fluid, crack sealing time and other parameters in the conditions.

Realistic conditions around the Nobeoka Thrust
Temperature and pressure in the host rock at 10 km depth are set to 250 °C and 260 MPa 3,4 , respectively. The host rock of sandstone 3 is composed of 52% (Qtz = 0.52) quartz 5 . The neutron porosity of the host rock () falls in the range 3% to >10%, as recorded by geophysical wireline logs across the Nobeoka Thrust in a borehole experiment in 2011 by SRED and Raax Co., Ltd. 6 ; we used 3% ( = 0.03) in this work because fewer cracks existed before the earthquake. The permeability of the host rock () falls within the range 1 × 10 −20 to 1 × 10 −18 m 2 as measured by a triaxial pressure apparatus 7 ; we used an average value,  = 1.0  10 −19 m 2 , in this study. The rate constant for silica precipitation via quartz growth on quartz surfaces (k-) is 7.4  10 −6 s −1 at 250 °C, following the Arrhenius relation of the experiment 2 . The SiO2 concentration in the fluid in the host rock (CSiO2) is 6.0  10 2 mg/kg(H2O) as a result of pure water saturated via quartz at 250 °C and 260 MPa (Fig. 3) 8 . The pore fluid pressure in a crack is assumed to decrease close to the maximum hydrostatic pressure. The hydrostatic pressure at 10 km depth is 98 MPa. Therefore, the difference in pressure (ΔP) is ~160 MPa (Fig. 3). At 250 °C, silica precipitation occurs only in the liquid phase because boiling occurs at ≤4 MPa. The specific volume and viscosity of water is Vsp = 1.1  10 0 1.2  10 0 cm 3 /g 8 and  = 1.5  10 -4 1.1  10 -4 Pa 9 , respectively, depending on ΔP = 0.01160 MPa. When pore fluid pressure decreases to hydrostatic (ΔP = 160 MPa), the solubility of quartz in pure water decreases to 5.0  10 2 mg/kg(H2O) (Fig. 3) 8 . Under these conditions, the quartz solubility in seawater (3 wt% NaCl) is 5.5  10 2 mg/kg(H2O) and 4.8  10 2 mg/kg(H2O) at lithostatic and near-hydrostatic pressure (ΔP = 160 MPa), respectively 8 . The maximum difference in quartz solubility in pure water, 9.4  10 1 mg/kg(H2O), is larger than that in sea water, 7.8  10 1 mg/kg(H2O). The quartz dissolution rate is strongly dependent on pH rather than salinity 10 . Fluid inclusion analysis in the Shimanto Belt of Kyushu, Japan revealed that typical fluid salinity is up to 5.5 wt% NaCl and is mainly similar to or lower than that of seawater 11 . The stability fields of the vein minerals indicate that extension veins formed from relatively oxidized, locally derived pore fluids of neutral pH 12 , in the case of a dissolution rate in sea water of about one order higher than that in pure water 10 . Therefore, the sealing time of a crack would be similar to or shorter than that estimated in this study (Figs. 4-6 in the main manuscript).

Diffusion model
In the diffusion model of ref. 13, the sealing time of crack by diffusion (tD) of a crack is written as follows: where Vq is the molar volume of quartz (2.2  10 −5 m 3 /mol),  is matrix porosity (porosity of the host rock, 3%),  is tortuosity (1), Df is the pore fluid diffusion coefficient (1  10 −8 m 2 /s), and Fd is the volume fraction of quartz dissolved out of the depletion zone (0.40; i.e., 40%). The length of the quartz fibre (l) in m is the same as the mode aperture width of the quartz vein (wv = 5.2  10 1 μm) in this study. The local equilibrium concentrations at the sites of dissolution and growth (C1, C2) in mol/m 3 are the same as the concentrations of SiO2 in the fluid at the initial lithostatic pressure (CSiO2) and after fluid pressure drop (CSiO2,Qtz,eq) at 250 °C, respectively. Since the fibrous overgrowths by fluid diffusion are limited by dissolution rather than precipitation 10 , the linear dissolution rate constant (kl) in m/s in equation (16), is determined by using molar volume of quartz (Vq), equilibrium constant and precipitation rate constant (k-) 2 , based on the quartz-water reaction 1 . The dissolution rate constant (kl) is in range from 1.6  10 −12 m/s to 1.1  10 −12 m/s at ΔP from 0.1 MPa to 160 MPa, respectively.

Size variations in quartz veins and changes in the relative porosity index
The presence probability (or frequency) of extensional quartz veins along the Nobeoka Thrust (n,m) is evaluated as follows: where xn and ym are the nth and mth relative frequency values of aperture width (wv,n) and length (lv,m), which are determined at each log-0.1 step in the ranges 1  10 1 -1  10 3 μm and 1  10 0 -1  10 2 cm in the histograms, respectively (Fig. 2a,b). In these classifications, the values of both n and m are 1-21. The total number of possible n-m combinations is 441 (21  21), 225 of which are considered presence-probable conditions (n,m > 0). To estimate the porosity change, we determined the relative porosity index (t) from 0 to 1 at time tp in year. t = 0 is just before an earthquake, where the porosity is the initial value of the host rock without any cracks; t = 1 is just after earthquake (tp = 0), when all cracks of presence-probable size have opened ( Fig. 6). First, the volume of the disk-shaped crack (n, m) (Vcrack,n,m) and its sealing time (ts,n,m) are calculated in the model of this study. Next, the contribution ratio of the crack (n, m) to the total porosity (n,m) is weighted and normalized using the presence probability of the crack size as follows: n,m = (Vcrack,n,m  n,m)/∑∑(Vcrack,n,m  n,m).
Based on the assumption that a crack sealed in less time contributes to porosity healing occurring earlier, we sorted the crack sealing times (n, m) (ts,n,m = ts,i), where i is the sort index ranging from 1 to 225, corresponding to shorter to longer times, respectively. The relative contribution ratio of crack size to porosity was also sorted by sealing time, as n,m = i. Here, we assume that the precipitation rate is constant until a crack is sealed completely. As a result, the contribution ratio of the crack (n, m), rewritten as the crack sort index i on relative porosity i at time tp, can be determined as follows: i = i,  (t/ts,i), when tp < ts,i.
The sealed crack index at time tp is the maximum value of i that satisfies equation (19), and the total change in relative porosity index at time tp is estimated as follows: where 0 = 1 (tp = 0) when all cracks of presence-probable size open. In addition, the cumulative number index of sealed crack (t = 0-1) is a rate between a cumulative number of cracks completely sealed (i.e., satisfy the conditions of equation (19)) and total number of possible length-aperture width combinations (225).