Geological constraints on dynamic changes of fluid pressure in seismic cycles

Fluid pressure along faults plays a significant role in fault behaviors in seismic cycles in subduction zones. When a thermal pressurization event occurs, the fluid pressure rises; conversely, when a fault-valve behavior event occurs, the fluid pressure falls. The stress state changes with seismic cycles from a reverse fault regime to a normal fault regime, as observed in both geophysical observations and geological records. Fluid pressure has been estimated for both modern accretionary prisms and exhumed accretionary complexes. However, changes in fluid pressure on seismogenic faults have not been connected to seismic cycles. Here, we quantitatively show the dynamic change in fluid pressure in a seismogenic fault with geological evidence from an exhumed accretionary complex. We found extensional veins related to seismogenic fault records that exchanged stress states the during seismic cycles. We also constrained the fluid pressure quantitatively, both at an increasing stage during an event and at a decreasing stage after an event. In this procedure, we propose new methods to constrain the magnitude of vertical stress and rock tensile strength.

Variations in fluid pressure along seismogenic faults are important in understanding earthquake mechanisms 1 . The shear strength along a seismogenic fault changes with fluid pressure during seismic cycles, and the variations in the shear strength are time-dependent 1 . Recent experimental studies have also suggested that fluid pressure can affect the frictional stability of fault zones 2,3 .
There are two opposing models for the shift in fluid pressure caused by earthquakes. One is the thermal pressurization model, in which the fluid pressure rises due to frictional heating during a seismic slip. The increase in fluid pressure reduces the shear strength and promotes slip along a seismogenic fault [4][5][6] . The other is the faultvalve model, in which the fluid pressure decreases after an earthquake due to fluid escape along newly developed fractures or faults along a subducting plate interface, enhancing the shear strength along the fault. It has been shown that variations in shear strength due to changes in fluid pressure during an earthquake cycle, as well as stress changes, may be important in controlling the onset and recurrence of rupture on seismogenic faults along subduction plate interfaces 7 . Quantitative constraints on fluid pressure variations are therefore important because these two conflicting models can coexist at different timescales during an event, potentially leading to a better understanding of seismic mechanisms and recurrence timescales.
The change in the stress state from reverse to normal fault stress regimes was observed during the Tohoku earthquake 8,9 . The stress change from the reverse to normal fault stress regime occurred at least one day after the Tohoku earthquake 8 . A geophysical study has reported that the recovery of the stress state to almost the same stress state as that of the pre-earthquake period took six years following the Tohoku earthquake 10 . Similar changes in the stress state were also identified in on-land accretionary complexes and in the cores of the Taiwan Chelung-pu Fault Drilling Project (TCDP) using a micro-fault inverse method 11,12 . However, it is difficult to geologically estimate the time scale of the changes in the stress state.
In a previous study, a poro-elastic model was used to constrain the fluid pressure decrease during seismic cycles, as demonstrated by the extensional veins of the Nobeoka thrust in southwest Japan 13 . However, because the study only focused on extensional veins formed during normal fault stress regime, it did not fully capture the change in fluid pressure during seismic cycles.
In this study, we focus on the underplating thrust zone in a mélange zone within an exhumed accretionary complex. The mélange adjacent to the thrust zone has extensional veins showing network texture, which is related to the activity of the underplating thrust zone because of their occurrences 14 . We examined the change in the stress state with seismic cycles using mixed Bingham statistics 15,16 for extensional veins. Then, the rock failure theory 13,17 and fluid pressure obtained from fluid inclusion microthermometry for the veins 18 instead of

Geological setting
The study area is an exhumed accretionary complex, the Mugi mélange in the Late Cretaceous to Early Paleogene Shimanto accretionary complex, southwest Japan 19,20 (Fig. 1a). The Mugi mélange is composed of tectonic mélanges, mainly sandstone blocks and shale matrices, with traces of basalt, chert, tuff, and red shale. Oceanic materials represent the ocean floor stratigraphy, showing a map-scale duplex structure indicating that the Mugi mélange was formed by underplating 19,20 . The Mugi mélange is divided into lower and upper sections, with maximum paleotemperatures of 130-150 °C and 170-200 °C, respectively, separated by an out-of-sequence thrust 19,21,22 (Fig. 1b). The lower section contains a sequence of three thrust imbricates (Units 1-3). The depositional age is approximately 60.5-62.3 Ma 20 . In addition, pseudotachylytes have been observed in the upper unit of the Mugi mélange, which provides strong evidence that the mélange was located within the seismogenic zone around a subduction plate interface 21,23 .
The fault zone we studied is related to underplating, comprising cataclasites and ultracataclasites mainly of basalts and a small number of cherts that steeply dip to the south or north and strike east-northeast direction 19 (Fig. 1c). The average strike and dip of the shear zone (major shear surface) is N68° E 58° S. The thickness of the cataclasites and ultracataclasites (shear zone in Fig. 1c) ranges from about 3-7 m, and the total thickness of the fault zone, including the fractured damage zone, is about 20 m (Fig. 1c). Geological evidence of thermal pressurization, such as fluidization and injection of comminuted material, dilatational brecciation in extensional jogs, and the stretching of fluid inclusions in calcite veins, has been identified in the ultracataclasites in the fault zone [24][25][26] , suggesting that the shear zone represents a fossil seismic fault. Riedel shears within the shear zone indicate a dextral sense of shear in the current setting ( Supplementary Fig. S1). The dextral sense of shears contradicts plate reconstruction showing a sinistral lateral convergence component in East Asia after the Late Cretaceous-Paleogene 27 . Therefore, the inconsistency in the sense of shear suggests that the shear zone is inverted and that there is possibly a mélange adjacent to the shear zone facing south.
The mélange immediately above the fault zone is composed of terrigenous sandstone and shale matrices. The extensional veins concentrate the boundary between the mélange and the shear zone, but only on the mélange side, and show a network texture (Fig. 1d). Because network veins grew solely along the shear zone, their occurrence must be linked to underplating after the mélange formation 14  The network veins formed after the mélange formation because they cut the mélange fabrics. Each vein in the network crosscuts with others without any specific orientation of veins, suggesting that multiple stress states were repeatedly exchanged (Fig. 1e). The microstructures of the vein minerals represent blocky and dendritic textures, suggesting that the network veins were formed by rapid precipitation from an oversaturated fluid 28,29 . Estimations of paleo-stress and driving fluid pressure. We used mixed Bingham statistics 15,16 to predict the stress state and driving fluid pressure ( P * ) for the network veins (see the Methods section). As a result, three stress states, Stress 1, 2, and 3, were identified, and P * for each stress state were estimated (  Table S1). The stress ratios (Φ) (see the "Methods" section) for each stress state were 0.1000, 0.1545, and 0.1291, respectively. Stress 1 shows a strike-slip fault stress regime with sub-horizontal σ 1 in the NS direction and sub-horizontal σ 3 in the EW direction. Stress 2 represents a normal fault stress regime with subvertical σ 1 and sub-horizontal σ 3 in the NE-SW direction. Stress 3 exhibits a normal fault stress regime with sub-vertical σ 1 and sub-horizontal σ 3 in the NNW-SSE direction. In addition, the P * values for each stress state were 0.30, 0.15, and 0.24, respectively. P * is defined as the maximum fluid overpressure normalized by the differential stress for the extension veins, and it varies between 0 and 1 13 . A P * of 0 indicates that the minimum fluid pressure during extensional vein formation is equivalent to the rock failure theory conditions 17 . Our results show P * values that are all greater than 0, ranging from 0.15 to 0.30. This suggests that the extensional veins formed with an overpressure fluid larger than that predicted by the rock failure theory 17 in the mélange just above the fault zone (see the "Methods" section).
Change in stress state with seismic cycles. The shear zone was assumed as horizontal to the decollement plane at the time the extensional veins were developed because no thermal gap at the shear zone was observed; this suggests that the shear zone was parallel to the horizontal thermal structure 30 . This is the same interpretation as that for the Minami-awa fault with pseudotachylytes in the Mugi mélange 21 . Therefore, the shear plane must be rotated horizontally to estimate the stress state in the event. Aaccrodingly, the estimated principal stress orientations were rotated by 58° to the south with a rotation axis of N68° E based on the average strike and dip of the fault zone (N68° E 58° S  (Fig. 2). Stress 1 shows a normal fault stress regime with sub-vertical σ 1 and a WSW-ENE sub-horizontal σ 3 . Stress 2 indicates a non-Andersonian stress regime with principal stresses, NW, intermediately plunged σ 1 , NE, intermediately plunged σ 2 and SW sub-horizontal σ 3 . Stress 3 is a reverse fault stress regime with NNW-SSE sub-horizontal σ 1 and sub-vertical σ 3 . Although Stress 2 shows a non-Ander-   Firstly, we used a -�σ diagram ( is the fluid pressure ratio, = P f /σ v ) to constrain the fluid pressure ratio for extensional veins 7,31,32 . In the diagram, negative and positive �σ values represent the normal and reverse stress regimes, respectively (Fig. 3). After Eqs. (6), (8), (13), and (15) are divided by the vertical stress, the relationships between the maximum or minimum fluid pressure ratios and differential stress for the extensional veins in reverse and normal fault stress regimes can be drawn in the -�σ diagram (Fig. 3). The lines from Eqs. (6) and (8) indicate the lower limit of the fluid pressure ratio, whereas the lines in Eqs. (13) and (15), including P * , represent the upper bound of the fluid pressure ratio (Fig. 3). However, the vertical stress and rock tensile strength are still not constrained.
In order to constrain these parameters, we used previous estimates of fluid pressure based on fluid inclusion analysis from the network veins, which estimated that the fluid pressure was between 107(+ 8/− 4.1)-149.4 MPa 18 . As the previous study did not classify the extensional veins based on stress regime, the range of fluid pressures from the fluid inclusion was expected to be included in different stress regimes because the extensional veins were collected randomly in the previous study. When the differential stress is the maximum value of 4T s (see the "Methods" section), the fluid pressure of 149.4 MPa must be the maximum in the reverse fault stress regime at point A in Fig. 3. However, the minimum fluid pressure of 107(+ 8/− 4.1) MPa was in the normal fault stress regime at point B, as shown in Fig. 3. When T s varies from 1 to 20 as reasonable values 33,34 , and using Eqs. (7) and (16), the relationships between T s and σ v at points A and B are constrained by the maximum differential stresses, �σ = 4T s (Fig. 4). Vertical stresses are calculated using σ v = ρgz (where ρ is the rock density, g is the gravitational acceleration, and z represents depth), meaning that σ v can be converted to a depth ( ρ is 2.600 kg/ m 3 and g is 9.8 m/s 2 in the conversion). The cross-cutting relationship in the network veins suggests that the extensional veins formed at similar depths and rock tensile strengths in the normal and reverse fault stress regimes. Therefore, the cross point of the two slopes in the relationships between T s and depth for the normal and reverse fault regimes indicates the condition for extensional vein formation (Fig. 4). Consequently, T s and depth were constrained to 6.94-9.38 MPa and 5.14-5.33 km, respectively. In addition, the geothermal gradient was estimated to be approximately 24.39-29.17 °C/km based on the maximum paleo-temperature in the lower section of Mugi Mélange, which ranges from 130 to 150 °C 19 . The constrained values for tensile strength, depth and geothermal gradient are reasonable, which supports our expectation that extensional veins include both reverse and normal fault stress regimes.
Quantifications of dynamic change in fluid pressure along a seismogenic fault. The difference between the lower and upper limits for fluid pressure in the -�σ diagram is �σ P * , and the maximum value of  Figure 3. Relationship between the fluid pressure ratio ( ) and differential stress ( �σ ). Differential stress equals the differences between horizontal stress ( σ h ) and vertical stress ( σ v ). Positive and negative differential stresses correspond to reverse and normal fault stress regimes. Lines from Eqs. (6), (8), (13), and (15)  www.nature.com/scientificreports/ the difference is 4T s P * , as shown by the differences between Eqs. (6) and (13) for the normal fault stress regime as well as Eqs. (8) and (15) for the reverse fault stress regime. The maximum differences are quantified to be 6.7-9.0 MPa for the normal fault stress regime and 8.3-11.3 MPa for the reverse fault stress regime when we use the constrained T s described above.
The absolute value of the fluid pressure ratio was also estimated using the constrained σ v . The fluid pressure ratio in the reverse fault stress regime was approximately 1.12 at the maximum point A (Fig. 3), which exceeded the lithostatic pressure ( > 1 ). This is a requirement for an extensional vein in the reverse fault stress regime, as shown in Eqs. (8) and (15). Because the condition with > 1 is unstable, we can interpret that a large fluid pressure was achieved in dynamic events. Previous studies have shown that the studied fault zones include geological evidence of thermal pressurization, which promoted seismic slip in this fault zone [24][25][26] . In our study, the dynamic fluid pressure reaching an over-lithostatic pressure indicated thermal pressurization during a seismic event. Our quantitative constraints indicate that fluid overpressures of about 6.7-11.3 MPa above the minimum requirement for the formation of extensional veins in a reverse fault stress regime were achieved at the time of the seismic event. The unstable condition with > 1 also suggests that the fluid pressure increased by thermal pressurization, overcoming the dilatancy effect which tends to decrease the fluid pressure [35][36][37] .
However, our findings show that the fluid pressures in the reverse and normal fault stress regimes differ. The fluid pressure in the reverse fault stress regime was higher than that in the normal fault stress regime (Fig. 3). The smaller fluid pressure in the normal fault stress regime may indicate that fluid pressure decreases related to a fault-valve model, because the change in stress state between normal and reverse fault stress regimes can be related to seismic cycles, as discussed above. A decrease in fluid pressure may stabilize faults in a normal fault stress regime [35][36][37] . Our results also quantified the fluid pressure drop from the maximum to minimum fluid pressure after a seismic event, which was approximately 42 MPa from point A to point B in Fig. 3.
In our study, we constrained fluid pressure variations only in the extensional veins along a seismic fault, but we cannot rule out the possibility that the fluid pressure can decrease further after an earthquake to a level where extensional veins cannot form but instead promote the formation of shear veins 7,31 when the differential stress becomes larger than 4T s ; this means that the fluid pressure decrease can be larger than 42 MPa. It must be noted that the constrained fluid pressure range in this study applies only to the extensional veins, as was also true in a previous study 13 . Shear veins were observed close to the shear zone within the mélange in this study area, although the relationship between the shear veins in the mélange and the shear zone below the mélange needs to be examined carefully.
In summary we proposed new methods to constrain the depth, rock tensile strength, and dynamic fluid pressure variations from geological data. Our study successfully constrained the dynamic change in fluid pressure quantitatively in a seismic cycle, both during the increase stage caused by thermal pressurization and during the decrease stage caused by fault-valve behavior in a seismogenic fault zone. To fully evaluate fluid pressure fluctuations in seismogenic fault zones, investigations on other modes of failure, such as shear veins and hybrid extensional shear veins, must be expanded 7,31,32 .   (Supplementary Table S1). We can approximate stress ratio (Φ) by the ratio of the concentration parameters, � = κ 2 /κ 1 38 . The stress ratio is defined as follows: where σ 1 , σ 2 , and σ 3 are the maximum, intermediate, and minimum principal stresses, respectively. The number of stress states (K) was set by the user as 1-5. The minimal value of the Bayesian information criterion (BIC) for each result with variables K was used to determine the optimal number of stress states in the mixed dataset 15,16 . Consequently, the optimal number of stress states in our dataset was three ( Supplementary Fig. S2). The method also classifies the extensional veins that fit each stress state. The fitness of each extensional vein to a stress state is expressed by "membership" 15 . The membership indicates the probability that each extensional vein belongs to a stress state, ranging from 0 to 1. For each stress, the normal stresses normalized by the differential stress were calculated on each extensional vein. With the likelihood, the membership decreases as normalized normal stress increases. A previous study 39 proposed the 95th percentile point of the probability, termed the driving pressure index (DPI), for the driving pressure ratio ( P * ) under fluctuated pressure conditions. DPI is defined as, where κ 1 and σ are the concentration parameter and normalized normal stress.
Solving Eq. (2), we obtain We can use DPI as a representative value of P * for the Bingham distribution model 39 . P * is defined as 13 , where P f and T s are fluid pressure and rock tensile strength. As three stress states were estimated from our dataset in this study, three P * were obtained for each stress state. However, in this study, we focused on two stress states for the normal and reverse fault stress regimes. The P * contains the lower bound of the maximum fluid pressure for the extensional vein formation 38 .
Constraints on fluid pressure using rock failure theory. Brittle-rock failure is controlled by coefficients of internal friction, differential stress, and rock tensile strength ( T s ) 17,40,41 . Extensional cracks formed under a low differential stress of �σ < 4T s 17 . Extensional veins are formed when the fluid pressure exceeds the sum of the normal stress on the extensional vein (minimum principal stress) and T s 17,42 . Hence, extensional veins were formed when In the case of the normal fault stress regime, the maximum principal stress is equal to the vertical stress. Therefore, the minimum fluid pressure and vertical stress during the extensional vein formation under a normal fault stress regime can be approximated as follows: where σ v and �σ are vertical stress and differential stress. Therefore, Because the minimum principal stress in a reverse fault stress regime is equal to the vertical stress, the fluid pressure during extensional vein formation can be expressed as The fluid overpressure ( P O ) is the excess fluid pressure greater than the minimum fluid pressure for the extensional veins, which is the sum of the minimum principal stress and rock tensile strength from Eq. (5). The fluid overpressure 13 is defined as: (4) P * = P f − σ 3 − T s σ 1 − σ 3 , (5) P f > σ 3 + T s .
(8) P f _reverse = σ v + T s , (9) σ v = P f _reverse − T s . www.nature.com/scientificreports/ The driving fluid pressure ratio ( P * ) is defined in Eq. (4), and it is rewritten as the P O normalized by �σ 13 as follows: We can thus express P f arranging Eq. (11) as follows In the case of a normal fault stress regime, the minimum principal stress can be substituted with σ v − �σ . Therefore, in the normal fault stress regime, the fluid pressure and vertical stress can be calculated using Eq. (12), as follows: In a reverse fault stress regime, because the minimum principal stress is equal to the vertical stress, the fluid pressure during extensional vein formation can be expressed using Eq. (12) as

Data availability
The data supporting the findings of this study are available in Supplementary information.  (12) P f = σ 3 + T s + �σ P * .