Pore Pressure Regime and Fluid Flow Processes in the Shallow Nankai Trough Subduction Zone Based on Experimental and Modeling Results from IODP Site C0023

Pore fluid pressure is a critical parameter controlling the dynamics of subduction zones and related geological hazards, but quantitative constraints on pore pressure are sparse. Here we provide robust estimates of pore pressure in the shallow Nankai Trough subduction zone based on consolidation experiments and numerical models. Our experimental data reveal excess pore pressures with nearly constant ratios of pore pressure to lithostatic stress for the accreted and underthrust sediments, implying the décollement is not a fluid barrier as previously inferred. Also, our research suggests repeated lithostatic pore pressure transients in the last 3 kyrs, probably associated with the propagation of the frontal thrust. The coupled pore pressure dissipation and chemical transport model results let us conclude the updip migration fluid originates from an overpressured horizon nearby, confirming the view of a heterogeneous plate boundary with overpressured and permeable patches along the fault surface that shift in time and space. The resulting high excess pore pressure leads to the shear strength of Nankai subduction thrust as low as 0.4–0.8 MPa. This provides quantitative constraints on the weakness of subduction thrusts at the toe of accretionary prisms and may explain why seismic slip can propagate all the way to the trench in some earthquakes (e.g., 2011 Tohoku event).

. Underthrust sediments show higher porosity than expected, which has been interpreted as an indicator of elevated pore pressure inherited from rapid vertical loading by the overriding wedge as well as inhibited drainage across the subduction thrust (e.g., Gamage & Screaton, 2006; G. F. . In accreted sediments better drainage and additional tectonic compaction are supposed to cause lower porosities in comparison to reference sites outside the subduction system. The change of stress field compared to the references seaward of the accretionary prism, however, complicates pore pressure estimation. Most existing pore pressure studies circumnavigate this problem by calculating fluid sources from an assumed porosity field for the accreted sediments that in turn drive pore pressure in fluid flow models (e.g., Bekins & Dreiss, 1992;Saffer & Bekins, 1998;Screaton et al., 1990). In a recent study, Flemings and Saffer (2018) proposed a generalized porosity-based approach to predict pore pressure in accreted sediments by combining a critical state soil model and the knowledge of stress path of sediments entering an accretionary prism. This method takes the porosity loss driven by shear stress into consideration and provides a more constitutive relationship between porosity and stress. All these previous approaches have been used with success to understand the overall time-averaged pore pressure pattern in subduction zones. But they provide only limited insight into detailed pore pressure transients because of assumed priori stress states and porosity fields. The need for a better understanding of transient pore pressure becomes apparent by the discrepancy between predicted near lithostatic pore pressure and moderate excess pore pressure measured in sealed ocean drilling boreholes along the shallow plate interface of subduction zones (also called décollement). For example, pore pressure in Barbados accretionary complex is only one-third of the lithostatic stress (K. Becker et al., 1997).
Here, we focus on the Nankai Trough subduction zone where five M w ≈ 8 earthquakes in the last four centuries caused tsunamis (Ando, 1975). Numerous studies have investigated excess pore pressure in the shallow portion, but were either based on geophysical data (Bangs et al., 2004;Tsuji et al., 2008), time-averaged fluid flow models (e.g., Saffer & Bekins, 1998) or hydrogeological models focusing on the underthrust sequence (Screaton et al., 2002;Skarbek & Saffer, 2009). In a recent study, Park et al. (2014) showed a large portion of the shallow décollement in the Nankai Trough is characterized by impedance-decreases ( Figure 1a). They explained the decreasing acoustic impedance with (1) a high-impedance accretionary prism versus a low-impedance underthrust section or (2) increase of pore pressure within the fault zone, which causes decreases of seismic velocity and thus reversal of acoustic impedance across the fault, as discovered in the Oregon accretionary prism (J. C. Moore et al., 1995;Tobin et al., 1994).
To test if the shallow plate boundary in the Nankai Trough is overpressured, we collected samples at the top of the major reverse fault zone and conducted a multi-methodological approach, which takes into account the stress state of the incoming sediment, its stress path during burial and the timing of pore pressure build up. In addition to samples close to the décollement, we also sampled accreted and underthrust sediments to characterize the pore pressure regime in the whole sediment succession. For the pore pressure estimates, we conducted uniaxial consolidation experiments, which is a widely and successfully utilized method to derive in situ pore pressures (e.g., Saffer, 2003Saffer, , 2007. The analysis of the data accounts for the compressive strengthening from secondary consolidation (creep) and diagenesis observed in the incoming sediments and the change of stress field in the prism sediments. We compared the inferred pore pressure estimates with those predicted from the recently developed critical state soil mechanics model by Flemings and Saffer (2018). Through a one-dimensional coupled pore pressure dissipation and chemical transport model, we also addressed the transient nature of pore pressure build up along the shallow décollement and constrained the timing of the fault activities.

Geological Setting
The central Nankai Trough subduction zone forearc, where the Philippine Sea plate is subducting to the northwest beneath the Eurasian plate (Figure 1a), has been recently sampled at Site C0023 by the International Ocean Discovery Program (IODP) Expedition 370. The drill site is located on the frontal accretionary prism off Cape Muroto, ∼2 km landward of the trench, and was established to investigate the limits of life in the deep subseafloor biosphere (Heuer et al., 2017). The sediment succession of Site C0023 is divided into 4 main lithostratigraphic units and consists of sediments offscraped from the Philippine Sea plate (Heuer et al., 2017). The lithostratigraphy is similar to that of nearby Ocean Drilling Program (ODP) Sites 1174 and 808, with a ∼38 m thick décollement (Figure 2a) located between 758 and 796 mbsf (meters below seafloor) in bioturbated hemipelagic mudstones of the Lower Shikoku Basin facies (Heuer et al., 2017).
Shipboard measured porosity and pore water geochemistry at Site C0023 (Figures 2b and 2c) provide qualitative evidence for excess pore pressure and active fluid flow. According to the prevailing model of the shallow subduction zone hydrogeology, the porosity offset to higher values across the décollement is supposed to reflect high excess pore pressure in low permeable clay-rich underthrust sequence from rapid loading and inhibited fluid expulsion across the subduction thrust (e.g., Bray & Karig, 1988;Gamage & Screaton, 2006;G. F. Moore et al., 2001). The general chloride profile is similar to that of ODP Sites 1174 and 808 with a broad low-chloride anomaly, reflecting fluid freshening by in situ clay dehydration and possibly an additional contribution of freshened fluids by updip migration (Saffer & Bekins, 1998;Saffer & McKiernan, 2009 Heuer et al. (2017). Solid dots mark ODP Sites 808, 1173 and 1174, as well as IODP Sites C0011 and C0023. Thick blue bars mark décollements with reverse polarity reflections and red bars mark décollements with normal polarity reflections after Park et al. (2014). Insets: The location of seismic line 332 and tectonic background of Nankai Trough. (b) Pre-stack depth migrated 2D seismic section of Inline 332 in the area of Site C0023 after Heuer et al. (2017). Solid reverse triangles mark the position of Site C0023 and projected positions of Sites 808 and 1174. The dashed line marks the proto thrust zone. IODP, International Ocean Discovery Program; ODP, Ocean Drilling Program. migration of fluids formed by clay mineral dehydration at greater depth (Saffer & Bekins, 1998). In situ clay dehydration as an alternative explanation can be excluded because no significant increase in clay content is observed in the horizon (Heuer et al., 2017).

Theory of Consolidation Experiment and Sample Materials
Consolidation is the gradual volume reduction of fully saturated sediments with increasing effective stress (Knappett & Craig, 2012;Terzaghi, 1943). The consolidation characteristics of sediments can be assessed through an uniaxial consolidation experiment, for example the incremental loading oedometer test, the constant rate of strain loading test, and the one-dimensional consolidation test in a triaxial apparatus. For the incremental loading oedometer test, the specimen is held inside a confining ring and lies between two porous stones (Figure 3a). A loading cap sits on the upper porous stone and passes the applied load onto the specimen. The whole assembly is placed into an oedometer cell, which is filled with water and allows the water inside the specimen to drain out. With the increase of the applied load, the specimen is compressed, leading to the upper porous stone and loading cap moving down. The resulting displacement is recorded by a linear variable displacement transducer (LVDT). Because the confining ring prevents the specimen from lateral deformation, the strain is uniaxial and can be calculated by the vertical displacement. On the semi-logarithmic stress-strain plot, a yield point which separates the elastic from elasto-plastic deformation behavior can be identified (Figure 3b). This yield point is commonly regarded as the maximum vertical effective stress that the sediment has sustained and is termed preconsolidation stress pc   (e.g., Casagrande, 1936). The sediment is said to be normally consolidated when the in situ vertical effective stress ( v   ) equals pc   while it is called overconsolidated if pc   is larger than v   (e.g., Knappett & Craig, 2012).
For sediments that have only experienced uniaxial consolidation (K 0 consolidation) and no unloading ever ZHANG ET AL.  happened, pc   is commonly interpreted as v   (Saffer, 2003). Hence, assuming pc   is equal to v   , the pore pressure P can be computed based on Terzaghi's effective stress theory: where v  is the lithostatic stress calculated by shipboard measured bulk density data (e.g., Dugan & Germaine, 2008).
For this study, a total of 11 specimens sampled from 450 to 1022 mbsf at Site C0023 (Table S1) were tested by incremental loading consolidation experiments following ASTM D2435/D2435M-11 (ASTM International, 2011). In order to maintain the natural water content, all samples were sealed and stored in the core liner at 4°C before testing. The specimens were trimmed into cylinders with a diameter of 25.4 mm and a height of ∼10 mm and were kept saturated in a fixed ring oedometer using an artificial seawater brine (3.5%). The inner surface of the confining ring was highly polished and coated with silicone grease to minimize the effects of friction between the specimen and confining ring. In order to avoid intrusion of sediments into the pores of the porous stone, filter screens were placed between the specimen and the porous stones. During our consolidation experiments, we first doubled the load and then used small increments in the vicinity of pc   so that the transition between the recompression and virgin compression lines could be determined with more precision (e.g., Brumund et al., 1976). Although pc   was unknown before the experiments, we utilized hydrostatic vertical effective stress ( vh   ) which was calculated by subtracting hydrostatic pore pressure from v  as a reference. For each load, a constant load increment duration of 24 h was performed to allow complete dissipation of excess pore pressure. Vertical load and sample height were continuously recorded throughout the experiments.
In this study, we determined pc   using a modification of the widely used Casagrande method (Casagrande, 1936;McNulty et al., 1978). However, due to the scale dependence of pc   by graphical construction (Clementino, 2005), we also used the Pacheco Silva method (1970) to estimate pc   for comparison since it requires only modest subjective interpretation (Grozic et al., 2005). Moreover, we determined the confident range of pc   by estimating the maximum and minimum possible pc   values. The maximum pc   was acquired from effective stress on the stress-strain plot at which the strain begins to deviate from the virgin compression line (e.g., Saffer, 2003). Similarly, the minimum pc   can be given by the first point deviating from the recompression line. However, this method produces over-conservative estimates of minimum pc   because the starting stress of transition from recompression to virgin compression is significantly lower than pc   . For example, Jacobsen (1992) proposed pc   = 2.5 k   , where k   is the point with maximum curvature in Casagrande method and also larger than the starting stress of transition from recompression to virgin compression. Therefore, we plotted our consolidation data with the work-stress method (D. E. Becker et al., 1987)   (e.g., Dugan & Germaine, 2008;Santagata & Germaine, 2002). The estimating procedures for the modified Casagrande method and all laboratory data are shown in the supporting information (Text S1, Table S2, and Figures S1 and S2) and available in Zhang (2020).

Correction for Secondary Consolidation and Diagenesis
Although Equation 1 was successfully used to infer pore pressure (e.g., Dugan & Germaine, 2008;Saffer, 2003), it has been recognized that secondary consolidation and diagenesis could shift pc   to higher values (e.g., Bjerrum, 1973;Morgan & Ask, 2004). In this study, we accounted for these effects and corrected pc   . The quantitative relation between pc   and v   is given by the overconsolidation ratio (OCR): For normally consolidated sediments OCR equals unity, whereas for sediments subjected to unloading, through erosion for example, OCR is larger than unity (e.g., Knappett & Craig, 2012). OCR > 1 is also possible when secondary consolidation (creep), cementation or thixotropy causes a strengthening of the sediment (e.g., Bjerrum, 1973;Locat & Lefebvre, 1986). Strictly speaking, OCR can never be less than one (e.g., Knappett & Craig, 2012). However, many authors replaced v   by vh   in Equation 2 to calculate OCR owing to difficulty in obtaining v   and presented OCR<1 (e.g., Kitajima & Saffer, 2012;Saffer, 2003). Here we used the strict definition of OCR (Equation 2).
Previous consolidation tests on sediments from reference Sites 1173 and C0011 and the underthrust section of Sites 808 and 1174 (Morgan & Ask, 2004;Bellew, 2004;Hüpers & Kopf, 2012;Guo & Underwood, 2014;Kitajima & Saffer, 2012) have revealed OCR values mainly ranging from 1.2 to 2.3 in the Nankai Trough ( Figure 4 and Table S3). The v   values at reference Sites 1173 and C0011 were inferred from vh   since low sedimentation rates prevented pore pressure in excess of hydrostatic (e.g., Kitajima & Saffer, 2012;Screaton et al., 2002), whereas at Sites 808 and 1174, v   (2004), Bellew (2004), Hüpers and Kopf (2012), Guo and Underwood (2014) and Kitajima and Saffer (2012). OCR (overconsolidation ratio) = pc Dashed line marks the depth where an obvious increase in the measured percentages of illite in illite/smectite mixed-layer clays is observed (Steurer & Underwood, 2003), which may strengthen the sediment. Note that the preconsolidation stress generally increases linearly with depth except for the two data points from Site 808 which are obviously underestimated (Morgan and Ask, 2004). secondary consolidation and diagenesis have been invoked to explain the OCRs >1 (Hüpers & Kopf, 2009;Morgan & Ask, 2004). Secondary consolidation results in continuous settlement under constant effective stress, which enhances the stiffness of sediments (Karig & Ask, 2003) and shifts the yield point in consolidation test to higher effective stress (Figure 5). For sediments subject to diagenesis, for example, cementation, the intergranular bonding strengthens the sediment matrix and further increases pc   (Morgan & Ask, 2004). In contrast, overconsolidation driven by lateral tectonic force is negligible because the underthrust section in subduction zones and reference sites is widely considered to deform in uniaxial compression (e.g., Screaton et al., 2002;Skarbek & Saffer, 2009). To investigate the mechanism for the strength evolution with increasing burial, Perret et al. (1995) compiled a set of parameters with idealized subseafloor profiles. In addition to OCR, they used v   and pc   profiles as well as the overconsolidation difference (OCD) which is defined as Olsen et al., 1986). According to those authors, secondary consolidation leads to a constant OCR with depth, while pc   and OCD increase linearly versus depth with gradients higher than those for normally consolidated sediment. We note that the Nankai inputs generally have similar character to this pattern with an average OCR of 1.76, although the data show a considerable scatter. The latter can be explained by the formation and breakdown of structural links built by cementation, which have been observed by previous studies (e.g., Hüpers et al., 2015;Sample et al., 2017;Steurer & Underwood, 2003). The slight increase of OCR with depth may also be related to the progressive diagenesis at greater depth (Morgan & Ask, 2004).
Our assessment of the available preconsolidation stresses for the Nankai inputs shows that the assumption of normally consolidated sediments does not apply ( Figure 4). Secondary consolidation and diagenesis have shifted the yield point in consolidation tests to higher stresses, which has to be addressed in the calculation of pore pressure from pc   . We therefore modified Equation 1 to account for the increase in pc   by employing Equation 2 to calculate v   : Because a solid relationship between OCR and depth is absent, we divided the Muroto transect into two sections. The boundary between the two sections was set at 800 mbsf where an obvious increase is observed in the measured percentages of illite in illite/smectite mixed-layer clays (Steurer & Underwood, 2003). The growth of authigenic illite (e.g., Morgan & Ask, 2004) and possible precipitation of quartz (e.g., Thyberg & Jahren, 2011) during enhanced smectite-to-illite transformation are expected to strengthen the sediment. Therefore an average OCR of 1.54 from reference Site 1173 at similar depth (Bellew, 2004;Morgan & Ask, 2004) was used for the upper section, whereas for the lower section an higher average OCR of 2.14 was utilized based on subducted samples from nearby Sites 1174 and 808 (Morgan & Ask, 2004). To understand the effects of OCR values on pore pressure estimates, we calculated the error bar by employing a minimum and maximum OCR value for both upper and lower sections, respectively. Due to the complexity of OCR values, the minimum and maximum OCR (1.25 and 2.12) were determined at values covering 90% of all published data for the upper section ( Figure 4b). However, only several published data were available in the lower section so that we used the minimum and maximum OCR values (1.75 and 2.65) directly for the lower section (Figure 4b).

Correction for Stress Field
For the analysis of the samples from the accreted section, an additional correction is required because the lateral tectonic force changes the in situ triaxial stress state in the accretionary prism (e.g., a horizontal orientation and the overburden corresponds to the minimum principal stress (Figure 6), which would lead to an underestimation of the pore pressure. To overcome this problem, J. C. Moore and Tobin (1997) directly interpreted pc   as the horizontal maximum principal effective stress in the prism. In contrast, building on the observation that porosity loss is caused by the mean effective stress increase ( Figure 6), we calculated the mean effective stress at yield point using pc   . Because the consolidation test is widely regarded as K 0 consolidation and has a constant stress ratio (K 0 ) of horizontal to vertical effective stress, the mean effective stress can be estimated by: where mpc   is the mean effective stress at the yield point. If we assume the intermediate principal effective stress equals the minor principal effective stress in accretionary prism, the in situ mean effective stress can be calculated by: where m   is the in situ mean effective stress and K in is the in situ stress ratio of minor to maximum principal effective stress. Assuming the mean effective stress at the yield point in the consolidation tests ( mpc   ) is equal to the in situ mean effective stress in accretionary prism ( m   ), the following expression is obtained when v   is replaced by σ v −P: Then the pore pressure in accretionary prism is estimated by transforming Equation 6 In our calculations, we used a K 0 = 0.81 for consolidation tests (Morgan & Ask, 2004) and K in = 0.5 for in situ conditions  that were experimentally determined for accretionary prism samples from Sites 1174 and 808, respectively. Both K 0 = 0.81 andK in = 0.5 are similar to the calculated stress ratios (K 0 = 0.90, K in = 0.59) based on the numerical model by Flemings and Saffer (2018) using a sediment friction angle of 15° (Casey et al., 2016). To evaluate the effect of stress ratios on pore pressure estimates, we also calculated the maximum and minimum K 0 and K in values with a sediment friction angle of 5-30° after Flemings and Saffer (2018) by the numerical model. Stress ratios of K 0 = 0.64-0.99 and K in = 0.33-0.84 were used to quantify the uncertainty of estimated pore pressures.

Numerical Estimates of Pore Pressure
To get a continuous insight into pore pressure variations at Site C0023, we compared our results with a recently developed method for pore pressure prediction by Flemings and Saffer (2018). In contrast to previous approaches, this method is built on the Modified Cam Clay (MCC) critical state model under plane stain conditions (Roscoe & Burland, 1968) and allows also the prediction of pore pressure in the prism sediments.
The method contains two models for the accreting sediments: an intact sediment model and a weak fault model. Both models postulate that the maximum principal stress is horizontal whereas the minimum principal stress is vertical. The intact sediment model further assumes the accretionary prism is in a state of Coulomb failure and the stress state is fully controlled by the material's inherent properties, that is sediment friction angle and the porosity-mean effective stress relationship. Based on critical state soil mechanics, the pore pressure can be estimated by: where ϕ is the friction angle of sediment material, e is Euler's number, e in is the void ratio, λ is the slope on the   ln m   −e in plot for accreting sediments at Coulomb failure, and e λτ is the intercept at m   = 1 MPa (see Flemings & Saffer (2018) for details). In contrast, the weak fault model hypotheses that the stress state is affected by the properties of faults such that the accretionary prism reaches Coulomb failure in advance due to preexisting weak faults or is not at Coulomb failure (Flemings & Saffer, 2018). Therefore, the friction angle of weak faults is introduced in the weak fault model and the pore pressure can be predicted by: where ϕ wf is the friction angle of weak fault, e λwf is the intercept at m   = 1 MPa on the   ln m   -e in plot for weak faults at Coulomb failure.
For the underthrust section, the method assumes the sediment is in a state of K 0 consolidation and it is converted into the widely used empirical logarithmic equation to calculate vertical effective stress by porosity (e.g., Skarbek & Saffer, 2009). The prediction equation with soil mechanic meaning is expressed as: where e λK0 is the intercept at m   = 1 MPa on the   ln m   −e in plot for underthrust sediments at K 0 consolidation state.
Assuming the elastic deformation is negligible, the yield surface in mean effective stress-maximum shear stress space represents the iso-void ratio stress path (Wood, 1990), which makes the stress variation predictable when the sediment deforms under undrained condition. The MCC critical state model shows how to derive the   ln m   −e in relationship at any stress state when that at one specific state is known (e.g., Azizi, 2000). In our calculations, we considered e λK0 = 1.16 and λ = 0.463 inferred from reference Site 1173 by Flemings and Saffer (2018) as known parameters. A sediment friction angle of 15° (Casey et al., 2016) and a weak fault friction angle of 5° (Flemings & Saffer, 2018) were used to derive K 0 , e λτ , and e λwf (see Flemings & Saffer (2018) for details). By transforming the shipboard porosity data at Site C0023 (Figure 2b) to void ratio, the pore pressure profile was obtained.

Coupled Pore Pressure Dissipation and Chemical Transport Model
A vertical one-dimensional coupled pore pressure dissipation and chemical transport model was established to constrain the pore pressure estimates along the major reverse faults. Specifically, we tested the hypothesis if the obtained pore pressure profile was affected by lateral fluid flow as indicated by the observed Cl − anomaly between 600 and 780 mbsf.
A constant pore pressure gradient equivalent to the averaged pore pressure ratio (λp = P/σ v ) of 0.71 as described in Section 4.1 was assigned to the model as initial pore pressure profile. For the initial Cl − concentration, the general decreasing trend (Figure 2c) was assigned. Based on the fault-valve model (Sibson, 1990;Sibson et al., 1988), the fault permeability increases immediately after fault activity due to the inherent roughness of natural rupture surfaces, which leads to the updip migration of overpressured fluids. The high pore pressure along faults, in turn, further enhances and maintains the fault permeability through reducing the effective stress (e.g., Fisher & Zwart, 1997). Under sufficient fluid source condition, the discharged fluids can be compensated by releasing fluid from the source and high pore pressure can exist for long periods of geologic time (Bredehoeft, 2009). In the Muroto transect, active smectite-to-illite transformation (Steurer & Underwood, 2003) in the downdip direction provides a potential fluid source for the maintenance of high pore pressure. Hence, we induced the duration of transience into our model, which represents the timescale of transient lateral fluid flow and was simulated by maintaining the lithostatic conditions over an extended period before we allowed pore pressure to dissipate. Modeling results have shown that pore pressures along the fault can rise rapidly to near lithostatic value (e.g., Bekins et al., 1995;Saffer & Bekins, 1998). We therefore simulated the transient lateral fluid flow by increasing pore pressures along faults to lithostatic and decreasing Cl − concentrations to certain low values. However, this does not mean the lithostatic pore pressure exists throughout the subduction thrust since the fault zone is heterogeneous (Brown et al., 1994). Four temporarily active fault zones were assigned in our model because (1) we observed four low Cl − concentration spikes in chlorinity profile, and (2) several reverse faults and shear fractures were found around each low Cl − concentration spike. The timing of the fault activities was assigned separately based on the width of the spikes, given that older fault activity has more time to allow advection and diffusion of solutes and thus has wider spikes (e.g., Saffer & Bekins, 1998). When the transient lateral fluid flow ceases, the lithostatic pore pressure along faults starts to dissipate. We used the concentration profile of the inert Cl − to constrain the pore pressure dissipation, given that the vertical fluid flow away from the fault zone drives the advection of solutes. Constant pore pressures and Cl − concentrations were specified for both the top and bottom boundaries of this model under the assumption that transient fluid flow can only affect limited horizons within the timescale covered by the model due to low permeability.
The pore pressure dissipation is quantified by the differential equation of consolidation (Terzaghi, 1943): where t is the time, c v is the coefficient of consolidation, and z is the depth. We used a c v = 0.0008 cm 2 /s for outer trench-wedge and trench-to-basin transitional facies, c v = 0.0028 cm 2 /s for Upper Shikoku Basin facies, c v = 0.0032 cm 2 /s for faulting sediments of Lower Shikoku Basin facies and c v = 0.0007 cm 2 /s for underthrust sediments of Lower Shikoku Basin facies, calculated from our experimental data.
For the variation of Cl − concentration in a water-saturated one-dimensional sediment column, we used the following differential equation (Boudreau, 1997): where C is interstitial solute concentration, that is, Cl − concentration, D is the effective chemical diffusion coefficient, and v is velocity. D is determined by D = D 0 /(1−log(n 2 )) (Boudreau, 1997), where D 0 is estimated after Boudreau (1997) using the temperature profile reported for Site C0023 in Heuer et al. (2017) and the porosity, n, is based on a robust polynomial regression of the porosity profile (n = -5.776 × 10 −15 ·z 5 + 1.539 × 10 −11 ·z 4 −1.410 × 10 −8 ·z 3 + 4.997 × 10 −6 ·z 2 −6.595 × 10 −4 ·z + 0.485, adjusted R 2 = 0.98). The velocity, v, is given by Darcy's law: where k is the vertical permeability, computed by porosity based on the empirical equation (i.e., log(k) = −20 + 5.5 n) by Saffer and Bekins (1998) and η w is the temperature-dependent viscosity of pore water after Huyakorn and Pinder (1978). The porosity and permeability profiles were assumed constant in our model because of the limited time interval of 3 kyr covered by the model (Section 4). We further neglected porosity rebound during transient unloading because of the small porosity change associated with it (e.g., Nobes et al., 1991). Also, the constant vertical permeability is consistent with the view that lateral fractures only affect horizontal permeability (Saffer, 2015).
Previous studies showed that the duration of transient lateral fluid flow varies from several years (e.g., Hiramatsu et al., 2005) to thousands of years (e.g., Henry, 2000;Howald et al., 2015), we therefore conducted sensitivity analysis using a duration of 10-1,000 years for the transient lateral fluid flow. Meanwhile, we also tested an influx fluid with Cl − concentrations from 350 to 430 mM. The best parameters which predict our experimental pore pressure estimates were determined and constraints on the timing of the fault activities were acquired by matching the modeled results to the observed data after a certain time step.

Experimental Results and Numerical Pore Pressure Estimates
The compression index, Cc, and recompression index, Cr, are the slopes of the virgin elasto-plastic compression and elastic recompression lines on the logarithmic mean effective stress versus void ratio plot (Knappett & Craig, 2012). Our consolidation tests yielded Cc and Cr values mainly in the range of 0.11-0.34 and 0.01-0.03, respectively.
In this study, the modified Casagrande method (McNulty et al., 1978) and Pacheco Silva method (1970) produced similar pore pressure estimates. The pore pressure difference between these two methods is lower than 0.7 MPa and the maximum difference as a percentage is 5.6% (Figure 7c). For the convenience of comparison with other studies, only results by the modified Casagrande method are described.
The consolidation tests conducted on four accreted sediment samples from Upper Shikoku Basin facies to outer trench-wedge facies (Unit III to IIB) yielded pore pressure in excess of hydrostatic (Figure 7c). The experimentally determined excess pore pressure varies from 1.3 MPa at 491 mbsf to 2.3 MPa at 568 mbsf. The pore pressures increase linearly with depth with a slope larger than hydrostatic, which corresponds to a narrow pore pressure ratio range of 0.69-0.76 (Figure 7d). In comparison, our numerical approaches predict similar pore pressures for the weak fault model, whereas in the absence of weak faults, the intact sediment model results in higher pore pressure ratios (λp = 0.71-0.85).
For the underthrust sequence, the experimental data also revealed pore pressure in excess of hydrostatic (Figure 7c). A total of five samples were tested for the subducted Lower Shikoku Basin facies with excess pore pressures ranging from 3.4 MPa to 4.2 MPa between 843 mbsf and 1,022 mbsf. The pore pressure ratios are in range of 0.71-0.75 (Figure 7d). Although the excess pore pressures for subducted sediments are higher than those for the accreted section, their corresponding pore pressures follow the same increasing trend with depth. Again, our modeling results show good agreement with the experimental data.
Further, our experimental results indicated high excess pore pressures of 5.7 MPa and 5.0 MPa for the samples from 714 and 721 mbsf, respectively (Figure 7c). These values correspond to pore pressures supporting 95% and 89% of the overburden stress. The depth of the high excess pore pressures coincides with that of low Cl − anomalies and the upper boundary of a depth interval with enhanced number of reverse faults and shear fractures (Figure 7b).

Coupled Pore Pressure Dissipation and Chemical Transport Model
Sensitivity analysis for our coupled pore pressure dissipation and chemical transport model shows that the maintenance of excess pore pressure is sensitive to the duration of transient lateral fluid flow ( Figure S3). Lithostatic pore pressure dissipates quickly when using a duration of 10 years, but can be nearly maintained for several thousand years if a duration of 1,000 years is employed. Along the major reverse faults, a duration of 10 years generates obviously lower pore pressure estimates compared to experimental data, while the duration of 1,000 years yields obviously larger values. Only the modeled pore pressure from a duration of 100 years best fits the experimental data.
In contrast, the Cl − concentration profile is controlled by both the duration of transient lateral fluid flow and the Cl − concentration of influx fluids ( Figure S3). A longer duration of transience generates more freshened fluid for the same Cl − concentration, whereas a higher Cl − concentration of influx fluids leads to less freshened water for the same duration of transience. For a duration of 10 years, the modeled Cl − concentrations are overall larger than the measured values. On the contrary, a duration of 1,000 years predicts obviously lower values than the observed Cl − concentrations. The modeled Cl − concentrations using a duration of 100 years are generally consistent with the measured values and reach best fit when a moderately freshened fluid with Cl − concentration = 400 mM is used.
In addition, the modeling results show that the assigned multiple fault zones with different timing of activity match the observable amplitudes of the Cl − concentration peaks (Figure 7e). Our model predicts both ZHANG ET AL.

10.1029/2020JB020248
12 of 19 the experimental pore pressure estimates and the Cl − concentration profile best when assuming (1) four separate fault zones which were active between 3 kyr and 0.5 kyr before present, (2) lithostatic pore pressure is maintained for 100 years before we allow pore pressure to dissipate, and (3) a moderately freshened fluid with Cl − concentration = 400 mM ( Figure S3).

Comparison with Previous Pore Pressure Estimates
Here, we provide new estimates for the excess pore pressure at the seaward edge of the Nankai Trough accretionary prism offshore Muroto peninsula based on a combination of uniaxial consolidation tests and numerical modeling. In contrast to many previous studies focusing on the underthrust sequence, we report on pore pressure in both the accreted and the underthrust sediments with regard to the consolidation state and stress path. The Cc and Cr values generated in our experiments are consistent with previous results of 0.32 ± 0.15 and 0.03 ± 0.02 for the sediments at the Muroto transect (Bellew, 2004;Feeser et al., 1993;Morgan & Ask, 2004), which might be a good hint for the quality of our experiments. Our experimental estimates of pore pressure are reliable considering the error bar is generally narrow.
Our experimental data reveal an overpressured underthrust sequence with excess pore pressures ranging from 3.4 MPa to 4.2 MPa (Figure 7c). The pore pressure results are consistent with previous estimates from uniaxial consolidation experiments using a triaxial cell (Morgan & Ask, 2004), corroborating our correction scheme for secondary consolidation and diagenesis. The laboratory-based estimates are also in good agreement with our own and other numerical modeling results that utilized Site 1173 as a normally consolidated reference site (e.g., Hüpers & Kopf, 2009;Screaton et al., 2002;Skarbek & Saffer, 2009;Tobin & Saffer, 2009). As discussed in Section 3.1.2, the porosity profile at Site 1173 reflects the combined effects of primary, secondary consolidation and diagenesis. Hence, those effects are accounted for in the vertical hydrostatic effective stress-porosity relationship of the reference site, which produced results similar to those obtained by our correction scheme. This confirms previous assumptions by Screaton et al. (2002) that the vertical hydrostatic effective stress-porosity relationship observed at Site 1173 continues to depths equivalent to those at Sites 1174 and 808.
Our consolidation tests also reveal excess pore pressure in the accreted Upper Shikoku Basin facies to outer trench-wedge facies (Unit III to IIB). The estimated pore pressures are higher than previous estimates at Sites 808 and 1174 by Tsuji et al. (2008) based on seismic velocity. Their mean effective stress-porosity relationship derived from underthrust sediments does not consider the different stress state for the accreted sediments and causes the underprediction of excess pore pressure. Our experimentally determined pore pressure ratios meet the theoretical limits for in situ stress at Site 808 inferred from borehole breakout by Huffman and Saffer (2016). In comparison, our numerical approaches predict similar pore pressures only for the weak fault model, in which the strength of the prism sediment is governed by low friction along preexisting faults or the accretionary prism is not at Coulomb failure. The intact sediment model results in higher pore pressure ratios (λp = 0.71-0.85) than in the underthrust sequence (Figure 7d). Using the same numerical approach with friction angles of 5°-30°, Flemings and Saffer (2018) obtained similar values (λp = 0.7-0.9) for Sites 1174 and 808 and proposed the excess pore pressure was attributed to rapid loading by burial and horizontal compression. Although uncertainty may exist in some of above estimates, all these data are much higher than the fluid overpressures of <0.2 MPa monitored in a sealed borehole observatory at Site 808 (Sawyer et al., 2008). Flemings and Saffer (2018) raised the possibility that small leaks in the observatory bleed excess pressure more rapidly than fluids can be delivered from the surrounding clay-rich wallrock.

Pore Pressure Regime in Shallow Subduction Zones
In general our data is consistent with the view of rapid loading by burial and horizontal compression causing elevated pore pressure in the shallow low permeable clay-rich subduction zone sediments (e.g., J. C. Moore & Vrolijk, 1992;Saffer & Tobin, 2011). However, our results provide new insights into pore pressure regime and timing of fluid flow processes in the shallow subduction zone, although we made some assumptions and simplifications in our calculations and modeling (Section 3).
The experimental results show a continuous pore pressure increase across the accreted and underthrust sequence with λp = 0.69-0.76, except for the two samples in the major reverse fault zone. This does not support previous views that the décollement acts as a fluid flow barrier and causes excess pore pressure in the underthrust sediments (e.g., Gamage & Screaton, 2006;Saffer & Tobin, 2011). Although previous studies have reported that fault-normal permeability can be several orders of magnitude lower than the wall-rock (e.g., Faulkner & Rutter, 1998), neither the permeability is low enough, nor is the major reverse fault zone sufficiently thick to seal the fluid pressure according to the theoretical fluid pressure barrier model (Deming, 1994). This is also corroborated by the measured permeability of sheared sediment samples at Site 1174 (Ikari & Saffer, 2012), which is 1∼2 orders of magnitude higher than the estimated permeability required for a fluid barrier (Gamage & Screaton, 2006). Therefore, the porosity offset across the shallow décollement reflects different consolidation processes above and below the décollement, instead of a fluid barrier. Our experimental work therefore supports the view of Flemings and Saffer (2018) that the additional porosity loss in the accreted section results from a higher mean effective stress and shear-induced compaction caused by the tectonic loading.
Our experimental data do not support previous speculations by Flemings and Saffer (2018) that pore pressure in the accreted section is higher than in the underthrust sediment as predicted by the intact sediment model ( Figure 7c). However, if near lithostatic pore fluids dissipate from the major reverse fault zones, the pore pressure can rise up in the wall rock and exceed our estimates based on consolidation experiments, as shown in our coupled pore pressure dissipation and chemical transport model ( Figure 7c). Therefore, our experimental estimates represent the lower bound of in situ pore pressure. On the contrary, the intact sediment model may predict the upper bound because it assumes constant friction angle and compressibility (λ). Recent studies have shown that the friction angle (Casey et al., 2016) and compressibility (Casey, 2014) decrease with increasing stress. These two overestimated values at high stress will lead to overestimation of pore pressure. The weak fault model shows a good agreement with experimental data, hinting that the stress state in the prism might be controlled by the presence of preexisting weak faults or that the accreted intact sediment is not at Coulomb failure (Flemings & Saffer, 2018;Huffman & Saffer, 2016). An outer prism not in Coulomb failure would also be consistent with the dynamic Coulomb wedge theory that predicts stress relaxation of the prism in the post-seismic interval in which the sediment is brought back to a stable elastic regime (Wang & Hu, 2006).

Near Lithostatic Pore Pressure along Décollement
For the top of the major reverse fault zone at 714 and 721 mbsf, near lithostatic pore pressure is revealed from the experimental tests, which we relate to fault activity as discussed in Section 3.3. However, unloading by transient lithostatic pore pressure should not be reflected in the pore pressure estimates by consolidation tests. We speculate that the disturbance from deformation in the major reverse fault zone destroys the fabric of sediment and causes memory loss of past maximum effective stress. This speculation is supported by the optical microscopic observation from Maltman et al. (1993) that intense fabric reorientation is visible in the shear zone at Site 808. Furthermore, these authors also observed hydraulically fragmented breccia at 788 mbsf, which provides evidence that lithostatic pore pressure ever occurred in this horizon. If the hydraulic brecciation is present in the major reverse fault zone at Site C0023, the fabric of sediment might be further destroyed. With the most conservative estimation of excess pore pressures without any correction for secondary consolidation and diagenesis (i.e., using Equation 7 with OCR = 1), excess pore pressures of 4.2 and 5.4 MPa are still higher than that of surrounding sediments and correspond to 84% and 93% of the lithostatic stress, respectively (Figure 7c). The high excess pore pressures are also supported by the low Cl − anomalies at corresponding depths, which are common signals for the transient and localized updip migration along a temporarily dilatant fault in response of lithostatic pore pressure (Bolton et al., 1998;Tobin et al., 2001).
In addition, based on boundary conditions and assumptions described in Section 3.3, our coupled pore pressure dissipation and chemical transport model predicted the observed Cl − anomaly ( Figure 7e) and the remaining present pore pressure was similarly high to experimental estimates (Figure 7c), which we also see as a prove for the reliability of the fault zone pore pressure estimates. Along the décollement, the density change during unloading is small (e.g., Nobes et al., 1991), but the high pore pressure may cause a decrease in the seismic velocity (e.g., Tobin et al., 1994) and lead to a reverse-polarity of seismic reflections (e.g., J. C. Moore et al., 1995). Therefore, our data support that the impedance-decreasing décollement found by Park et al. (2014) could be at least partially attributed to the increase of pore pressure within the fault zone. The reflection of décollement in the seismic profile (Figure 1b) might represent a highly overpressured horizon, which also marks the migration path of freshened pore fluids.
Further, the coupled pore pressure dissipation and chemical transport model shows that the high excess pore pressures can be maintained when a duration of 100 years for the transient lateral fluid flow is employed ( Figure S3). This duration, in which both the lithostatic pore pressure and low Cl − concentration are maintained at the fault, is similar to the fault sealing rate (70 years) used by Gratier et al. (2003). Instead, it is considerable shorter than previous estimates for the Nankai margins of 80 kyrs (Saffer & Bekins, 1999). We explain the difference with the number and the ∼100 m wide distribution of the fault strands as well as relatively small excursions assumed to be attributable to updip fluid migration (Saffer & McKiernan, 2009) in our model. This greatly reduces the diffusion time compared to previous models (Saffer & Bekins, 1998, 1999 which attribute all the geochemical anomalies to episodic fluid flow and consider a single thick fault strand. Also, our results of frequently active fault strands are in line with studies on biomarker thermal maturity indicators that identified seismic faults in drill cores recovered from the Japan Trench subduction zone (Rabinowitz et al., 2020). These results show that even fault zones that have hosted earthquakes with displacement ≥10 m ruptured repeatedly during a period of 36 kyrs.
The view of transient fluid flow along the subduction thrust is not at odds with the proposed continuous pore pressure increase across the accreted and underthrust sequence, because excess pore pressure along the faults will eventually reduce to the background pore pressure in times of inactivity. Due to the temporal proximity of the last fault activity, the lithostatic pore pressure may not have fully dissipated to the background values. Such an overpressure is not captured by previous laboratory data and any numerical methods applied in this study or previous ones, because they do not account for porosity change by unloading through transient pore pressure increase. Based on the best fit of our coupled pore pressure dissipation and chemical transport model, excess pore pressure could be high in the depth interval between 600 and 780 mbsf.
Setting a constant pore pressure gradient (equivalent to λp = 0.71) and the general decreasing trend of chlorinity ( Figure 2c) as initial conditions, the best fit of the Cl − profile in our coupled pore pressure dissipation and chemical transport model only requires a moderate fluid freshening. This led us speculate on the fluid source depth. While low Cl − values are commonly associated with deep fluids, typical geochemical signatures for a deep-seated origin, such as elevated B, Li or CH 4 concentrations (e.g., Kastner et al., 1991;You et al., 1996), are missing in the respective horizons at site C0023 (Heuer et al., 2017). We therefore relate the fluid origin to an overpressured zone along the décollement downdip at a depth of advanced clay dehydration but below peak formation temperature of thermogenic hydrocarbons (100°C-150°C). These assumption contrasts the previous view of episodic fluid flow along the entire fault surface after earthquake rupture but supports the emerging picture of a heterogeneous plate boundary with overpressured and permeable patches along the fault surface that shift in time and space (Saffer, 2015). A possible trigger for the spatially limited fluid flow could be the propagation of the frontal thrust which is rooted to the subduction thrust close to Site C0023 (Figure 1b).

Mechanical Strength of Décollement and Its Implication
Because pore fluid pressure is a critical parameter that determines the mechanical strength and slip behavior of faults (Saffer & Tobin, 2011), we computed the shear strength using the effective stress at 714 and 721 mbsf and a friction coefficient of 0.27 for Muroto décollement reported by Ikari and Kopf (2017).
The resulting shear strength is as low as 0.4-0.8 MPa, lower than the 2∼3 MPa previously estimated by Tobin and Saffer (2009) and also lower than the ∼2 MPa suggested by Kopf and Brown (2003) as well as Skarbek and Saffer (2009). Instead, it is similar to the inferred coseismic shear strength (0.54 MPa and 0.22-1.32 MPa) in the 2011 Tohoku earthquake from temperature anomaly after the event by Fulton et al. (2013) and from high-velocity friction experiments without excess pore pressure by Ujiie et al. (2013), respectively. Although we do not have pore pressure estimates for the décollement, the excess pore pressure there should be equally high because lithostatic pore pressure transients migrated along that horizon (Figure 7c). Faulkner et al. (2011) proposed that propagation through material with very low shear strength is energetically favorable as it provides little resistance and may explain how a large rupture is not arrested from propagating to the trench as in the Tohoku earthquake. Moreover, overpressured shallow décollement zones with low effective stress are conditionally stable and can be easily destabilized by a sudden dynamic slip (Park et al., 2014). Laboratory data by Saffer and Marone (2003) further show that under low effective stress conditions the tendency for velocity strengthening behavior diminishes. Therefore, if the reverse-polarity of seismic reflections along the décollement horizon offshore Shikoku Island (Park et al., 2014) represents similarly high pore pressure as revealed at Site C0023, seismic slip could be propagated to the trench and cause tsunami hazard offshore SW Japan when a large earthquake occurs.

Conclusions
We characterized the pore pressure regime and fluid flow processes in the shallow Nankai Trough subduction zone based on robust pore pressure estimates from consolidation experiments and numerical models. From our results, excess pore pressures increase continuously across the accreted and underthrust sediments, with nearly constant pore pressure ratios. We suggest that both accreted and underthrust sediments constitute one hydrogeological system and the décollement is not a fluid barrier as previously believed. Near lithostatic pore pressure along the subduction thrust seems to be attributable to the propagation of the frontal thrust which causes repeated lithostatic pore pressure transients migrated along faults in the last 3 kyrs. The resulting extremely low shear strength is similar to the inferred coseismic shear strength in the 2011 Tohoku earthquake, which triggered massive near-trench displacement and a catastrophic tsunami.

Data Availability Statement
Data reported in this study are available in supporting information and PANGAEA (https://doi.pangaea. de/10.1594/PANGAEA.918818). Data used in Figures 2 and 7 are available in Heuer et al. (2017). Data used in Figure 4 are available in Morgan and Ask (2004), Bellew (2004), Hüpers and Kopf (2012), Guo and Underwood (2014), and Kitajima and Saffer (2012), and have been collected in supporting information.

Acknowledgments
The data and samples used in this research were provided by the IODP. The manuscript benefitted from the discussion with Demian Saffer and Matt Ikari. We greatly appreciate constructive comments by Editor Douglas Schmitt and two anonymous reviewers. Junli Zhang thanks China Scholarship Council for supporting his research at University of Bremen. Open access funding enabled and organized by Projekt DEAL.