Evolution of the geological structure and mechanical properties due to the collision of multiple basement topographic highs in a forearc accretionary wedge: insights from numerical simulations

We propose a conceptual geological model for the collision of multiple basement topographic highs (BTHs; e.g., seamounts, ridges, and horsts) with a forearc accretionary wedge. Even though there are many BTHs on an oceanic plate, there are few examples of modeling the collision of multiple BTHs. We conducted numerical simulations using the discrete element method to examine the effects of three BTH collisions with forearcs. The typical geological structure associated with a BTH collision was reproduced during the collision of the first BTH, and multiple BTH collisions create a cycle of formation of BTH collisional structures. Each BTH forces the basal décollement to move up to the roof décollement, and the roof décollement becomes inactive after the passage of the BTH, and then the décollement moves down to the base. As the active décollement position changes, the sequences of underthrust sediments and uplifted imbricate thrusts are sandwiched between the décollements and incorporated into the wedge. At a low horizontal compressive stress, a “shadow zone” is formed behind (i.e., seaward of) the BTH. When the next BTH collides, the horizontal compressive stress increases and tectonic compaction progresses, which reduce the porosity in the underthrust sediments. Heterogeneous evolution of the geological and porosity structure can generate a distinctive pore pressure pattern. The underthrust sediments retain fluid in the “shadow” of the BTH. Under the strong horizontal compressive stresses associated with the next BTH collision, pore pressure increases along with a rapid reduction of porosity in the underthrust sediments. The distinctive structural features observed in our model are comparable to the large faults in the Kumano transect of the Nankai Trough, Japan, where a splay fault branches from the plate boundary and there are old and active décollements. A low-velocity and high-pore-pressure zone is located at the bottom of the accretionary wedge and in front (i.e., landward) of the subducting ridge in the Kumano transect. This suggests that strong horizontal compressive stresses associated with the current BTH collision has increased the pore pressure within the underthrust sediments associated with previous BTHs.

subduction zones. Seamounts are typical BTHs and common features on the global ocean seafloor (Wessel et al. 2010). Faulting along trenches creates horsts and grabens, which occur along most global subduction zones, as indicated by bathymetric data (Kobayashi et al. 1998;Masson 1991;von Huene and Ranero 2003). Subduction of BTHs has been reported in several subduction zones, such as in the Nankai Trough (e.g., Park et al. 2003;Bangs et al. 2006), New Zealand (e.g., Barker et al. 2009; Barnes et al. 2010), Costa Rica and Nicaragua (e.g., von Huene et al. 2000;Ranero and von Huene 2000) Sumatra (e.g., Singh et al. 2011), and northeast Japan (Tsuru et al. 2000). The deformation of the overriding plate caused by the BTHs is recorded in the seafloor bathymetry, such as the indentation and collapse of slopes, and uplift of accretionary wedges (e.g., Lallemand et al. 1989;Park et al. 1999;Pedley et al. 2010). Subducted BTHs may lead to enhanced sediment subduction, which can locally elevate pore pressures, or result in km-scale spatial variations in composition along plate boundary faults that can produce variable fault slip behavior . The fault slip may involve large megathrust earthquakes, slow slip events (SSEs), low-frequency earthquakes (LFEs), very-low-frequency earthquakes (VLFEs), and tremors (e.g., Cloos 1992;Bell et al. 2010;Ikari et al. 2013;Wang and Bilek 2014).
In this study, we focus on the effects of collisions of multiple BTHs, particularly small-scale BTHs that are buried in sediment in accretionary margins, and not erosional sediment-starved trenches. Collision of multiple BTHs can lead to complex geodynamic processes that reconstruct and overprint existing geological structures formed by a previously subducted BTH. Furthermore, even small-scale seamounts buried in sediment can generate distinctive structures in overriding forearcs (Morgan and Bangs 2017). Smaller seamounts are more common on oceanic plates (i.e., there are 25 million small-scale seamounts that are > 100 m in height worldwide; Wessel et al. 2010). This implies that smaller seamounts frequently collide with subduction zones. In fact, multiple collisions of small-scale BTHs buried in sediment have been observed in natural subduction zones, such as the Hikurangi subduction margin Pedley et al. 2010) and Nankai Trough (Tsuji et al. 2013). Therefore, multiple collisions of BTHs are thought to be quite common, particularly in the case of smallscale BTHs.
Analogue and numerical experiments have revealed the long-term tectonic evolution of an overriding plate affected by a subducting BTH, which can involve subduction erosion, margin uplift, and formation of a frontal accretionary wedge and trench re-entrant (e.g., Dominguez et al. 2000; Morgan and Bangs 2017;Ruh et al. 2016). The heterogeneous mechanical and hydrological conditions obtained from numerical simulations, particularly those associated with the contrast between the leading flank and wake of the seamount, can explain observed patterns of megathrust slip (Ellis et al. 2015;Sun et al. 2020). However, these studies have only examined subduction of a single seamount, and there are few studies that have investigated the effects of subduction of multiple seamounts. Analogue modeling experiments of sequential subduction of two seamounts have demonstrated that a distinctive duplex structure develops along the leading edge of the subducting seamount (Wang et al. 2021). However, even in that study, the effects of subduction of the second seamount on the structure created by the first seamount were not evaluated. Consequently, the geological structures formed by the repeated subduction of BTHs are still unclear.
Here, we present a conceptual geological model for the collision of multiple BTHs, based on numerical simulations, and we then consider the mechanical and hydrological conditions in the collision zone. We modeled multiple BTH collisions that involve three BTHs, and documented changes in mechanical properties, such as stress and porosity. Our study highlights how multiple BTH collisions reform the pre-existing geological structure and mechanical properties, which can subsequently control variations in fault slip behavior in the wedge.

Model settings
We used the discrete element method (DEM), which is a particle-based numerical simulation (Cundall and Strack 1979). DEM modeling has reproduced the complex geological structures of subduction zones (e.g., Miyakawa et al. 2010Miyakawa et al. , 2019Morgan 2015). Our two-dimensional experimental set-up follows that of Miyakawa et al. (2010). The initial model conditions represent a stack of coherent sediments on an oceanic plate (Fig. 1). The thickness of the sediments was set to 1000 m without water loading (i.e., a free surface). In this simulation, only solid deformation is calculated, and no fluid was incorporated into our model. We used two sets of particle parameters (Table 1). These parameters were selected from biaxial compression tests that simulate the internal friction angle of strong material (39°) and weak material (27°). The radii of the particles were variable: 12-30 m for particle A and 6-10 m for particle B. Layer A is the upper part of the sediment and consists of strong particles (particle A) that represent sediments and/or sedimentary rocks, whereas Layer B is the lower part and consists of a thin layer (100 m) of weak and small particles (particle B) that represent a weak layer and/or décollement where deformation is localized over a thickness of Page 3 of 13 Miyakawa et al. Progress in Earth and Planetary Science (2022) 9:1 several tens of meters (e.g., Mascle et al. 1988;Taira et al. 1991). Given that the initial model layers reflect incoming sediments, which are not fully cemented, we did not consider inter-particle bonding in our simulations. Rigid walls were used to confine the model boundaries, and friction between the walls and particles was the same as the inter-particle friction for particle A. The simulation was run using PFC2D software (Itasca Consulting Group Inc 2014). A moving wall was used to compress the sediments and simulate the accretion process. The displacement rate was small enough to approximate the deformation as quasi-static and, therefore, was set to 0.0009 m (0.09 cm) during each calculation cycle. The apparent deformation rates in the numerical simulation were derived from the shortening rate of the subduction zone and moving wall. The shortening rate is 4 cm/yr in the Nankai Trough, southwest Japan (Seno et al. 1993), whereas the shortening rate of the moving wall was 0.09 cm per time-step in our simulations. Therefore, to calibrate the shortening rate of the model with that observed in nature, the timescale was assumed to be 0.0225 yr/step (~ 8 d/step). The moving wall shortened the sediment layers by 5000 m, which is equivalent to a time interval of 1.25 Myr.
We performed simulations using two types of model settings: (1) a reference model; (2) a BTH collision model. The reference model contains no BTHs (Fig. 1a), but three BTHs were incorporated into the BTH model (Fig. 1b). For simplicity, we used a single size and shape for the BTHs (150 m in height and 300 m in width), and the spacing between the BTHs was constant (20,000 m). We modeled small BTHs that are buried in thick sediments (1,000 m). These BTHs were constructed by fixing the movement (i.e., their velocity and rotation) of the particles already set in the model. Given that the particles were fixed to the bottom wall of the model, they penetrated the moving wall during shortening of the model, and finally escaped from the modeled system.

Measurement of porosity and stress in the DEM model
We investigated the mechanical and hydrological evolution of the wedge by measuring the stress and porosity in the simulated models. These physical quantities were estimated from particle assemblies within a circular measurement region (radius = 270 m). The size of the measurement circle was set to cover several tens of particles in order to robustly calculate the porosity and stress. The measurements were undertaken for the entire accretionary wedge model by setting the measurement circles at 90-m-spaced intervals, both in the vertical and horizontal directions. We calculated the porosity and stress from the simulations according to the following. The porosity, n, was defined as the ratio of total void volume (i.e., area in this 2D simulation), V void , within the measurement region to the total volume of the measurement region, V reg : where V mat is the volume of particles in the measurement region. The stress in a measurement region was computed by averaging the contact force and vectors that connect the centroids of adjacent contacting particles in the measurement region (Christoffersen et al. 1981). The average stress, σ ij , in a measurement region of volume V was calculated as: where N c is the number of contacts in the measurement region or at its boundary, F (c) is the contact force vector, L (c) is the branch vector joining the centroids of the two bodies in contact, ⊗ denotes the outer product, and the compressive stress is positive. These measurements were conducted with the "Measurement" command in the PFC2D software (Itasca Consulting Group Inc 2014).
We monitored the porosity and stress in the wedge by tracking the same particle in the simulation. The particle tracks observed in the DEM simulation imitate the tracks of the sediments in the accretionary wedge. The initial location of the monitoring particle was set in the strata behind the first BTH (i.e., ~ 22,000 m in the horizontal direction and ~ 450 m in the vertical direction from the bottom of the left-moving wall; Fig. 1). We then measured the porosity and stress within the measurement circle centered on the monitoring particle at every 200,000 steps (i.e., 180 m of shortening).
The stress state and deformation control the change in porosity in the model. The porosity in this study is the volume fraction of particles, and thus it is strongly affected by the mean stress. However, the porosity is not always identical to the mean stress, because porosity is also affected by changes in particle arrangement due to deformation. These processes lead to a porosity decrease due to compaction and a porosity increase due to faulting and fracturing.

Geological structures caused by multiple BTH collisions
A simulation without BTHs provided evidence that the model set-up can reproduce typical wedge structures in subduction zones (Fig. 2). The reference model reproduced the typical structure of a fold-and-thrust belt, including the self-similar wedge geometry and uniformly spaced foreland (left-side) dipping thrust that propagates from the basal décollement ( Fig. 2; Additional file 1: Movie S1). The position of the décollement remains at the base of the wedge throughout the whole simulation. The typical geological structure associated with a seamount collision, which has been generated in previous studies, was reproduced in the BTH collision model during the first BTH collision (BTH1) (e.g., Dominguez et al. 2000;Morgan and Bangs 2017;Noda et al. 2020;Fig. 3; Additional file 2: Movie S2). As BTH1 collides, it forces the basal décollement fault (F1) (shortening of 10,440 m in Fig. 3a) to move up to the boundary fault between the overlying thrust sheets and incoming underthrusted sediments that contain the monitoring particle (shortening of 15,480 m in Fig. 3b). The BTH produces a "shadow zone" that protects the strata in the wake of the BTH (Dominguez et al. 2000;Noda et al. 2020;Wang et al. 2021). The protected, undeformed strata are underthrusted beneath the overlying imbricate thrust sheets (e.g., shortening of 15,480 m in Fig. 3b). As the BTH subducts, the increasing load due to the overlying imbricate thrust sheets and the increase in friction between the overlying sheets and subducting strata lead to the formation of a new basal décollement at the base of the incoming sediments (Dominguez et al. 2000). In our simulation, after the passage of BTH1, the roof décollement becomes inactive (F1), and the most active fault moves down to the basal décollement (F2) (shortening of 15,840 m in Fig. 3c). The underthrust sediments that contain the monitoring particle remain at the bottom of the wedge (shortening of 23,220 m in Fig. 3d). The upward movement of the roof décollement and uplift of the imbricate thrust lead to a surface slope with a high angle (shortening of 15,480 m in Fig. 3b), but the downward movement of the basal décollement (i.e., the formation of fault F2) and formation of a new imbricate thrust return the surface slope to a low angle (shortening of 23,220 m in Fig. 3d).
Page 5 of 13 Miyakawa et al. Progress in Earth and Planetary Science (2022) Fig. 3b, e, i), and the roof décollement becomes inactive after the passage of the BTH, and then the active décollement moves down to the base of the model (shortening of 15,840 and 32,760 m in Fig. 3c, f ). Each BTH collision causes this change in the position of the active décollement. As the active décollement position changes, the sequences of underthrust sediments and uplifted imbricate thrusts sheets are sandwiched between the décollements and become incorporated into the wedge. We refer to the structural units bounded by the old roof and new basal décollements as BTH collisional units (Fig. 4). These BTH collisional units are stacked and preserved with each BTH collision. The décollement that forms the unit boundary is a long-active and large-offset fault. The angle of the former décollement becomes high as the accretionary wedge grows, and it is exposed at the seafloor as a unit boundary fault.

Mechanical development during multiple BTH collisions
The general trend of the stress and porosity evolution is cyclical and self-similar, as is the wedge growth, in the model without BTHs (Fig. 2; Additional file 1: Movie S1). The continuous horizontal compressional stress propagates from the moving wall to the toe of the wedge in the reference model (Fig. 2). This continuous horizontal compression forms uniformly spaced thrust faults. The formation of the thrust faults breaks the initial strata and, consequently, generates a high-porosity zone along the faults at their toe. The shallow part of the wedge has a high porosity due to the lower level of stress and deformation accompanying the faulting and fracturing. Compaction within the deep part of the wedge reduces the porosity as the wedge grows ( Fig. 5a; Additional file 3: Movie S3). The mean stress increases continuously around the monitoring particle, which is in the thrust sheet incorporated into the deep part of the wedge. As the stress level increases, the porosity decreases around the monitoring particle (shortening of 10,000 to 25,000 m in Stage PI in Fig. 5a). However, the porosity increases slightly when an existing fault is reactivated close to the monitoring particle (shortening of 30,000 to 35,000 m in Stage PI in Fig. 5a). The porosity increases abruptly and significantly when a minor out-of-sequence thrust is formed and reactivated around the monitoring particle at the end of the simulation (~ 40,770 m of shortening in Stage PII in Fig. 5a). Both the stress and porosity structure are strongly affected by the multiple BTH collisions (Fig. 3; Additional file 2: Movie S2). Until the BTH approaches, the stress field in the imbricate thrust sheet is characterized by horizontal compression similar to the reference model. The porosity becomes high along the imbricate thrusts, but remains low within the thrust sheets (Fig. 3).   Fig. 3b, e, i). This uplift-induced deformation is also observed in analogue models (Dominguez et al. 2000).
The changes in stress and porosity in the underthrust sediments around the monitoring particle are as follows. A high horizontal compressive stress is formed in front (i.e., landward) of BTH1, and a low horizontal compressive stress (i.e., a "shadow zone") is formed behind (i.e., seaward of ) BTH1 (shortening of 15,480 m in Fig. 3b). With an increase in the overburden due to the overlying imbricate thrust sheets, the mean stress increases even behind BTH1 (shortening of 15,480 m in Fig. 3b; Stage BI in Fig. 5b). The porosity in the underthrust sediments then decreases due to compaction under the high stress levels, even in the stress "shadow zone" (i.e., the porosity decreases from 0.15 to 0.14 or less; Stage BI in Fig. 5b). As BTH1 continues to subduct and the active fault moves down to the basal décollement (F2), the horizontal compressive stresses increase in the underthrust sediments and a low porosity is maintained (Stage BII in Fig. 5b). When BTH2 collides, the horizontal compressive stress increases and tectonic compaction progresses, which further reduce the porosity (0.135 to 0.130; Stage BIII in Fig. 5b). The porosity increases (from 0.130 to 0.140) at the end of Stage BIII, because the high horizontal stresses and associated high deviatoric stresses result in the formation of minor faults and fractures in the underthrust sediments. The horizontal stress in the underthrust sediments decreases dramatically as BTH2 passes through the model (shortening of 36,000 m in Fig. 5b; Stage BIV). In this stage, the small amount of uplift and fracturing associated with the passage of BTH2 causes the porosity to increase (from 0.140 to 0.150). After the passage of BTH2, the underthrust sediments again experience compressive stress at intermediate depths in the accretionary wedge (Stage BV). The compaction then gradually progresses and the porosity decreases (i.e., from 0.150 to 0.140; Stage BV). These changes in the mechanical properties of the underthrust sediments are repeated every time a BTH collides ( Fig. 3; Additional file 2: Movie S2). The repeated evolution of the mechanical properties in the underthrust sediments is shown in the schematic cross-section of the evolution of a wedge with multiple BTH collisions (Fig. 4).
The porosity increase due to a change in particle arrangement induced by deformation is larger than the porosity decrease due to the increase in mean stress. In the transition from Stage PI to PII in the model without BTHs, the porosity increases while the mean stress increases. At the same time, the mean stress and deviatoric stress increase, which leads to deformation around the marker. In fact, during the transition from Stage PI to PII, fault reactivation occurs around the marker, and the resulting change in particle alignment causes fractures and a significant increase in porosity. The same is true for the BTH model. Up to Stage BI-BII and early in Stage BIII, the porosity decreases as the mean stress increases. However, at the end of Stage BIII, even though the mean stress remains high, the porosity increases. At this time,  The stress state (mean stress, deviatoric stress, and differential of the horizontal to vertical stress) and porosity around the monitoring particles. a The reference model without BTHs (Fig. 2) and b the model for multiple BTH collisions (Fig. 3) Page 9 of 13 Miyakawa et al. Progress in Earth and Planetary Science (2022) 9:1 the deviatoric stress increases, which indicates that a fracture has formed around the marker and the porosity increase is due to the change in particle arrangement.

Hydrological effects of multiple BTH collisions
We investigated the hydrological conditions within the wedge based on the changes in porosity and stress field. Geological processes in the outermost forearc are dominated by disequilibrium compaction associated with rapid burial and lateral tectonic loading associated with plate subduction (Saffer and Tobin 2011). This compaction (i.e., porosity reduction) discharges pore fluid or increases the pore pressure. The high permeability of fractures and fault zones exerts a primary control on drainage patterns (Moore 1989;Carson and Screaton 1998). In contrast, the low permeability of the sediments undergoing compaction limits the ability of the fluids to access permeable conduits and facilitates the development of overpressures (e.g., Neuzil 1995;Saffer and Bekins 2002). Evolution of the porosity and formation of fractures and faults inferred from our model can provide insights into the hydrological conditions in the wedge. The heterogeneous development of geological and porosity structures in the BTH model can generate a distinctive pore pressure pattern. The pore fluid in the imbricate thrust sheet at the tip of the wedge is easily expelled along fractures caused by faulting and/or uplift-induced fractures that are associated with the seamount collision (Dominguez et al. 2000). In contrast, the underthrust sediments behind the BTHs can store pore fluid. The BTH carries underthrust sediments that are protected by the stress shadow into the deep part of the accretionary wedge (Stage BI in Fig. 4). The underthrust sediments do not undergo faulting or fracturing that would form dissipative flow paths for the fluid (Stages BI and BII). Therefore, the underthrust sediments retain fluids, which increases the pore pressure while the porosity is reduced (Stages BI to BII in Fig. 4). Under the strong horizontal compressive stress conditions associated with the next BTH collision, the pore pressure increases dramatically and there is a rapid reduction in porosity (Stage BIII in Fig. 4). The pore pressure in the underthrust sediments becomes highest at this stage. From the end of Stages III to IV, the collision of the BTH with the upper plate causes minor faulting and fracturing, which is expected to locally enhance permeability and lead to dewatering and fluid escape (Sun et al. 2020). After the passage of the BTH, the porosity increases and minor faults and fractures should dissipate the pore fluid and reduce the pore pressure (Stage BIV in Fig. 4). After the pore fluid is dissipated, the pore pressure does not increase when the strata are again subjected to the compressive stress field (Stage BIV in Fig. 5). This hydrological evolution is repeated by the multiple BTH collisions that occur each time another BTH is subducted. The pore pressure distribution is not only controlled by the presence of the subducting BTH, but also by the geological and hydrological record of previous BTH collisions.

A comparison of the simulation results with natural subduction zones
The structures caused by multiple BTH collisions in the numerical simulations are observed in natural subduction zones. The Kumano transect obtained across the Nankai Trough off the Kii Peninsula, Japan, is one of the highest quality 3-D seismic reflection datasets in the world, but the nature of the complex geological structures along the transect is controversial (e.g., Moore et al. 2009;Shiraishi et al. 2019). In the Kumano transect, several BTHs (i.e., ridges) associated with faulting within the subducting oceanic crust have been identified (Tsuji et al. 2013). The collision of these multiple BTHs has affected the complex geological structure of the Kumano transect (Fig. 6).
The most distinctive structure in our simulations of multiple BTH collisions is the old and new décollements that are sandwiched between the BTH collisional units (Fig. 4). These structural features are comparable to the large faults in the Kumano transect (Fig. 6), which are old and active décollements and a splay fault that extends upward from the plate boundary. The ridge forces the active décollement to extend out to the seafloor, although the décollement extends to the bottom of the wedge and is parallel to the top of the oceanic crust landward of the ridge. The old décollements are located above the active décollement and connect with the fault that extends out to the seafloor. The shape of the units sandwiched between the active and old décollements (Fig. 6) is similar to the BTH collisional units observed in our model (Fig. 4). Furthermore, the splay fault is imaged as a boundary between the low seismic reflectivity and highly deformed zone of the hanging wall and the coherent seismic reflectivity layer of the footwall (Shiraishi et al. 2019;Fig. 6). The adjacent highly deformed and coherent units across the splay fault in the seismic profile correspond to the adjacent highly deformed thrust sheets and underthrust sediments across the old décollements in our model (Fig. 4). Therefore, the splay fault may be an older décollement and the unit sandwiched between the splay fault and old décollement is a BTH collisional unit, which formed by previous collision of a ridge that has now been subducted.
The velocity structure in the Kumano transect can also be explained by the hydrological evolution inferred from our model. At the base of the wedge, a significant Page 10 of 13 Miyakawa et al. Progress in Earth and Planetary Science (2022) 9:1 low-velocity and high-pore-pressure zone developed along the old and new décollements (Park et al. 2010;Tsuji et al. 2014). The low-velocity and high-pore-pressure zone are located at the base of the BTH collisional unit and in front (i.e., landward) of the current ridge (Fig. 6). This scenario corresponds to Stage BIII in our model, where the strong horizontal compressive stresses associated with the current BTH collision increase the pore pressure within the underthrust sediments that were transported to depth by the previous BTH collisions (Fig. 4). The strata that are the pore fluid reservoirs were transported behind the BTHs that had previously undergone collisions, and the high pore pressure anomaly observed in the present wedge is caused by compression due to the current BTH collision.

Effects of multiple BTH collisions on earthquakes
We speculate that our results also have implications for the distribution of slow earthquakes. Pore fluid accumulation and, in particular, elevated fluid pressures have been linked to fault slip, including SSEs, tremors, LFEs, and VLFEs (Saffer and Tobin 2011). The mechanical and hydrological differences between the leading flank and stress shadow of a BTH undergoing collision could explain the observed pattern of megathrust slip, with earthquakes and microseismicity being concentrated at the down-dip edge of seamounts and aseismic or slow slip occurring in the up-dip stress shadow (Sun et al. 2020). However, tremors and SSEs occur in front of BTHs in some subduction zones, such as the Nankai Trough (Tsuji et al. 2014;Shiraishi et al. 2020) and Hikurangi margin (Barker et al. 2018).
By comparison with our model, we infer that BTHs have previously been subducted under areas where tremors and SSEs are distributed along the leading flank of a currently colliding BTH. The underthrust sediments transported in behind previous BTHs can supply a large amount of fluid and generate a high pore pressure in front of the actively colliding BTH (Stage BIII in Fig. 4). The presence of large amounts of fluid and high pore pressure on the leading flank of the actively colliding BTH can inhibit coseismic activity and causes SSEs. In fact, in the Kumano transect, VLFEs occur in front of the ridge (Nakano et al. 2018;Shiraishi et al. 2020), which is the area we interpreted to be underthrust sediments (Fig. 6). Past seamount subduction has also been identified at the Hikurangi margin, where tremors occur in front of the actively colliding seamount (Gase et al. 2021). Thus, understanding the variety of earthquakes associated with BTH subduction requires consideration of not only the actively colliding BTH, but also of previous BTHs.

Limitations of our fluid-absent model
No fluid was incorporated into our model, and thus this assumes that the fluid pressure acting on the solid matrix is negligible. However, pore pressure would decrease the effective stress, leading to decreased shear strength on faults, which then affects the mechanics of the accretionary wedge (e.g., Davis et al. 1983;Saffer and Bekins 2002). Therefore, the assumption of negligible pore pressure is a major limitation of the deformation modeling. To address this problem, advanced numerical techniques, which can simulate intrinsic coupling between fluid percolation and deformation of a solid matrix, need to be applied in future studies (e.g., Omlin et al. 2018).

Conclusions
We undertook numerical simulations to examine how multiple BTH collisions affect the geological structure of an accretionary wedge, and how mechanical properties, such as stress and porosity, evolve in the wedge. Each BTH collision forces the basal décollement to move up to the roof décollement, and the roof  Fig. 6 Interpretation of the seismic section extracted from 3-D seismic data along the Nankai Trough-Kumano transect from Moore et al. (2009). Colored units (gray, brown, green, and blue) represent each interpreted BTH collisional unit. The location of a very-low-frequency earthquake (VLFE) cluster is marked according to Shiraishi et al. (2020). The original locations of the VLFEs were determined by Nakano et al. (2018) Page 11 of 13 Miyakawa et al. Progress in Earth and Planetary Science (2022) 9:1 décollement becomes inactive after the passage of the BTH, and then the active décollement moves downward to the base of the wedge. As the active décollement position changes, the underthrust sediments and uplifted imbricate thrust sheets are sandwiched between the décollements and incorporated into the wedge as BTH collisional units. This distinctive structural feature is comparable to the old and active décollements and splay fault in the Nankai Trough, Japan. The stress state in the underthrust sediments changes from isotropic compression in the "shadow zone" behind the BTH to strong horizontal compression in the leading flank of the next BTH to collide. This stress state evolution causes a porosity reduction under undrained conditions and leads to a high pore pressure in the underthrust sediments. The high pore pressure on the leading flank of the colliding BTH can inhibit coseismic activity and causes SSEs. In fact, in the Nankai Trough, VLFEs occur in front of the ridge. Thus, understanding the earthquakes associated with the subduction of BTHs requires consideration of subduction of previous BTHs. The complexity of geological structures and types of earthquakes in subduction zones may be caused by not only the actively colliding BTH, but also the subduction of previous BTHs. These considerations are limited for the solid-based model and will be further investigated by solid-fluid coupled simulations.