Simple topographic parameter reveals the along-trench distribution of frictional properties on shallow plate boundary fault

The critical taper model best describes the first-order mechanics of subduction zone wedges. The wedge geometry, which is conventionally defined by two parameters, slope angle and basal dip angle, accounts for the strength of megathrust. By applying this theoretical model, fault frictional properties and earthquake occurrences can be compared among subduction zones, and within a single subduction zone, and the spatial distribution or temporal change of fault strength can be investigated. Slope angle can be accurately estimated from bathymetry data, but basal dip angle must be inferred from subsurface structure, which requires highly accurate depth-migrated seismic reflection profiles. Thus, application of the critical taper model is often limited by an insufficient number of highly accurate profiles, and the spatial distribution of frictional coefficients must be inferred from relatively few data. To improve this situation, we revisited the theoretical formula of the critical taper model. We found that the effect of basal dip angle on the critical taper model is small, and slope angle can be a proxy for the effective friction when the pore fluid pressure ratio is high, internal friction is small, or both. These conditions are met in many subduction zones. The validity of the approximation can be checked with a parameter newly introduced in this study. Therefore, this finding allows use of variations in slope angle, which could be obtained accurately from only the bathymetry as an approximation for relative variations in the effective coefficient of basal friction, if the targeted subduction meets the validity. We applied this approximation to the Japan Trench and estimated the variations in the friction coefficient distribution on the shallow plate boundary fault from 71 data points. We found that the area where the friction coefficient was smaller than the mean corresponded to a segment, where a large coseismic shallow rupture occurred during the 2011 Tohoku-Oki earthquake (Mw 9.0). Thus, by approximating tapered wedge geometry with a simple topographic parameter that can be obtained from existing global bathymetry, we can quickly estimate the distribution of frictional properties on a plate boundary fault along a trench and related seismic activity.


Introduction
The first-order mechanics of a subduction zone wedge, a representative feature of fold-and-thrust belts, can be clearly explained by the critical taper model (e.g., Dahlen 1990). This geomechanical model, which is based on the Mohr-Coulomb failure criterion, allows frictional properties on a plate boundary fault to be determined (Dahlen 1990). This model is a key method for understanding megathrust earthquake mechanisms, because direct measurements of plate boundary fault strength require drilling into the deep décollement to obtain samples and are quite rare (e.g., Chester et al. 2013;Ujiie et al. 2013). According to this model, the tapered wedge geometry (slope angle α and basal dip angle β) is determined by the strengths of the wedge materials and the effective friction on the megathrust fault (μ b ′) (Fig. 1). Thus, the critical taper model allows the geomechanical condition of a subduction wedge to be determined. This information can be used to compare geomechanical conditions among different subduction zones (e.g., Dahlen 1990) or to examine their spatial distribution within a single subduction zone (e.g., Fagereng 2011; Koge et al. 2014) or temporal changes along a single profile (e.g., Wang et al. 2010;Wang and Hu 2006). Slope angle α can be calculated from the bathymetry above a subduction wedge, which is typically observed by a multi-beam echosounder onboard a ship. In general, the bathymetry is obtained with a vertical error on the order of several meters. Even at a depth of 9000 m in the deepest trench, the errors are only about 0.5% (~ 45 m) (Nakanishi 2011). As the scale of subduction wedges is several tens of kilometers, the accuracy with which α is determined is sufficient for characterizing wedge geometry. On the other hand, the subsurface geometry parameter used in critical taper model calculations, namely, the basal dip angle β, requires further consideration. In a depth-converted profile of seismic reflection data, the depth to the plate boundary fault depends strongly on the velocity model used, and the accuracy of the depthconversion process affects the value of β. The frontal wedge of the Japan Trench, for example, is estimated to have P-wave velocity (V p ) of 1.7-3.4 km/s (Tsuru et al. 2002). We could estimate the value of β (height ÷ distance) to be 0.34 s ÷ 12 km from Nakamura et al. (2013). If we convert the second unit to kilometers with the existing minimum V p = 1.7 km/s, the height be 0.34 s × 1.7 km/s = 0.58 km. If with the maximum V p , the height be 0.34 s × 3.4 km/s = 1.16 km. Thus, the height has a width of 0.58 km. Considering this 12 km Slope angle estimated from bathymetry can be used as a proxy for friction of a plate boundary fault.
We can quickly estimate the distribution of frictional properties along a trench. α Applied to Japan Trench: The low frictional segment is consistent with the slip distribution of the 2011 Tohokuoki Earthquake.  Kimura et al. 2012). α: slope angle, β: basal dip angle, μ b ′ effective friction on the megathrust fault. The frontal wedge area is between the blue broken lines. B Diagram showing the self-similar growth of a bulldozer wedge (modified from Dahlen 1990). σ 1 : compressive principal stress, ψ 0 : the angle between σ 1 and the upper slope, ψ b : the angle between σ 1 and the basement wedge situation, the error range of height in a 20 km wedge, what we want to know, is ~ 1 km. This vertical error range is much larger than the error range of bathymetry, which is only a few dozen meters. As a simple case to consider the error scale of the value of α and β, we imagine a 20 km horizontal line with one side endpoint fixed (named A), and the other side endpoint (named B) can move vertically. If point B moves vertically by tens of meters, the change in slope angle of line AB only affects the sixth decimal place. On the other hand, for a vertical move of 1 km, the slope angle will change by up to 2.8°. These slope angle errors correspond to the error ranges of α and β, respectively. This simple exercise shows that, unless only highly accurate depth-migrated profiles are used to calculate the critical taper model parameter, comparisons within and among wedges are likely to be unreliable. On a scale of several kilometers, we can use pre-stack depth migration data or, at larger scale, a cross section of the velocity structure combined with refraction data to accurately determine β for the critical taper model. However, highly accurate pre-stack depth migration data on seismic reflection profiles survey are often not available, because they require long offset seismic data to improve the accuracy of the velocity model in deeper portions. As a result, the number of accurate cross sections available for a critical taper analysis is often insufficient to reveal detailed along-strike variations of frictional properties in subduction zones.

Revisiting, validating, and simplify critical taper theory
We first review formulations of Coulomb wedge/critical taper theory. All of the formulas are based on a noncohesive wedge model, which assumes non-viscosity (Dahlen 1984). According to the Mohr-Coulomb failure criterion, τ = σ · tanφ + C , where τ is shear stress, σ is vertical stress, φ is the internal friction coefficient (also expressed as μ, the coefficient of internal friction averaged over the wedge), and C is the cohesion force. Because internal friction forces are proportional to vertical stress, whereas cohesion forces are independent of vertical stress, the cohesion term can be neglected when considering huge geological structures with large σ. Thus, the noncohesive critical taper model is valid in the entire wedge.
Next, we theoretically verify the effect of the basal dip angle β on the calculation of effective friction μ b ′ and show that the effect of β becomes small when the pore fluid pressure in the subduction zone is high, which is the mean condition in trenches (Lallemand et al. 1994). Hence, we propose that basal friction in subduction zones can be reliably estimated from only the slope angle α determined from the bathymetry.

Revisiting critical taper theory: Overview of the critical taper model to obtain the effective coefficient of basal friction
Critical taper theory (Davis et al. 1983;Dahlen 1984;Lehner 1986) is a geomechanical model based on the Mohr-Coulomb failure criterion according to which the wedge geometry (α and β) is constrained by the balance between wedge strength and effective friction μ b ′ (e.g., Adam and Reuther 2000).
In critical taper theory, we obtain μ b ′ and the pore fluid pressure ratio (λ) in a wedge by drawing cross plots between λ and μ b ′, as explained below (e.g., Adam and Reuther 2000;Wang and Hu 2006;Wang et al. 2010Wang et al. , 2019 (Fig. 2). In the critical taper theory formulation, the modified slope angle αʹ under subaerial conditions is formulated as where α is a parameter obtained from the bathymetry/ seismic profile, ρ is wedge sediment density, ρ w is fluid density, and λ is the pore fluid pressure ratio. Then, the uniform angle between the most compressive principal stress axis σ 1 and the upper slope, ψ 0 (see Fig. 1), is calculated as where φ is the angle of internal friction within the wedge. Because along-strike stresses are not considered in the

Fig. 2
Cross plot between the pore fluid pressure ratio λ and basal friction μ b ' in the wedge. All extensionally critical states form the left limb of the critical state curve, and all compressively critical states form the right limb. The stable region is under the curve (white). The straight-line intersecting the critical state curve represents constant λ critical taper model, the following simple geometric relation is applicable ( Fig. 1): where ψ b is the angle between σ 1 and the basal plane. Then, the effective coefficient of basal friction (μ b ′) is obtained from the Mohr-Coulomb failure criterion τ and the stress balance of the basal condition as To draw the limb of the cross plot between μ b ′ and λ, we set λ to range from 0 to 1 (Fig. 2). The left limb of the critical state curve represents extensionally critical states, and the right limb represents compressively critical states (Table 1).
The effective coefficient of friction is obtained by creating an intersection point on this curve. In most subduction zones, fixing the fluid pressure ratio λ as a constant provides the constraint line for the intersection. Instead of such constant λ, in tectonic erosion, we could provide the constraint line μ b ′ = μ(1 − λ), which assumes that the basal fault is relatively strong and the wedge is in a compressively critical state (e.g., Wang et al. 2010;Wang and Hu 2006;Koge et al. 2014). However, the tectonic erosion equation also assumes constant φ. Therefore, for all these constraints on the intersection, it is difficult to discuss whether the spatial variation in effective frictional coefficients obtained by only the critical taper is due to the friction coefficient or pore fluid pressure.
For example, in the calculation with a constant λ, if we assume the mean wedge parameters ρ = 2700 kg/ m 3 , ρ w = 1000 kg/m 3 , internal friction angle φ = 34°, and λ = 0.88 (Lallemand et al. 1994) (Fig. 2), then μ b ′ can be determined from the intersection of λ with the critical state curve. Because the prism wedge in subduction zones should be in a constant compressively critical state just before failure, we focus only on the intersection with the right limb (representing compressively critical states). Thus, in this example, we obtain μ b ′ = 0.06. For more details than are provided in this simple review, please see the cited studies.

Validation and simplification: Effects of the geometric parameters on μ b ′
In this section, we examine the sensitivity of the calculated μ b ′ to the assumed α and β values to investigate how their accuracy affects the estimation of μ b ′. First, we used the mean wedge parameters presented in Sect. 2.1 and tested the sensitivity of μ b ′ to α and β. Afterward, we also tested a wider range of subduction zone parameters to explore the applicability of our approximation. First, we investigate variations in μ b in response to variations in α and β with the mean subduction zone parameter values. The states of the frontal wedge with α = 5° and β = 1° or β = 5° are shown in Fig. 3A; in Fig. 3B, both α and β are varied from 1° to 5°. We allowed β to range from 1° to 5°, because that range includes the basal dip angle of most subduction zones (Lallemand et al. 1994;Wang and Hu 2006). Here, since μ b ′ cannot be determined when α = β = 0, those results were removed.
We found that β has little influence on the estimation of μ b ′ when the pore fluid pressure is high, a consideration that had been neglected in previous studies. The change in β (from 1° to 5°) dominantly accounts for the change in the width of the critical state curve (i.e., the angle between its limbs) between the two states illustrated in Fig. 3A. When λ is high, the intersection between λ and the right limb of the critical state curve is near the curve peak. Therefore, the change in the width due to a change in β only slightly affects μ b ′. In typical subduction zones, λ is high (~ 0.88) and φ is 34º (Lallemand et al. 1994), so the effect of β should be regarded as small.
This finding is also favorable in terms of the accuracy of μ b ′ obtained by applying the critical taper theory, because the seismic profile depth used to calculate β depends on the velocity model/structure of the seismic profile, which is often not obtained with high accuracy. Moreover, the number of available profiles is also important to obtain the distribution of frictional properties by applying critical taper theory. Thus, because β must be obtained from depth-converted profiles with low accuracy, the resulting error is large (Koge et al. 2014). In contrast, α can be determined with negligible error. The seafloor depth, which is used to calculate α, is mostly based on multibeam data and the sound velocity profile (SVP) of seawater. Basically, the SVP can be obtained with high accuracy by conductivity/temperature/depth (CTD) measuring systems or by deploying expendable bathythermograph (XBT) or expendable CTD instruments. The observed SVP can be used to correct the water depth derived from two-way time to the seafloor; thus, SVP can be used to suppress the vertical error to on the order of several or a few dozen meters. The more significant factor influencing the error in α is whether the obtained cross section is aligned with the direction of maximum slope. However, if the deviation from the maximum slope direction is no more than 18°, the error in α will be less than 5% (see Additional file 1). Therefore, α can reliably be obtained with relatively high accuracy.
To evaluate the validity of neglecting β and approximation α-μ b ′ for use with the critical taper theory, we propose a parameter, which we call "weight of alpha" (WOA). To determine the sensitivity of the effective friction coefficient μ b ′ to α or β, we used linear multiple regression analysis, a statistical method that can be used to predict the value of a variable (the response variable) from the value of two or more other variables (explanatory variables). Here, we used α and β as explanatory variables and μ b ′ as the response variable. We conducted this analysis using the stats-model API in Python (Seabold and Perktold 2010). The obtained regression equation is characterized by the coefficients of α and β (A and B, respectively) and the coefficient of multiple determination R 2 . To evaluate the relative influence of α and β on μ b ′, we defined WOA as (WOA) = A/(A + B). For the mean subduction zone, we obtained R 2 = 0.99 and WOA = 78.02%; thus, the goodness of fit of the regression equation and the contribution of α to μ b ′ were both high. Although the contribution of β still remained 21.98%, which is not small enough to ignore when discussing each individual friction coefficient, it could be sufficiently small to discuss relative variations in friction coefficients along the trench.
Next, to determine whether the effect of β could be small in most subduction zones, we expand the parameters λ and φ to plausible ranges in nature (λ = 0-1, φ = 20-39°), and then calculated WOA and R 2 for all combinations of λ and φ in these ranges (Fig. 4A, raw data in Additional file 2). For example, we used mean subduction wedge values of 21 representative trenches (Lallemand et al. 1994), (λ, φ) = (0.88, 34°), and calculated a WOA = 78.02%. Of course, in some exceptional subduction zones, such as at the Makran Trench and Nankai-inner accretionary prism, λ is small, so the application of the α-µb' approximation might be inappropriate. Taking the toe of the Japan Trench as an example, the assumption of a bulk density of 2700 kg/m 3 for the wedge material is too high for sedimentary material. The bulk V p of the wedge material near the trench varies from 1.7 to 4.2 km/s (Koge et al. 2014). Given standard V p -to-density empirical relationships (e.g., Brocher 2005), this range would imply a bulk density of < 2000 kg/ m 3 . Furthermore, the range of α here is wider than those in typical subduction zones, taking values from 1° to 10°. Seismic refraction studies estimated β to be 3° to 6° in the shallowest part of the Japan Trench (e.g., Miura et al. 2003;Takahashi et al. 2004;Ito et al. 2005). The range of β does not exceed the range of α. A WOA diagram of the Japan Trench was calculated with the above parameters (ρ = 2000 kg/m 3 , ρ w = 1000 kg/m 3 φ = 20°-39°, α, β = 1°-10°) (Fig. 4B). At first glance, there are almost no changes in the two diagrams in Fig. 4. This is because, in the curve of critical taper as shown in Fig. 2, the change of the density from 2700 to 2000 kg/m 3 is only a slight parallel translation of the top of a trajectory in the lower right direction, resulting in only small changes in the calculated WOA values. Furthermore, note that the range of α and β is also larger in Fig. 4B, although the effect of these changes in ranges on the WOA diagram is small. Kimura et al. (2012) and Wang et al. (2019) assumed the following parameters for the Japan Trench: (λ, φ) = (0.95, 36°) and (0.8, 21.8°), respectively: under these conditions in the toe of the Japan Trench, WOA = 83.15% and 78.19%, respectively (Fig. 4B). Therefore, the variation of α could be a reliable proxy of the frictional variation along the Japan Trench.
Although we showed that our approximation of α-μ b ′ is applicable for a wide range of parameter values, the validity of the approximation should be checked carefully; for the sake of reliability, the method should be applied by creating a WOA diagram based on the parameters of each subduction zone (ρ, ρ w , α-range, β-range) and plotting the known λ and φ on this diagram.

Application example of the α-μ b ' approximation for the Japan Trench
Through our review and verification of the critical taper model, we found that α can be used as a proxy for the effective friction along the plate interface when WOA is high owing to either high λ or low φ (Fig. 4A). Thus, using our approximation, the relative distribution of the friction coefficient within a single subduction zone can be investigated by revealing spatial variations in the slope angle α estimated from only bathymetry data. Here, we apply this approximation to the Japan Trench, which is suitable for the application of this method, because, as shown in Sect. 2.2, WOA is high in the toe portion of the wedge and the friction coefficient along the plate interface can thus be estimated from only α. Furthermore, a coseismic rupture occurred near the trench axis during the 2011 Tohoku-Oki earthquake (Mw 9.0); therefore, by comparing the distribution of α, which should reflect that of μ b ′ in the Japan Trench, we can address the question of why such a coseismic rupture occurred in this specific area.
To apply our approximation model to the Japan Trench, we used the 250-m grid bathymetry data of Kishimoto (2000), focusing on the shallow portion within a horizontal distance of 25 km from the trench axis, and obtained the distribution of α, interpreted as the relative friction distribution, on the shallow megathrust (Fig. 5, Table 2). Using bathymetric profiles and slope angles obtained by . The closer WOA is to 100%, the more friction can be considered in terms of seafloor topography alone, because the effective basal friction μ b ' can be determined from α alone. A WOA diagram of the mean trenches (Lallemand et al. 1994). The white circles indicate the average value of these 21 subduction zones. The white boxes indicate the individual subduction zones, where the values were listed (Wang and Hu 2006). Conditions in the Japan Trench according to Kimura et al. (2012) and Wang et al. (2019) are shown by the two white squares refer the script as described by Wessel and Luis (2017), we accurately obtained the along-strike distribution of µ b ' on the shallow plate boundary fault at 71 points. This is a much larger data set than those used in typical studies that apply critical taper theory (typically, only a dozen points or fewer). The low-α segment (low-μ b ′ segment), located at approximately 37°-39°N, corresponds to the coseismic slip distribution of the 2011 Tohoku-Oki earthquake (Fig. 5;Lay. 2018). According to Lay (2018), the north end of the coseismic slip is still under discussion, but for the south end, many recent studies suggested that the large shallow slip did not reach 36°N; in particular, those studies using tsunami data indicated that the slip distribution was mostly limited to north of 37°N. In our study, we use the tsunami inversion of Saito et al. (2011), because it is most likely to be related to the deformation of the trench axis. The results suggest that large fault rupture occurred in the low-α segment, causing the slip to propagate to the shallow portion of the plate boundary fault because of the low friction there, resulting in the generation of the devastating tsunami (e.g., Ide et al. 2011). In contrast, the south and north ends of the coseismic slip zone are relatively high-friction areas. As a result, the slip could not propagate to these other segments, because the high friction there acted as barriers. Therefore, the low friction in the shallow area can be considered to be a major Comparison of the spatial variation of slope angle α in the Japan Trench and the coseismic slip distribution. A Compiled coseismic slip along the Japan Trench (red contours). The epicenter of the Tohoku-Oki earthquake is shown by a yellow star, and the red lines show the positions of bathymetric profiles used to obtain α. B Distribution of the slope angle α along the Japan Trench. The red area corresponds to the coseismic slip area in the map in A (Saito et al. 2011). The orange bar indicates the peak of the α distribution histogram. C North-south variation of the effective coefficient of basal friction from the seismic profiles (Koge et al. 2014). The red dots were determined by Tsuru et al. (2002) and the black dots were determined by Tsuji et al. (2011), and the JAMSTEC (Japan Agency for Marine-Earth Science and Technology) crustal structural database site facilitator of the huge tsunami. In addition, the low-α segment identified here approximately corresponds to the central segment along the Japan Trench (~ 37°N-39°N) inferred from the distribution of seismic activity detected by the S-net ocean-bottom seismograph network (Nishikawa et al. 2019). This area is the large slip region of the 2011 Tohoku-Oki earthquake, in which tremors and very low frequency earthquakes (VLFEs) have not been observed since 2011. Matsuzawa et al. (2015) suggested that this area had VLFEs clusters before the 2011 earthquake but became less active after the event. Low friction might have contributed to the non-occurrence of such earthquakes after the main shock, but the detailed mechanism requires more study. Low-friction conditions might prevail generally along the Japan Trench margin, except in the regions of high friction caused by the recent subduction of a seamount (Mochizuki et al. 2008) or the presence of petit-spot volcanoes (Hirano et al. 2006). Because μ b ′ depends on both λ and φ, it is not possible to determine whether variation in α (i.e., relative μ b ′) is due to a change in λ, φ, or both. Although here we cannot separate the effect on α of physical properties from that of pore fluid pressure, both effects are reflected in the strength of the megathrust. Koge et al. (2014) analyzed the Japan Trench with sparse datapoints considering non-approximated critical taper theory. They used seismic profiles by Tsuru et al. (2002), Tsuji et al. (2011), and the JAMSTEC (Japan Agency for Marine-Earth Science and Technology) crustal structural database site (https:// www. jamst ec. go. jp/ jamst ec-e/ IFREE_ center/ index-e. html) to estimate the basal dip angle β ( Fig. 5; red dots were determined by Tsuru et al. (2002), others are black dots). For the JAM-STEC data, Koge et al. (2014) converted the depth with the P-wave velocity of Tsuru et al. (2002). They found three locations with high frictional strength: the north and south ends of the coseismic slip area of the 2011 earthquake and locations near the epicenter. Our study increased the resolution of along-strike variations in the frictional strength with the approximated theory. Our updated results show that the overall effective frictional coefficient in the Japan Trench is low, although that beyond the north and south ends of the coseismic slip area of the 2011 earthquake is high, as documented by Koge et al. (2014). However, the high frictional strength near the epicenter of the 2011 earthquake documented by Koge et al. (2014) was not observed in our results. Although we estimated β for the north and south ends of the coseismic slip area of the 2011 earthquake based on Tsuru et al. (2002), we estimated β near the epicenter of the 2011 earthquake based on Tsuji et al. (2011) due to lack of relevant data in Tsuru et al. (2002). Tsuji et al. (2011) and JAMSTEC used different processing methods than Tsuru et al. (2002), and uncertainty in depth conversions of seismic profiles may cause biases in estimating frictional strength, as we noted in the introduction. In particular, using the migrated-profile D6 of Park et al. (2021), which is the nearest transect to that of Tsuji et al. (2011), we obtained α and β values of 4.69° and 3.26°; we obtained μ b ' = 0.09 in D6, with the critical taper based on the parameters of Koge et al. (2014). Now, our updated results for the Japan Trench highlight the usefulness of our α-μ b ' approximation.

The available data sources for using the α-μ b ' approximation
Bathymetry data are usually acquired by research vessels. However, even when such data are not available, bathymetry can be obtained from ETOPO1 (Amante and Eakins 2009) or other global data sets from satellites, most of which are freely accessible, and with vertical errors on the order of 150 m (Varga and Bašić 2015); such errors are only slightly larger than the errors in data acquired from research vessels. For higher resolution global terrain data, SRTM15 + V2.0 or GEBCO2021 have become available (Tozer et al. 2019; GEBCO Compilation Group 2021). However, in general, the terrain data combined with direct and indirect measurements such as SRTM15 + V2.0 or GEBCO2021 could include artifacts, where each data overlaps. Therefore, before using such satellite data, artifacts, such as patches or line structures, should be checked. Since the landward slope angle α corresponds to large wavelengths in the topography, low-resolution data should be sufficient for our new approximation. Therefore, whichever of these data sets is used, the frictional properties now can be easily evaluated.

Conclusions
We found a simplification of the critical taper model that provides reasonably accurate results under certain conditions. The spatial variation of the coefficient of the plate boundary faults could be approximated to the variation of the landward slope angle, provided that λ is high, φ is low, or both. This approximation method of the critical taper model intriguingly has the potential to advance our ability to characterize basal friction along shallow plate interfaces in subduction zones.
First, we reviewed the critical taper model formulas used for calculating the effective coefficient of basal friction μ b '. We found that in most subduction zones, the effect of β on basal friction can be regarded as slight. And more, we proposed a parameter WOA to show the validity of our α-μ b ' approximation. Our approximation seems largely valid when WOA is high, which occurs when λ is high, φ is low, or both. The