Bridging the gap between subduction dynamics and the long-term strength of the Sunda megathrust

Domain: earth_science
Published: 2025-11-28
Editorial outcome: minor_revision

Abstract
The link between great earthquakes and subduction dynamics across a wide range of scales remains a crucial, yet elusive, tenet of the seismotectonics of convergent margins. Here, we show high-performance computational simulations of the three-dimensional Sunda subduction zone dynamics matching plate motions, tectonics and current deformation, thereby providing insights into the stress regime of the Java-Sumatra-Andaman margin. Testing various tectonic forces reveal the primary control of the Java slab pull and induced mantle flow on the whole margin tectonics, driving northward-increasing oblique convergence, Sumatran trench advance and increased tectonic coupling along this segment. The modelled deviatoric stresses reproduce the geodetically-constrained interseismic compression in Sumatra and megathrust stress orientations notably consistent with the seismic P-axes, with magnitude comparable to static stresses on seismogenic faults. We map the frictional strength along the Sumatra-Andaman megathrust, identifying critical thresholds at key seismogenic depths, which remarkably correlate with the location of great earthquakes. Our outcomes show how subduction dynamics may critically prime the conditions for seismicity along the Sunda margin.

Introduction
The most devastating seismic events, classified as great earthquakes (M ≥8.0), occur along the plate interface of subduction zones, called the megathrust1. These earthquakes are explained within the framework of the subduction paradigm, which provides a context for stress distributions across scales, from the long-term tectonic loading over millions of years to the short-term stress released during seismic ruptures and its subsequent relaxation over years2. Although interconnected, these processes operate across vastly different spatio-temporal scales, making it a significant challenge for the scientific community, and society at large, to bridge our understanding across scales effectively.

The connection between subduction zones dynamics and their seismicity is exemplified by the megathrust strength, that is, the resistance of megathrusts to the action of tectonic stresses. The strength is the balance between the applied lithostatic and tectonic stresses, due to the mass distribution and the subducting slab-mantle interactions, and the intrinsic properties of the megathrust interface, due to fluids, temperature, lithology, and fault geometry, resisting (balancing) the applied stress. Constraining the strength of megathrusts requires integrating tectonics, structural features, geodetic and seismic strain across disparate scales. Geodetic imaging of the seismic and aseismic strain at major convergent zones shows that the strain accrued during the seismic and postseismic periods aligns with that accumulated between major earthquakes2, showing that the viscous behaviour of the lithosphere-mantle system controls the strain across timescales3, once the transient deformation due to earthquakes is relaxed. Therefore, several works have tested the controls subduction dynamics have on megathrusts seismicity, addressed within the context of tectonic coupling4, i.e., the tectonic stress across the plate interface. Statistical analyses have highlighted correlations between seismicity and slab age, that is, their inferred buoyancy5, upper plate strain6, a proxy for tectonic coupling, and mantle tractions that may enhance interface stress7. Similar analyses also propose the statistical relevance of physical properties of megathrusts acquired during subduction, such as low megathrust curvature8, increased roughness9 and thick sediment accumulation10. These works reveal the key role of the balance between boundary stresses acting on megathrusts and the internal stresses, although constraining these processes beyond statistical correlations has been proven difficult.

Additional constraints are provided by tests of plate kinematics and static force balance, although these do not necessarily constrain the tectonic stress acting on megathrusts. The plate kinematics frame characterising the seismic coupling does not describe the stress in the interseismic period, when this can be nonuniquely associated with either vanishing stresses or high megathrust strength4. This frame has been called into question, with long-held views on plate motions-relations confuted by recent seismicity11 and the controversial potential to lead to unrealistic high seismic coupling estimates12. Megathrust stresses and tectonic coupling, have been addressed directly, from static force balance calculations13 and its extension across stages of the seismic cycles14, and from thermal and mechanical properties constrained by heat flow15 and focal mechanisms16 at subduction zones. These processes are different and occur at different stages of the seismic cycle and, although providing useful insights, they do not support a unified physical model for the stress loading and release at subduction zones17.

The Sunda subduction zone serves as a pivotal case in seismotectonics, highlighting the challenges in understanding the complex interplay between subduction dynamics, long-term tectonics, and seismicity. This active convergent margin exhibits deformation consistent over the geological and seismic timescales18, making it a valuable case for the test proposed here. The Sunda margin, between the Indo-Australian and Eurasian plates, comprises three major segments, namely, from north to south, Andaman, Sumatra, Java and Banda. The largest seismic energy has been released along the Sumatran segment, which hosted three great earthquakes of magnitude Mw ≥8.5 within three years, among others and has highest coseismic strain of the region19. Notably, the 2004 Mw 9.3 Sumatra earthquake triggered a massive rupture stretching roughly 1300 km1,20, from northern Sumatra to the Andaman Islands (Fig. 1a), the longest along this margin21 (Fig. 1b). The estimated seismic coupling coefficient, i.e., the ratio between seismic and tectonic strain, in this area is close to unity22, suggesting that the deformation integrated over the historical catalogue and that due to convergence must be comparable. The compatibility between short- to long-term deformation is shown by the relative thickening of the Sumatran backarc, ongoing since ~15 Ma18, the advancing motions of this trench (Fig. 1a, red arrows), as opposed to the retreat in the other segments23, and matching GPS-constrained interseismic contraction in Sumatra19, that is, net of transients due to large earthquakes. Further evidence of this regime comes from flexural24 and seismotectonics analysis25 along the Sumatran segment, illustrating near-trench compression in the downgoing plate, extending this regime across the margin. The deformation accrued over million years and that constrained geodetically, show that the strain pattern is consistent across timescales, with the largest compression localised in the Sumatran segment, while further away from the trench, part of the deformation is accommodated by transcurrent faulting21,26.

The deformation along the Sunda margin squarely challenges the subduction paradigm. The direction of plate motions in this area varies progressively northward, from trench-perpendicular to nearly trench-parallel (Fig. 1a). This implies that the fastest convergence is along Java and the slowest along Andaman, which is at odds with the deformation distribution along the margin, with major compression and strike-slip in Sumatra and Andaman, respectively, where lower convergence does not reconcile with intense tectonics. How this can be explained by the subduction paradigm remains unclear: the strongest earthquakes and compressional tectonics are located where the subduction forces are inferred to be the weakest, challenging the idea of driving slab pull forces. Beneath central Sumatra, the upper mantle slab shallows up abruptly from >600 km deep, beneath Java, to only ~300 to ~400 km, then becoming <200 km deep beneath Andaman27–29. The age of the subducting lithosphere near Sumatra is the youngest, ~46 Ma, making it the thinnest, and the slab has the shallowest dip (Fig. 1a, b), while in Java it is ~150 Ma, and ~100 Ma in Andaman. For such a short, thin slab the lowest slab pull and ridge push forces18,30 maybe invoked, and yet this is where most intense seismic and tectonic activities occur. Additionally, the slab has a step beneath Sumatra, and an older and thicker slab is found beneath Java (Fig. 1a), reaching ~600 km depths, extending further east to the Banda Sea, where subduction is disrupted by the arrival of the continental Australian lithosphere31. Yet, the stronger slab pull expected from the Java slab is associated with reduced great earthquakes occurrences, that is, seismic energy release, and tectonic shortening. This enigmatic association, also observed in other margins where similar slab steps are present32, warrants further investigation.

Here, we propose a computational approach to characterising the stress and strain at the Sunda convergent margin and plate interface, the megathrust, as a result of the subduction force balance, providing insights into the seismotectonics of this area. We setup a numerical model of the Asian margin, from the Indian collision to the Sunda subduction zone, assimilating geophysical and geological observations from the Indo-Australian and Asian plates and the upper and lower mantle structures, to reproduce realistic terms of the force balance (Fig. 1). This approach has been proven to reconcile subduction dynamics with self-consistent stress and strain distribution across scales33, their relation with plate margin motions34 and magnitude of earthquakes35,36. We devise different setups to systematically test different tectonic forces, from that of the Sunda slab to that of the Indian mid-oceanic ridge, and from the force of the Indian collision to that of the lower mantle mass and flow. Our approach is consistent with previous work focused on seismic energy and centroid moment tensors37, plate motions and mantle flow38,39 and is here extended to include a self-consistent plate interface constitutive behaviour, critical to the realistic simulation of stress at convergent zones34,40. The models are validated against observed plate motions, upper plate tectonics, geodetically-constrained strain rates and compressional (P) axes from earthquakes centroid moment tensors catalogue, proving the consistency of the computed stresses. We focus on the deviatoric stress components, that is, the stress magnitude at which creep or frictional sliding occurs (Methods section). This represents the stress that is eventually accommodated over the long-term, net of the seismic transient, and allows the discussion of the long-term stress on the Sunda megathrust. Our outcomes support the discussion of the conditions for rupture and seismicity of the southeast Asian margin, bridging the gap between subduction dynamics, seismotectonics and the strength of megathrusts.

Results
Modelling the geodynamics of the Sunda subduction zone
We provide a comprehensive test of how fundamental tectonic forces control margin stress and strain by solving the force balance of subduction, encompassing the whole mantle and plate interactions, through high-performance computational modelling.

The fundamental subduction force balance in this region arises from the negative buoyancy of the Indo-Australian slab in the upper mantle, its variation due to the strong depth and thickness gradient along the margin, the resistance of the viscous upper mantle and the interactions with the Sunda-Asian upper plate. These features are embedded in the reference Model 1 using the oceanic crust age41 to infer age-dependent lithosphere thickness and ridge push (Fig. 1c), while the buoyancy of the continental lithosphere is reduced (see Methods). At depth, the slab thickness and extent in the upper mantle are reconstructed using the Slab2 dataset, embedding Benioff zone, seismic reflection and tomography data27. Including the southeast Asian margin only and a passive lower mantle from 660 to 1000 km depth (Fig. 1d and Supplementary Figs. S1 and S2) the Model 1 serves as a reference to assess the most fundamental force balance of slab pull and mantle tractions only and their gradients along the margin.

Additional forces, external to the subduction zone, may act on convergent margins altering the fundamental force balance. These far-field forces are the mid-Indian Ocean ridge push and the lower mantle drag, due to the sinking slab mass imaged by tomography in the region. The contribution of individual tectonic forces to the margin stresses are addressed by extending the modelling space to include different plate boundaries and mantle depths of the region in several setups (Fig. 1c). Model 2 extends to the east to the Banda Sea and tests the role of the subducting Australian continent increased buoyancy and mantle flow, other parameters are unchanged. Model 3 includes the whole tectonic plate from the Indian Ocean ridge to the India-Sunda trench, combining the full ridge push and Indian collisional forces with the pull of the southeast Asian and Indian slabs in the upper mantle. The boundaries between Indian, Capricorn and Australian plates are not implemented due to lack of observational constraints, although our lithosphere models embed the thinning of the lithosphere in the Capricorn zone, due to the varying oceanic crustal age (Fig. 1c and Supplementary Figs. S1 and S2). Model 3 extends to the core-mantle boundary, 2891 km deep, with a passive mantle in this setup, that is, has no density anomaly to drive flow, and extends to include the Indian continental and oceanic lithosphere, accounting for the Indian collision and ridge push force, respectively. In additional models, we vary the properties inside the Model 3 space and refer to these with a letter added to the numeral (3a to 3c, here). In these models we add the mass in the lower mantle, not included in the Slab2 dataset, using tomographic P- and S-wave models, MITP08 and Smean42,43, in Model 3a and 3b, respectively. These bracket end-members velocity anomaly distributions and smearing44, whereas the tomography’s anomalies in the upper mantle are neglected38. While different tomographic models exist, the models chosen have validated conversions to density contrast45, which critically reduces the complexity of our modelling (see Methods). Model 3c is the same of Model 3a, but the density of the Himalayan slab is increased to match the excess inferred driving forces on the Indian plate46 and lithosphere tractions arising from mantle flow39. We refer to these as active lower mantle models.

The modelled forces are balanced solving the conservation equations under the assumption of Stokes flow of an incompressible fluid, e.g., refs. 38,39. The force balance of the Stokes flow indicates that stress primarily depends on buoyancy contrasts, therefore gradients of stress along the margin are the result of the three-dimensional geometries and distribution of, most notably, slab mass in the mantle, while the material constitutive laws do not vary within the model space.

With the distribution of mass and the constitutive laws constrained, solving the conservation equations for Stokes flow yields self-consistent plate motions, mantle flow, and distributions of stress and strain rate (Fig. 1d). The models run for a small number of time steps, until subduction velocities reach stability, and we show results from a single time step, corresponding to a time of ~700 kyr. Under the assumption of the instantaneous solution, temperature variations due to advection and diffusion are neglected.

In our computational method, Lagrangian particles embedded in the finite element mesh improve the solution of the shearing along the margin, tracking accurately stress and deformation of different materials in the model domain. Additionally, a test of composite constitutive laws for the mantle shows that the stresses along the interplate margin, the focus of our work, are mainly dependent on plasticity, and less dependent on mantle creep. We have tested Newtonian and a range of non-Newtonian mantle creep laws (Methods), showing that the patterns of convergence velocities are mostly similar with all rheologies tested, while the magnitude varies only by a factor ~1.4 (Supplementary Fig. S3). Similarly, stress and strain rates around the plate margin keep consistent patterns, although the magnitude may vary up to a factor 4 and 10 for the stress and strain rates, respectively. The calculated values are compatible with the observables, although further constraining these values is limited by the uncertainties in the constitutive laws47. In the reminder of this work, we only show the results of Newtonian mantle models. This approach is commonly followed and capture the fundamental force balance involved34,39,40,48, while minimising computational costs. Tests the free-slip top boundary condition against a pseudo-free surface49, show negligible differences in our model set up (Supplementary Fig. S4). Also, varying the modelling space size imposes an inherent constrain on the large-scale mantle flow. However, this only alters the flow locally to the far walls and we show later that the impact on the flow around the slab is negligible.

Controls on strain and stress of the Sunda subduction zone
We present the six models described above to investigate on the role of the forces of subduction and mantle flow in the current deformation and tectonic coupling of the Sunda margin. Our results show that the models reproduce plate convergence, trench motions, and strain consistent with what observed along the Sunda margin. We show that the strain and stress fields at the subduction zone are primarily controlled by the attached slab mass distribution and the upper mantle flow, while lower mantle flow and far-field forces affect mostly plate motions and upper plate deformation patterns away from the trench.

We show in Fig. 2 the result of Model 1, where only the slab and upper mantle are modelled, and illustrate the complex flow arising from the gradients in slab depth, increasing the tractions beneath Sumatra. In the Java segment, the flow is coupled with the downgoing slab and it is mostly entirely poloidal. Instead, further north, the flow is complicated by the varying slab depth, and the mantle is forced around the slab edge, upwelling from the sub-slab to the mantle wedge region beneath central Sumatra (Fig. 2a). This mantle flow pattern is compatible with shear wave splitting and receiver functions analysis, showing the upwelling mantle from the Sumatran sub-slab to the mantle wedge50. We show that the escaping mantle from the Java sub-slab region result in the largest vertical velocity and, therefore, in excess trench-ward tractions beneath the Sumatran upper plate segment (Fig. 2b). Additional tests show that other forces, namely, continental subduction along the Banda segment and mantle flow around the slab (Model 2), ridge push and passive lower mantle (Model 3), and lower mantle slab mass and flow (Model 3a, b, c), do not alter the increased upwelling and tractions beneath Sumatra, and only affect the large-scale flow field, far from the margin (Supplementary Fig. S5).

A critical test of the models consists in the ability of the modelled subduction dynamics to reproduce realistic plate and trench motions. All models reproduce the characteristic northwards-increasing oblique convergence along the Sunda margin and oblique motions along the fore-arc12 (Fig. 3a, Supplementary Fig. S8). The variable error is minimised in models with lower mantle mass, in different reference frames, with Model 3b and 3c best aligned with plate motions in the GSRM1.2_SU, that is Sunda-fixed reference frame (Supplementary Fig. S9), showing the anchoring effect of the sinking slabs in the lower mantle. The models’ trench motions are compared to the observed trench velocities in the seafloor spreading-aligned51 (SA) and in the No-Net Rotation23 (NNR) reference frames, which bracket the end members proposed for realistic trench motions including and excluding net-rotation, respectively. The trench motions in both reference frames agree on the fastest advance of the northern Sumatra trench, up to ~25 mm yr–1, progressively decreasing outward to trench retreat of –20 mm yr–1 in Andaman, whereas bracket Java trench motions between +18 and –18 mm yr–1 (Fig. 3b). This pattern is coherently reproduced by all models, although the magnitude varies. Java trench retreat velocities are compatible with all models’, except Model 3, where the passive lower mantle likely enhances trench retreat. The modelled trench velocities align well with observed values in Andaman, whereas the best matching models for the Sumatran advance are those with an active lower mantle, constrained by tomography, Model 3a, b, c, while Models 1 and 2, reproduce the pattern but underestimate the trench advance by ~50%. Further constraining the magnitudes of the trench motions may be prevented by the modelling simplifications adopted. This is mostly due to the lack of a full torque around plates52 and a global mantle flow38, and the choice of boundary conditions.

The stress solutions in our models define regimes compatible with the Sunda margin tectonics (Fig. 3c). Although the upper plate rheology does not include plasticity, tectonic regimes can be defined by the stress tensor’s eigenvectors in the absence of explicit faults, indicating thrusting regimes, under which reverse faulting occurs (blue), and transcurrent regime (orange), under which strike-slip faulting occurs. Our results reveal persistent thrusting in the forearc, a characteristic feature of subduction margins, in agreement with the observed thickened lithosphere and crust18,53, geodetic strain19 and thrusting in western Sumatra54 and Java55. The region of transcurrent regime in our models extends from the Andaman back-arc to central Sumatra, in agreement with the extent of a major strike-slip fault system, which has been estimated to accommodate up to ~30% of the total deformation21,26.

The strain rates in the upper plate along the subduction zone, further illustrate the first-order control of the deep Java slab on the increased tractions along the Sumatran segment. The interactions between the slab’s negative buoyancy and upper mantle flow, only components of the force balance in Model 1 (Fig. 4a), explain the contraction along the trench and forearc, peaking at ~–10–14 s–1 in the Java-southern Sumatra segment, decreasing further north. These strain rates follow the gradient in the slab mass along the margin, with largest slab pull in the south, where slab reaches ~600 km (Fig. 4a, labelled and colored contours). Contraction rapidly vanishes in the back-arc in the southern and northern domains of the convergent margin, however, the contraction remains strong in Sumatra, O(–10–15 s–1), extending to the downgoing plate, in the forebulge area. This further supports the role of far-field forces, that is, the additional pull of the neighbouring Java slab and the increased mantle tractions. Note that the choice of different mantle and upper plate creep laws may affect the magnitude of the strain rates, although do not alter their distribution.

Far-field forces have negligible impact on the strain patterns at the margin and only influence the upper plate strain far from the margin. Model 2 shows minor variations in strain rates in the Java back-arc compared to Model 1 (Fig. 4b), following the inclusion of the Banda arc slab segment and Australian continent subduction. Model 3 has a ridge push and a whole, yet passive, mantle, and shows similar strain pattern to Model 1, although the entire upper plate is under extension (Fig. 4c). This is reversed in Models 3a, 3b, and 3c where lower mantle density derived from seismic tomography drives active mantle flow, enhancing far-field trench-ward tractions (Fig. 4d-f). The MITP08 tomographic model in Model 3a results in the same overall pattern as Models 1 and 2, with values slightly higher in the back-arc and lower in the Indian ocean lithosphere (Fig. 4d). Similarly, the Smean model in Model 3b only induces increased contraction in the upper plate, far from the trench (Fig. 4e). Finally, the increased Indian slab pull in Model 3c does not reveal any significant difference (Fig. 4f). The strain pattern at the margin and the broad strain peak around Sumatra remains consistent in all models tested, as well as the extension in Andaman, except for Model 3b, while the deformation in the far-field in Java remains poorly constrained by observations.

The models’ dilatation patterns reproduce the observations constrained by geodetic measurements along the Sunda margin19, further attesting for consistency (Fig. 1b). The dilatation dataset excludes data from stations or periods affected by transient phenomena, reproducing the patterns that best represent the interseismic deformation19, this is the background deformation over the whole catalogue periods, comparable to the tectonics18. Our analysis focuses on major areas where modelled rates agree with the observed dilatation (Supplementary Fig. S6), extending from the trench to the back-arc in north Sumatra and east Java, east of 110˚E, and in small domains of the Andaman basin, while overestimated strain rates in the domains in between may be due to the incomplete nature of the data or to our simplified Sunda plate model. The dilatation patterns are comparable to the second invariant of the strain rates (Supplementary Fig. S7), providing a similar match to the observables. In a profile of the observed and modelled dilatation (Fig. 4g), we show dilatation rates of ~–10–15 s–1 along a profile at 400 km from the trench, matching to the first-order the observations, with highest contraction in the Sumatra and Java. Model 2 and 3, which incorporate the Australian continent subduction along the Banda segment and the ridge and passive lower mantle, respectively, predict extension north-east of Java, suggesting these must enable trench migrations. Conversely, all models except Model 3a and b capture the Andaman back-arc extension. The area under extension remains rather small in the Models 1 and 2 and is more pronounced in Models 3 and 3c, where slab subduction and ridge push and lower mantle, respectively, are included. Our models accurately capture the pronounced extension observed between the compressional regions of the Sumatra and Java segments, although whether this reflects a tectonic feature remains speculative.

The stress patterns at subduction zones further emphasise the primary control of the slab pull heterogeneities in the upper mantle, whereas additional far-field forces only impact the magnitude of the back-arc stresses. The horizontal stress invariant at the surface, τ II, in all models shows high values of ~150 MPa in the Java and southern Sumatra forearc and decrease further north. In the back-arc, τ II is relatively higher in the central Sumatra to Andaman compared to its southern parts (Fig. 5a-c). The addition of an active lower mantle, that is, when flow is sustained by the tomography-constrained densities, increases the τ II magnitude in the back-arc (Fig. 5d-f), as shown by the increased compressional strain rates. The regional stress patterns show increased stress around Sumatra in both the overriding and downgoing plates. The arc profile shows high τ II values of 100-120 MPa from Java to southern Sumatra, progressively decreasing northward to values <40 MPa (Fig. 5g). The increased stress in the downgoing plate support the role of Java slab pull forcing, through the rigid downgoing lithosphere, not easily explained by mantle tractions alone.

The higher τ II in the Java to southern Sumatra forearc is almost independent of the far-field forces and the lower mantle mass, as all models follow the pattern of Model 1, with the slab only. Instead, the back-arc profiles, located 400 km from the trench (Fig. 5h), show lowered stresses, <40 MPa, however, these are of comparable magnitude all along the margin, which is explained by the controls of the dynamics of the Java subduction in the upper mantle on the whole region, as shown for the strain rates. Additionally, the stress patterns are mostly similar, although the magnitude increases in Model 3a-c, revealing that the modelled lower mantle mass can increase the magnitude of the stresses, but has minor impact on their distribution, primarily set by the upper mantle slab buoyancy gradients.

Methods
We model subduction in a regional spherical space as the Stokes flow of an incompressible fluid within finite Prandtl number. Under these assumptions, the equations of conservation of mass and momentum reduce to:

∇·u=0 (1)

–∇p+η∇2u+ρg=0 (2)

where u is the velocity vector, p the pressure, η the viscosity, and ρg is the force owing to density variations and (downward) gravity acceleration, respectively. In the Supplementary Information we provide a derivation of the equation of momentum from the force balance of subduction. We neglect the energy conservation equation, as we solve for the instantaneous flow, which inhibits the effect of thermal diffusion. The continuity equations describing the Stokes flow in their non-dimensional form are resolved using the particle-in-cell finite element code Underworld71.

The stress tensor σ is defined as the sum of the isotropic and deviatoric components:

σ= –pI+σ0 (3)

where I and σ0 are the identity tensor and the deviatoric stress, respectively.

The constitutive law relating the stresses and strain rates is given by the Newtonian creep of the form:

σ0=2ηε_ = η (∇u+(∇u)T) (4)

where ε_ is the strain rate tensor.

The viscosity η is defined by constitutive law that embeds the creep and plastic flow. The creep law considers the effect of both diffusion and dislocation creep, that is, Newtonian and non-Newtonian creep:

η creep = η 0 (ε_ II / ε_ T) (1–n)/n (5)

where η 0 is the reference (Newtonian) viscosity, ε_ II is the (square root of the) second invariant of the strain rate tensor, n is the power-law exponent and ε_ T is the transition strain rate, i.e., the value at which the dominant deformation mechanism switches from linear to non-linear creep, i.e., diffusion to dislocation creep. We have tested various mantle creep laws from purely Newtonian, i.e., n=1, to non-Newtonian, i.e., n=3, and tested the role of different transition strain rates. The transition strain rate controls the area of linear vs. non-linear deformation. For lower values tested, the mantle is predominantly under non-Newtonian flow, accommodating stresses by dislocation creep. With increasing values of ε_ T, the distribution of dislocation creep flow reduces to high stresses areas in the mantle, whereas the rest is under diffusion creep (Supplementary Fig. S3c). Laboratory and geodynamics constraints suggest ε_ T ~10–15 s–1, in agreement with seismic evidence of a mantle predominantly deforming under diffusion creep72.

Plasticity is implemented in the models using a yielding criterion:

τ y = C0+μp(1–λ), at τ>τ y (6)

where C 0 is the cohesion and µ the friction at which yielding set forth in the viscous material and λ is the pore pressure fraction. To fit the common assumption of near-lithostatic pore pressure, λ is set to 1, and the yield stress reduces to:

τ y = C 0 (7)

The models assimilate the yielding criterion through a composite viscoplastic effective viscosity:

η = min (η creep , τ y / 2ε_ II) (8)

where ε_ II is the second invariant of strain rate tensor.

We have simulated the geodynamics of the Southeast Asian subduction zone in a part of a spherical shell. The model domain is discretised using a uniformly spaced mesh in the longitude and latitude direction and a non-uniformly spaced mesh in the radial direction. As a result, the element size is 15×15×10 km in the upper mantle and 15×15×37 km in the lower mantle in longitude, latitude, and radial direction. We embed in each element 20 Lagrangian particles, capitalizing on the code’s ability to track material properties and increase the resolution, allowing a more accurate recovery of different materials stress. All boundaries are set to free slip conditions. For full details on each model setup, refer to Supplementary Tables S1 and S2.

We treat the Indo-Australian subducting plate as a single plate. We rotate the boundary box of the numerical model by 16° clockwise to avoid the implementation of viscosity weakness along the Owen Fracture Zone (OFZ). The rotation angle is calculated by optimizing the angle between the boundary box and velocities of the Indo-Australian plate in the NNR reference frame (Supplementary Table S3). The slab geometry in the upper mantle, that is, the subducted lithosphere, is reconstructed using slab depth and thickness provided by the Slab2 dataset27. The thickness h of the oceanic plate, that is, the unsubducted lithosphere, is recovered using the half-space cooling model, defined by the 1300°C isotherm, and oceanic crustal age dataset41 as h = 2.32 sqrt(κt), where κ=10–6 m2 s–1 and age t in s, for t ≤ 100 Ma, and h = 130 km for older lithosphere. The lithospheric thickness is chosen to match the thickness from the dataset Slab2. This includes thinning around the Capricorn zone, Wharton Basin, however we do not address the putative deformation zone between Indian, Capricorn and Australian plates, to avoid the implementation of ill-posed viscosity weaknesses. The density contrast with the underlying mantle is 50 kg m–3 (refs. 30,73). The continental lithosphere is thinned to 90 km, implying a buoyancy increase of ~30%, while in Model 3c we additionally vary the density of the Indian slab (see Supplementary Table S2). The upper plate is modelled as a uniform lithosphere, with no internal boundaries and constant thickness of 90 km. The lower mantle density anomalies are scaled from the tomographic models using scaling factor of 0.3 and 0.6 for the S-wave and P-wave tomographies to convert velocity anomalies into density anomalies, respectively, while in the upper mantle we ignore density anomalies from tomographies due to their poor resolution, following ref. 45. The viscosity in the mantle varies only in the radial direction, we set a viscosity jump of a factor 30 between the upper and lower mantle.

The subducting lithosphere has a viscoplastic rheology in the top 30 km and adopt a linear viscous flow law in the rest of the slab. In the models, the brittle ductile transition is defined self-consistently, i.e., occurs at shallower depths, according to the coupling stress. The choice of the maximum viscoplastic layer is based on the observed maximum depth of the brittle to ductile transition in an oceanic plate74. The viscoplastic rheology decouples effectively the subducting and upper plate and has been tested against stress at subduction zones, tectonics and earthquake simulations34,35,40,75. For simplicity, the upper plate is treated as a uniform-thickness viscous plate to highlight the role of boundary stresses, as opposed to those resulting from geometrical singularities arising from upper plate structures.

Tests of viscosity, cohesion and numerical resolution, among others, are performed in a range of different models in Cartesian and spherical coordinates systems and two- and three-dimensional spaces, to reduce the computational cost and are publicly available76.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "In this paper, a three-dimensional viscoplastic model is constructed for the Andaman-Sumatra-Java subduction zone, and\nthe spatial changes in strain and stress fields associated with plate subduction are quantitatively calculated, taking into\naccount factors such as the balance of forces acting on the slab due to age changes in the direction parallel to the trench.\nFurthermore, the paper compares seismic activity, including megathrust earthquakes that have occurred in this region, with\nthe calculated effective friction coefficient, showing that there is a good correspondence between the spatial changes in\nseismic activity, making this a very interesting paper. However, there are concerns regarding the following points, which\nshould be properly examined.\nMajor comments\n1.In sections of “Introduction” and “Seismotectonics of the Sunda subduction zone”, nothing concrete is written, and it is\nunclear what and how much the authors want to clarify. For example, it is unclear how long a period is, what kind of data\nspecifically refers to geodetic data, and what kind of rheology model it is. In addition, similar technical terms are used without\ndefining the differences, and it is difficult to grasp their exact meaning. These points should be clearly written. Please refer to\nthe comments below for specific points.\n2. L.585 The data and period used should be written. In particular, the results will vary greatly depending on whether the\nperiod includes the megathrust earthquakes of 2004, 2005, and 2007. If the period includes these earthquakes or the period\nafter the earthquake, deformation due to afterslip and/or viscoelastic response should also be taken into account. The\ndilatation also varies depending on the density of observation points. For example, see the DDC by Shen et al. (1996, JGR) .\nThe data used, including the error, should be properly presented.\n3. In relation to the point 2 described above, the model in this study assumes viscoplasticity and shows the strain field in a\nsteady state as a long-term average. Therefore, as in the above point 2, it should be impossible to directly compare the strain\nfield obtained by the model with the strain field obtained by a specific short-term earthquake cycle, but there is no description\nof this. More specifically, in L. 246 -262, the results obtained from the model and the GNSS observation data are completely\ndifferent in terms of the time scale considered, so they should not be directly comparable. In particular, in areas such as\nIndonesia where megathrust earthquakes occur repeatedly, afterslip and/or viscoelastic response after a megathrust\nearthquake lasts for a long period of time, and it is difficult to understand the significance of comparing it with this model,\nwhich does not take into account the occurrence and healing of earthquakes.\n4. Is the viscosity η of L .490 a constant? Viscosity usually depends on temperature, strain rate, and creep law. So, if it is a\nconstant, it is an unrealistic assumption, and this assumption, including the value itself, should have a large effect on the\nresults. It is well known that the convection pattern of the mantle wedge is significantly different when the viscosity is\nconstant and when it is temperature dependent, and this effect should not be ignored in this study.\n5. The model in this study is a viscoplastic model, and the resulting strain and stress are originated from the rheology. On the\nother hand, earthquake phenomena are caused by the release of elastic strain, as is well known. Therefore, the two cannot\nbe compared. The former is described as an interseismic period in the text (L.404-405), but this understanding is incorrect.\nThe effective friction coefficient obtained in this study's model is the average effective friction coefficient after going through\nmany earthquake cycles, and there is no guarantee that it will match the value of the effective friction coefficient in a specific\ninterseismic period of one earthquake cycle. For this reason, it is incorrect to compare and discuss this model, which does\nnot take elasticity into account, with earthquake phenomena. If the authors want to compare it with earthquake phenomena,\nthey should construct a model that also takes elasticity into account and then discuss it.\nMinor comments\nL.70 The model by Wada and Wang (2009) assumes that the depth of coupling-decoupling on the upper surface of the slab\nis 70-80 km, and does not necessarily explain the heat flow data. So, the temperature structure model is questionable. It is\nrecommended to cite more recent reliable papers.\nL.109- L.110 What do seismic coupling and tectonic coupling mean? They should be clearly defined and used.\nL.109 constrained geodetically It should be clearly stated what kind of geodetic data it refers to and for what period.\nL.114 How long is long-term? It should be expressed in numbers.\nL.118-119 negligible further south This does not seem to be the case from the red arrow in Fig. 1a. The explanation in the\ntext should be written accurately based on the figure.\nL.119 tectonics and geodetic deformation The difference between the two should be clearly stated.\nL.148 What does “stress coupling” mean? It should be clearly defined and used.\nL .148 What does \" full force balance\" mean? You should clearly indicate what kind of force balance you have in mind\nusing an equation.\nL.149 The type of model (what is the input and what is the output) should be clearly described, including whether it is\nkinematic or dynamic .\nL.161 Lower mantle contribution to what?\nL.161 In Slab2 model, the thickness of the slab is given. How is the thickness of the slab in the lower mantle continuity or\ndiscontinuity based on seismic tomography?\nL.168 How and why?\nL.171-189 Even after reading this section, I cannot imagine the differences between Model 1, Model 2, and Model 3. An\neffort should be made to help the reader understand by using a table or drawing schematic figures.\nL.198 the stress What kind of stress is this ? There are six independent components. Which component do you mean ? Also,\ndoes the stress originated from elastic and/or viscous and/or plastic ? These points should be clearly stated.\nL .202 six models described above\nAren't there five models, Model 1, 2, 3a, 3b, and 3c ? Do you also include Model 3? In the latter case, the difference between\nModel 3 and Model 3a, 3 b and 3 c should be clearly stated. In this case, the naming of the Model is confusing. So, It should\nbe Model 1-Model 6 or Model 1, 2, 3, 4a, 4b, 4c .\nL.212 towards the, where This passage is not grammatically correct. Rewrite it in correct English.\nL.213, L. 214 toroidal, radial Define their directions before use them.\nL.214 stress and tractions have six and three components, respectively, but which component and which direction are you\nfocusing on? Please write it so that it is clear.\nL.215 What physical quantity is used to determine the “increasing the coupling” ? Please write it in a way that is easy to\nunderstand.\nL. 216 additional forces What kind of forces are added to what? Write it so that it is clear.\nL. 219 the dilatation strain rate pattern Which part of model are you referring to? Please clarify.\nL 280 The stresses Same as the comment for L. 214.\nL .307 seismic plate boundary This is not common terminology. Define the word.\nL .327 forearc or fore arc ? Consistently use the same expression throughout the paper.\nL.427-430 It is premature to draw such a conclusion based on earthquake data from only a limited period of time. In fact, the\nconsidered to be weak, for the first time in 1,000 years. The plate ages off the Java coast and off the Tohoku coast are about\nthe same, and their tectonic settings are also similar.\nL.479 Please clearly describe the elements of the forces acting on the plates considered in this study (ridge push, slab pull,\netc.). Also, how does the equation that their sum is 0 relate to the equation of motion f ? Please clearly show it in the form of\nan equation.\nL.515-518 In Slab2, the thickness of the slab is given, as written in L514-515. I don't understand the relationship between\nthat value and the value obtained by the formula written here. Please write it so that the readers can understand.\nL.521 What does “linearly viscous” mean?\nL.613 The upper end of Fig. 4a should be called \"contraction\" rather than \"compression\" .\nL.630 Is the “distance” in the vertical axis in Fig. 6b measured along the trench?\nL.630, L.631 Specify the period of data.\nL.631 In the Sumatra-Andaman region, there have been megathrust earthquakes of Mw 8.8-9.2 in 1833, Mw 8.3-8.5 in 1861,\nMw 7.9 in 1881, and Ms 8.1 in 1941. It is premature to make comparisons based on only short-term data.\nL.637-639 In Fig. 6e, the maximum value of apparent friction is 0.85, but in Fig. 6g, the maximum value is 2.0. Why are they\ndifferent? Shouldn't they both be on the same scale?", "reviewer_2": "I agreed to do this review even though I'm not a specialist in geodynamic modeling. My comments will therefore be very\ngeneral and will need to be supplemented by those of a numerical modeler more familiar with the methodology used.\nThe context of subduction in Indonesia is used here to demonstrate that high-performance 3D numerical computation\nincorporating long-term first-order driving and resisting forces enable the authors to account for short-term stress and\nseismogenic behavior in the vicinity of the subduction interface.\nOverall, I agree that the model performs well, despite the simplifications inherent in this type of approach.\nI am wondering whether the estimate of trench motions within an NNR reference frame rather than for example spreading\nalignment (Becker et al., 2015) to compare with the model is justified. Table S3 gives the values for trench advance/retreat\nwithin 3 reference frames that are quite old. Why not using improved reference frames like NNR or those of Becker et al. or\nGMHRF-1Ma, or T25M within MORVEL56 rather than NUVEL1A (which doesn't include the Sunda plate) plateset ?\nUsing NNR-MORVEL56 (Argus et al., 2011) which includes the Sunda plate, I obtain trench advance values off Sumatra\n(-10 to -20 mm/yr), almost stationary trench off Java (-1 to -4 mm/yr), and trench retreat values (10 to 19 mm/yr) off Andaman.\nThis is still consistent with the description of the authors (p.4).\nThe authors say that the largest compression is located along the Sumatra segment (p.4 & p.8). I agree that several recent\nstudies (e.g., Widiyantoro et al., 2024) show evidences of compression there but they are still tiny for me with respect to the\ncompression recorded in the backarc of Java (Koulali et al., 2016, 2017).\nI may have missed a step, but in the boundary conditions imposed by the various models, I didn't see any mention of\nAustralia's continental subduction under the Sunda arc. According to me, the first order setting in this area is an arcuate\nsubduction (varying slab pull along strike ...) between two collision zones, a major and mature one in the Himalaya region\nand another younger one north of Australia. Could the authors comment on this ?\nAlso, what causes the diverging toroidal mantle flow in front of Sumatra slab ? Is it a consequence of the trench advance ? I\nhad in mind another toroidal flow from below the Tibet to the Sunda area (e.g., Flower et al., 1998 ...) but here the mantle\nflows from Thailand/Burma to the Tibet (fig.1d) which puzzles me ?\np7, line 29 : what does that mean : \"large compression in Sumatra extends as far as the Indonesian region\" ? Could the\nauthors be more precise ?\np.8 - lines 54-58 are a bit confusing because I am not sure to see backarc extension in all models except 3b in the Andaman\nback-arc region (pale blue in models 1,2, 3a, 3b)\np.12 : did the authors sort the seismic events (interface only) in their plot (fig.6) ?\np.12 - lines 95-96 : I don't understand this sentence. Could the authors better explain ?\nfew typo errors:\n- p.3 - line 91 : remove \"the\"\n- p.4 - line 8 : remove e in \"enucleated\"\n- p.7 - line 12 : toward the ..???\n- p.10 - line 25 : northern\np.13 - line 37 : these show ...???\np.14 - line 45 : megathrust\np.14 - line 63 : match\nThe text would be easier to read if it were a little more structured, with subheadings for example like \"Limitations\" p.14 -\nparagraph 46-56, \"conclusion\" p.14 - paragraph 57-71.\nThe fig.6 a,b,d,e are difficult to read, especially the rupture lengths and envelops whose color contrast makes it difficult to\nsee them properly.\nTo conclude, this study deserves to be published after minor revisions.\nSerge Lallemand", "reviewer_3": "This cross-disciplinary study demonstrates that stresses along convergent margin primarily are controlled by the down-going\nslab’s buoyancy and the subduction-induced flow in the upper mantle and are as such influenced by dynamics operating\nthroughout a 3-D space encompassing plates and mantle. The study nicely encompasses three model families to test\ndifferent aspects systematically and individually and relate them to various observables of the Sunda margin region to\nprovide these insights. The model has been tested in-depth, although detailed tests (such as a resolution test) are only\naccessible in a previous (but open-access) publication. The underlying open-source code, Underworld, is state-of-the-art\nand well tested. Finally, the study is nicely and clearly written, and the nice figures are helpful, accurate and accessible.\nThis study demonstrates that models like the ones presented here can, and here do, close the gap between subduction\ndynamics and the long-term (multi-earthquake cycle) strength of megathrusts. As such, I see it as a nice fit for publication in\nNature Communications. I only have some minor comments outlined below.\nSpeaking of inter-disciplinarity: I think the presentation of this study, as intended for an interdisciplinary journal like Nature\nCommunications, should get a little more attention to the text.\nAs a reviewer, I only comment on the scientifically relevant bits, as for example on what is “deep subduction”: For me this is\nsubduction at around the upper mantle transition zone (670 km depth). For others this is maybe everything below the\nsurface. And it is not made clear what this study’s “deep subduction” refers to here. Could this be clarified.\nModel setup:\nThe models apparently employ a free-slip top boundary (according to line 512)—which obviously caught my eye—meaning\nthat there is no surface topography, which causes the subduction zone to behave significantly different to a model with a free\nsurface (e.g., the natural state of such models is double sided subduction, rather than a realistic single-sided sinking – see\nCrameri et al., 2012) because the subducting plate is forced to stay flat and at the top boundary until the very location of the\nplate interface. Has this simplification been tested on how it affects the stress distribution at the plate interface? If not, could\nthis be discussed shortly for a fair representation of the model?\nCrameri, F., P. J. Tackley, I. Meilick, T. V. Gerya, and B. J. P. Kaus (2012), A free plate surface and weak oceanic crust\nproduce single-sided subduction on Earth, Geophys. Res. Lett., 39(3), L03306, doi:10.1029/2011GL050046\nLines 54-55: I don’t understand what the authors want to communicate with this sentence. e.g., what does “set the stage”\nmean here? Should it say: “integrating [stresses] across these scales”? – could you clarify this?\nLines 55-58: Same here, not sure that comes across as intended. “the loading of megathrusts by subducting plate motions—\nwhich increases plate interface friction” – The “motion” (also, why plural?) would rather reduce plate interface friction, from\nmy point of view (there is more friction if an interface is static). Or you write: “the loading [of stress] faces scrutiny from recent\nhigh-energy earthquakes”. Instead it should say something like: “Our current understanding about xxx faces scrutiny from\nrecent insights gained from xxx”.\nI think I get what you want to say here, but it is written too confusingly, so you might want to consider rewriting.\nLine 78: What does “large-scale” mean? – “Crust-scale”, “Lithosphere-scale”? “Upper-mantle scale”?\nLine 148: What does “deep” refer to here?\nLine 197: What stability? Numerical or otherwise?\nLine 198 and 376: “Kyrs” – Units don’t have a plural form. Use “kyr” for a scientific notion and to be in line with other\ninstances in the manuscript.\nLine 445: Typo in “Megathrust”\nLines 445-456: I appreciate this modelling/study limitation paragraph. I think that is very useful for any reader to understand\nwhat has been done and achieved in this study. Of course this could be extended, but I think works here as is.\nFigures:\nThe Scientific colour maps (Crameri 2018) are nicely put in use here. – That is great to see! Would you mind also citing\nthem, as they are under a MIT license.\nCrameri, F. (2018). Scientific colour maps. Zenodo. https://doi.org/10.5281/zenodo.1243862\nIf the authors want to make their figures perfect (here on in future manuscripts), avoiding red&green combinations at similar\nlightness values (as here with lines) could additionally be avoided. But this is a pro-level-suggestion, in particular in\ncomparison to other current publications.\nThank you for a very interesting study I enjoyed reading!\nFabio Crameri\n20.11.2024"}, "author_response": "We are grateful for the support and the improvement suggested. The level of detail provided in this\nreview will be of great benefit to this paper.\nMajor comments\n1.In sections of “Introduction” and “Seismotectonics of the Sunda subduction zone”,\nnothing concrete is written, and it is unclear what and how much the authors want to\nclarify. For example, it is unclear how long a period is, what kind of data specifically refers\nto geodetic data, and what kind of rheology model it is. In addition, similar technical terms\nare used without defining the differences, and it is difficult to grasp their exact meaning.\nThese points should be clearly written. Please refer to the comments below for specific\npoints.\nWe have rewritten the introduction to clarify the aims and the terminology. We have added more\ndetails on data and rheology used throughout the text and better explained our choices.\n2. L.585 The data and period used should be written. In particular, the results will vary\ngreatly depending on whether the period includes the megathrust earthquakes of 2004,\n2005, and 2007. If the period includes these earthquakes or the period after the\nearthquake, deformation due to afterslip and/or viscoelastic response should also be taken\ninto account. The dilatation also varies depending on the density of observation points. For\nexample, see the DDC by Shen et al. (1996, JGR) . The data used, including the error,\nshould be properly presented.\nFirst, we share the reviewer’s concern, however, this comment focuses on the speculative part of\nour work, therefore, we do not consider this to be a potential flaw. In fact, the dataset questioned by\nthe reviewer is only one of those used to validate the modelling proposed. We have made sure that\nthis is clearly understood as a speculative part of our work.\nThen, we use the dataset by Kreemer et alii, this is the only global dataset publicly available which\nis a compilation of published and available data. The Global Strain Rate Map avoids entirely the\nperiod around and stations close to the major earthquakes to minimise the effect of transients, as the\nreviewer points out. We also are aware of the shortcoming of the methodology used, which\nnevertheless is the standard methodology in this field, also addressed in the Kreemer et alii’s\npublications.\nThe impact of the use of this dataset is mitigated in our work by using other observables. Therefore,\nwe have extensively compared this dataset and other datasets across timescales, from the geodetic\nto the geological observations, in our previous study (also see next comment):\nZuhair, M., Gollapalli, T., Capitanio, F.A., Betts, P.G. and Graciosa, J.C., 2022. The role of slab\nsteps on tectonic loading along subduction zones: Inferences on the seismotectonics of the Sunda\nconvergent margin. Tectonics, 41(9), p.e2022TC007242.\nWe have clarified these aspects in our work and added a “limitations” section to clear this even\nfurther.\n3. In relation to the point 2 described above, the model in this study assumes\nviscoplasticity and shows the strain field in a steady state as a long-term average.\nTherefore, as in the above point 2, it should be impossible to directly compare the strain\nfield obtained by the model with the strain field obtained by a specific short-term\nearthquake cycle, but there is no description of this. More specifically, in L. 246 -262, the\nresults obtained from the model and the GNSS observation data are completely different in\nterms of the time scale considered, so they should not be directly comparable. In\nparticular, in areas such as Indonesia where megathrust earthquakes occur repeatedly,\nafterslip and/or viscoelastic response after a megathrust earthquake lasts for a long period\nof time, and it is difficult to understand the significance of comparing it with this model,\nwhich does not take into account the occurrence and healing of earthquakes.\nThis is a very interesting point, the very reason why we deem our work useful.\nThis comment goes back to the previous one. We agree with the reviewer, for this reason we have\nused the dataset by Kreemer et alii. This dataset explicitly excludes seismic strain (i.e., only GPS)\nand exclude stations and periods close to large earthquakes, to avoid the transient component due to\nearthquakes.\nWhile we agree this methodology may need to be improved, it is a significant dataset that cannot be\nsimply dismissed.\nThe shortcomings the reviewer points to are well known and addressed in the original work. The\nGlobal Strain Rate Map is a project supported by several consortia (UNAVCO and IUGG’s ILP)\nand has gone through several iterations, through the years. To the best of our knowledge this is the\nmost robust dataset available.\nAdditionally, we have chosen the Sunda subduction zone as it is likely an area where the short-term\n(interseismic) deformation is in agreement with the long-term, offering the unique opportunity to\nunderstand more how these two scales are linked.\nWe have amended the text to better convey these points.\n4. Is the viscosity η of L .490 a constant? Viscosity usually depends on temperature, strain\nrate, and creep law. So, if it is a constant, it is an unrealistic assumption, and this\nassumption, including the value itself, should have a large effect on the results. It is well\nknown that the convection pattern of the mantle wedge is significantly different when the\nviscosity is constant and when it is temperature dependent, and this effect should not be\nignored in this study.\nViscosity is not constant, but depends on stress non-linearly, i.e., plastic, in the crust, while the rest\nof the model space uses layered piecewise constant viscosities to reproduce the layering resulting\nfrom the temperature and pressure dependence of rocks. Additionally, because the model solution is\n“instantaneous” the effect of temperature diffusion is neglected. This is a common procedure in this\ntype of modelling, e.g., Faccenna et al., 2010; Straume et al., 2024.\nTo further provide convincing evidence, we have added in the supplementary information a test of\nthe models using different mantle rheologies, from Newtonian to non-Newtonian. We show that the\nquantities discussed, that is, the velocity, the second invariant of the stress and strain rate around the\nmargin show very similar patterns and, although the magnitudes are different, they do not vary\nmore than a factor ~2 in the velocity and the stress and up to a factor 5 in the strain rates.\nThis is an expected consequence of the choice of a dynamic system, where stresses are determined\nby mass and gravity, primarily, and the scale of the subduction system, dominated by Newtonian\nviscous drag, while non-linearities are confined to small regions.\nFaccenna, C. and Becker, T.W., 2010. Shaping mobile belts by small-scale\nconvection. Nature, 465(7298), pp.602-605.\nStraume, E.O., Steinberger, B., Becker, T.W. and Faccenna, C., 2024. Impact of mantle convection\nand dynamic topography on the Cenozoic paleogeography of Central Eurasia and the West Siberian\nSeaway. Earth and Planetary Science Letters, 630, p.118615.\nTo reply to the following, we prefer to break it down, addressing each point individually:\n5. The model in this study is a viscoplastic model, and the resulting strain and stress are\noriginated from the rheology.\nSee comment above\nOn the other hand, earthquake phenomena are caused by the release of elastic strain, as\nis well known. Therefore, the two cannot be compared.\nThis point is important as it emphasises the reason for this work.\nFirst, we have cleared that the comparison is speculative, indeed, although this does not affect the\noutcomes of our work.\nSecond, the comparison is nevertheless possible, in principle: it is acknowledged that the energetic\nbudget of tectonic processes implies that the long-term tectonic loading results into a stored energy,\nthen released during earthquakes. In this sense, elasticity acts as a buffer, storing stress loaded by\nsubduction, then released by plasticity.\nHere, we do not model the earthquake rupture mechanism, but only one term, the stress loaded by\nsubduction, then we compare it with the measured stress released by earthquake. Modelling\nelasticity and rupture mechanisms would likely lead to the same result, buffering between the\nloaded stresses and those released by seismic events, and we’d compare the magnitude of our\nsynthetic earthquakes with the observed.\nThen, the comparisons of the computed stresses and magnitudes is a viable test of the assumption\nabove with two possible outcomes: it either confutes or corroborates the postulated equivalence\nbetween loading stress and released stress.\nThe former outcome, that is, the (elastic) stresses released during earthquakes surpass the calculated\n(viscous) loading stresses, would confute our model. This is not the case.\nThe latter, the (elastic) stresses released during earthquakes are similar to, or lower than, the\ncalculated (viscous) loading stresses, does not prove our model correct, yet it does not confute it,\neither.\nMost importantly, there is an additional reason not to model elasticity and earthquakes rupture\nmechanisms.\nWe only have empirical laws for the description of stress release during earthquakes. The use of\nthese adds a degree of freedom to the model and leads to nonunique solutions. This is illustrated by\nthe implementations of physical models of rupture, e.g., Ulrich et al., 2022, Kammer et al., 2024.\nWhile these works provide excellent insights, they require an a priori definition of the stress acting\non faults and energy, respectively. This latter is what we have addressed in our work and is, in this\nsense, complementary.\nTherefore, although caution is in order, as highlighted by the reviewer, we believe that we have\nprovided several reasons (also, see below) to illustrate how the case studied here may provide\ninsights, and engaging in this thought exercise - or even considering it as a provocation - remains\nworthwhile.\nWe have clarified this further in the revised manuscript.\nKammer, D.S., McLaskey, G.C., Abercrombie, R.E., Ampuero, J.P., Cattania, C., Cocco, M., Dal\nZilio, L., Dresen, G., Gabriel, A.A., Ke, C.Y. and Marone, C., 2024. Earthquake energy dissipation\nin a fracture mechanics framework. Nature communications, 15(1), p.4736.\nUlrich, T., Gabriel, A.A. and Madden, E.H., 2022. Stress, rigidity and sediment strength control\nmegathrust earthquake and tsunami dynamics. Nature Geoscience, 15(1), pp.67-73.\nThe former is described as an interseismic period in the text (L.404-405), but this\nunderstanding is incorrect. The effective friction coefficient obtained in this study's model is\nthe average effective friction coefficient after going through many earthquake cycles, and\nthere is no guarantee that it will match the value of the effective friction coefficient in a\nspecific interseismic period of one earthquake cycle.\nIndeed, this is the speculative part of our work, and we have clarified further:\n“The averaged apparent friction is not necessarily the local frictional fault strength. This latter\nvaries over the duration of the single seismic event due to shear stress release, while the normal\nstress is unlikely to change. Considering the modest stress drop during earthquakes69, the difference\nbetween the averaged and time-dependent friction remains minor.”\nFor this reason, it is incorrect to compare and discuss this model, which does not take\nelasticity into account, with earthquake phenomena. If the authors want to compare it with\nearthquake phenomena, they should construct a model that also takes elasticity into\naccount and then discuss it.\nWe do not aim to model earthquakes phenomena, here, but only address the long-term stresses,\nthat is, those available to drive deformation at the margin, and if and how these are compatible with\nthose released in the seismotectonics of the area.\nThe elastic stresses are due to the loading during subduction and are “stored”, that is accounted for\nas stored energy (W = ∫ dV). Through time, the energy released sums up to the total energetic\nbudget of subduction, as discussed above. Our comparisons illustrate this point, testing whether the\nstresses loaded by subduction and those released by seismicity have any correlation.\nIndeed, for the reasons above, implementing elasticity remains strictly necessary to the modelling of\nearthquake mechanisms, although it is not the goal of our work. On the other hand, as highlighted\nabove, modelling elasticity and rupture mechanics may lead to nonunique solutions.\nMinor comments\nL.70 The model by Wada and Wang (2009) assumes that the depth of coupling-decoupling\non the upper surface of the slab is 70-80 km, and does not necessarily explain the heat\nflow data. So, the temperature structure model is questionable. It is recommended to cite\nmore recent reliable papers.\nThanks for highlighting this. We have removed this reference as it was redundant and not used\nelsewhere in the text.\nL.109- L.110 What do seismic coupling and tectonic coupling mean? They should be\nclearly defined and used.\nAmended and added the definition of seismic and tectonic coupling.\nL.109 constrained geodetically It should be clearly stated what kind of geodetic data it\nrefers to and for what period.\nAmended and moved the reference.\nL.114 How long is long-term? It should be expressed in numbers.\nRemoved, it was stated the sentence after.\nL.118-119 negligible further south This does not seem to be the case from the red arrow in\nFig. 1a. The explanation in the text should be written accurately based on the figure.\nAmended, that was a typo, apologies.\nL.119 tectonics and geodetic deformation The difference between the two should be\nclearly stated.\nAmended.\nL.148 What does “stress coupling” mean? It should be clearly defined and used.\nWe have defined this throughout the text. In this version, we have improved the description of the\nterminology used.\nL .148 What does \" full force balance\" mean? You should clearly indicate what kind of\nforce balance you have in mind using an equation.\nWe have rewritten the conservation equations to show how the stress balance and amended the text.\nAdditionally, we have added a sketch in the Supplementary Figure 1 to explain the different forces.\nWe have added references to other works, as the subduction force balance is now considered\nacknowledged we did not modify the main manuscript further.\nL.149 The type of model (what is the input and what is the output) should be clearly\ndescribed, including whether it is kinematic or dynamic .\nAdded.\nL.161 Lower mantle contribution to what?\nAmended, we meant the lower mantle mass.\nL.161 In Slab2 model, the thickness of the slab is given. How is the thickness of the slab in\nthe lower mantle continuity or discontinuity based on seismic tomography?\nWe have further cleared that the Slab 2 constrains the slab mass in the upper mantle only, as it\ncombines tomography and Benioff zone seismicity.\nL.168 How and why?\nThis was not too clear, agreed. We do not rotate the dataset, we rotate the numerical model walls to\nbe parallel to the convergence along Java.\nWe have amended the text.\nL.171-189 Even after reading this section, I cannot imagine the differences between Model\n1, Model 2, and Model 3. An effort should be made to help the reader understand by using\na table or drawing schematic figures.\nApologies. We have amended figure 1 to include the boundaries the different models and improved\nthe description in the manuscript.\nL.198 the stress What kind of stress is this ? There are six independent components.\nWhich component do you mean ? Also, does the stress originated from elastic and/or\nviscous and/or plastic ? These points should be clearly stated.\nWe have amended the text to clarify this. However, here we simply meant that we show the results\nfrom a single time step and amended the text accordingly.\nL .202 six models described above. Aren't there five models, Model 1, 2, 3a, 3b, and 3c ?\nDo you also include Model 3? In the latter case, the difference between Model 3 and\nModel 3a, 3 b and 3 c should be clearly stated. In this case, the naming of the Model is\nconfusing. So, It should be Model 1-Model 6 or Model 1, 2, 3, 4a, 4b, 4c .\nWe have clarified the reasoning. Models 1 to 3 reflect different upper mantle setups, models with\na,b,c modify the numeral model’s lower mantle.\nL.212 towards the, where This passage is not grammatically correct. Rewrite it in correct\nEnglish.\nAmended\nL.213, L. 214 toroidal, radial Define their directions before use them.\nWe have rewritten this sentence trying to avoid jargon:\n“Instead, further north, the mantle is complicated by the varying slab depth, and the mantle flow is\nforced around the northern shorter slab, upwelling from the sub-slab to the mantle wedge region\nbeneath central Sumatra, where the radial velocity is largest, that is, the direction parallel to the\nradius of the spherical model.”\nL.214 stress and tractions have six and three components, respectively, but which\ncomponent and which direction are you focusing on? Please write it so that it is clear.\nWe have amended the sentence: “This configuration results in excess pull due to the neighboring\nJava slab, extending to the mantle transition zone, and additional trench-ward tractions beneath the\nSumatran upper plate segment, increasing the in-plane stresses between the plates, i.e., the tectonic\ncoupling.”\nL.215 What physical quantity is used to determine the “increasing the coupling” ? Please\nwrite it in a way that is easy to understand.\nSee above\nL. 216 additional forces What kind of forces are added to what? Write it so that it is clear.\nAdded.\nL. 219 the dilatation strain rate pattern Which part of model are you referring to? Please\nclarify.\nAmended.\nL 280 The stresses Same as the comment for L. 214.\nThis is a sentence that introduces the concept that the stress and the strain rate tensors are\ncorrelated. Because the viscosity is isotropic, a Newtonian equivalence (linear) between the two can\nbe used, which extends to each component and their combinations.\nWe have specified in the following paragraph when we refer to a specific component, in this case,\nthe second invariant.\nL .307 seismic plate boundary This is not common terminology. Define the word.\nAmended.\nL .327 forearc or fore arc ? Consistently use the same expression throughout the paper.\nAmended.\nL.427-430 It is premature to draw such a conclusion based on earthquake data from only a\nlimited period of time. In fact, the 2011 Tohoku-oki earthquake with a magnitude of Mw 9.0\noccurred off the coast of Tohoku, where interplate coupling was considered to be weak, for\nthe first time in 1,000 years. The plate ages off the Java coast and off the Tohoku coast\nare about the same, and their tectonic settings are also similar.\nWe agree. Drawing conclusions on what causes earthquakes is out of the goal of this work, here we\nonly test whether the computed stresses are compatible with what released seismically.\nWe have now expanded the Limitations section to further clarify this aspect.\nL.479 Please clearly describe the elements of the forces acting on the plates considered in\nthis study (ridge push, slab pull, etc.). Also, how does the equation that their sum is 0\nrelate to the equation of motion f ? Please clearly show it in the form of an equation.\nWe have added a sketch and a description of the force balance.\nAlso, we have rewritten the equation solved to show how the momentum conservation equation\nsolves for the stress balance in the volume, due to static (pressure), viscous stresses due to motion\nand dynamics stresses, due to mass and gravity. We hope this address the reviewer’s concern.\nL.515-518 In Slab2, the thickness of the slab is given, as written in L514-515. I don't\nunderstand the relationship between that value and the value obtained by the formula\nwritten here. Please write it so that the readers can understand.\nWe have cleared this further and hope this is clearer, now. The slab is constrained by the Slab 2.0\nmodel, the plate is constrained by the crust age, two parts of the lithosphere.\nL.521 What does “linearly viscous” mean?\namended\nL.613 The upper end of Fig. 4a should be called \"contraction\" rather than \"compression\".\nApologies, amended.\nL.630 Is the “distance” in the vertical axis in Fig. 6b measured along the trench?\nYes, we have amended the figure.\nL.630, L.631 Specify the period of data.\nDone\nL.631 In the Sumatra-Andaman region, there have been megathrust earthquakes of Mw\n8.8-9.2 in 1833, Mw 8.3-8.5 in 1861, Mw 7.9 in 1881, and Ms 8.1 in 1941. It is premature\nto make comparisons based on only short-term data.\nAs explained before, the dataset used, although short, excludes the strain associated with large\nearthquakes, providing the interseismic strain available, ideally representing the tectonic steady\ndeformation, comparable with our models. We have better cleared this in the main text, and we\nwelcome suggestions.\nL.637-639 In Fig. 6e, the maximum value of apparent friction is 0.85, but in Fig. 6g, the\nmaximum value is 2.0. Why are they different? Shouldn't they both be on the same scale?\nThese two figures use the same data, but we use them to show something different.\nThe saturation of the colour scale 6d, e, is necessary to better visualise the data. The threshold value\nused is under the assumption that any value >0.85 indicates conditions for rupture (above the yield\nstrength envelope).\nWe have amended the text to make this clearer.\nReviewer #2 (Remarks to the Author):\nI agreed to do this review even though I'm not a specialist in geodynamic modeling. My\ncomments will therefore be very general and will need to be supplemented by those of a\nnumerical modeler more familiar with the methodology used.\nThe context of subduction in Indonesia is used here to demonstrate that high-performance\n3D numerical computation incorporating long-term first-order driving and resisting forces\nenable the authors to account for short-term stress and seismogenic behavior in the\nvicinity of the subduction interface.\nOverall, I agree that the model performs well, despite the simplifications inherent in this\ntype of approach.\nI am wondering whether the estimate of trench motions within an NNR reference frame\nrather than for example spreading alignment (Becker et al., 2015) to compare with the\nmodel is justified. Table S3 gives the values for trench advance/retreat within 3 reference\nframes that are quite old. Why not using improved reference frames like NNR or those of\nBecker et al. or GMHRF-1Ma, or T25M within MORVEL56 rather than NUVEL1A (which\ndoesn't include the Sunda plate) plateset ?\nUsing NNR-MORVEL56 (Argus et al., 2011) which includes the Sunda plate, I obtain\ntrench advance values off Sumatra (-10 to -20 mm/yr), almost stationary trench off Java (-\n1 to -4 mm/yr), and trench retreat values (10 to 19 mm/yr) off Andaman. This is still\nconsistent with the description of the authors (p.4).\nThis is a very good point. We have added the Becker et al. 2015’s dataset, considering this to be the\nmost recent and redrafted figure 4.\nThe authors say that the largest compression is located along the Sumatra segment (p.4 &\np.8). I agree that several recent studies (e.g., Widiyantoro et al., 2024) show evidences of\ncompression there but they are still tiny for me with respect to the compression recorded in\nthe backarc of Java (Koulali et al., 2016, 2017).\nIn the submitted version we do refer to compression in these areas, in fact, in lines 358-360:\n“The horizontal stress invariant at the surface,  , in all models shows high values of ~150 MPa in\nII\nthe Java and southern Sumatra forearc and decrease further north.”\nHowever, we have now improved the discussion of this key point raised by the reviewer and added\na reference to the work mentioned.\nI may have missed a step, but in the boundary conditions imposed by the various models, I\ndidn't see any mention of Australia's continental subduction under the Sunda arc.\nAccording to me, the first order setting in this area is an arcuate subduction (varying slab\npull along strike ...) between two collision zones, a major and mature one in the Himalaya\nregion and another younger one north of Australia. Could the authors comment on this ?\nOur apologies, in fact, Model 2 includes the subduction of the Australian continental lithosphere in\nthe Sunda margin. We have clarified this better in the text, where we now refer explicitly to the\nsubduction of the Australian continent along the Banda subduction zone segment.\nFor sake of simplicity, we have removed the complexities associated with the collision of the\nAustralian shelf in PNG and all the other subduction zones. This area is too complex to be properly\ncaptured. We have discussed the limitation of this approach in the dedicated section.\nAlso, what causes the diverging toroidal mantle flow in front of Sumatra slab ? Is it a\nconsequence of the trench advance ? I had in mind another toroidal flow from below the\nTibet to the Sunda area (e.g., Flower et al., 1998 ...) but here the mantle flows from\nThailand/Burma to the Tibet (fig.1d) which puzzles me ?\nWe have added more information in the text.\nFirst, to avoid confusion, we have reworded the flow as upwelling. The idea of “toroidal” flow is\nusually associated with the slab edge and retreat, as the reviewer correctly point out. In our context,\nit is the simple result of the sinking of the slab and a slab depth step, that is, an edge and the flow\naround it. This upwelling flow naturally arises at any edge of the slab, as the slab sinks, therefore,\ndoes not necessarily relate to the trench motions or slab edge. Hope this clarify.\nThen, the evidence for such upwelling flow in this area has been illustrated by shear wave splitting\n(Hu et al. 2023; Kong et al., 2020a) constraining the mantle upwelling from the sub-slab to the\nwedge, beneath southern Sumatra and adjacent areas. Here, receiver function investigations (Kong\net al., 2020b) also show that the 410 km and 660 km discontinuities are uplifted beneath Sumatra\nand depressed beneath Java, inferred to be related to the upwelling and downwelling flow,\nrespectively.\nWe have limited our references to only one of these papers.\nHu, H., Zhao, D., Lin, J. and Pilia, S., 2023. A Slab Window Beneath North Sumatra Revealed by\nP‐Wave Mantle Tomography. Journal of Geophysical Research: Solid Earth, 128(6),\np.e2022JB025976.\nKong, F., Gao, S.S., Liu, K.H., Zhang, J. and Li, J., 2020. Seismic anisotropy and mantle flow in\nthe Sumatra subduction zone constrained by shear wave splitting and receiver function\nanalyses. Geochemistry, Geophysics, Geosystems, 21(2), p.e2019GC008766.\nKong, F., Gao, S.S., Liu, K.H., Ding, W. and Li, J., 2020. Slab dehydration and mantle upwelling in\nthe vicinity of the Sumatra subduction zone: Evidence from receiver function imaging of mantle\ntransition zone discontinuities. Journal of Geophysical Research: Solid Earth, 125(9),\np.e2020JB019381\np7, line 29 : what does that mean : \"large compression in Sumatra extends as far as the\nIndonesian region\" ? Could the authors be more precise ?\nWe have improved the wording here and provide the actual values.\np.8 - lines 54-58 are a bit confusing because I am not sure to see backarc extension in all\nmodels except 3b in the Andaman back-arc region (pale blue in models 1,2, 3a, 3b)\nIt was difficult to see in the figure, because of the colour scale, the pale green and pale blue are too\nclose. We have amended the colour scale and expanded the text, to better specify the reduced area\nof the dilating region.\np.12 : did the authors sort the seismic events (interface only) in their plot (fig.6) ?\nYes. We have cleared this, now.\np.12 - lines 95-96 : I don't understand this sentence. Could the authors better explain ?\nRewritten.\nfew typo errors:\n- p.3 - line 91 : remove \"the\"\n- p.4 - line 8 : remove e in \"enucleated\"\n- p.7 - line 12 : toward the ..???\n- p.10 - line 25 : northern\np.13 - line 37 : these show ...???\np.14 - line 45 : megathrust\np.14 - line 63 : match\nWe have amended these all. Apologies.\nThe text would be easier to read if it were a little more structured, with subheadings for\nexample like \"Limitations\" p.14 - paragraph 46-56, \"conclusion\" p.14 - paragraph 57-71.\nThis, of course, was to respect length constraints, and we are happy to amend it.\nWe have now broken the text down into sections and briefly expanded the limitation section.\nThe fig.6 a,b,d,e are difficult to read, especially the rupture lengths and envelops whose\ncolor contrast makes it difficult to see them properly.\nWe have redrafted figure 6 to improve readability.\nTo conclude, this study deserves to be published after minor revisions.\nSerge Lallemand\nWe thank Serge Lallemand for the supportive and thorough review of our work.\nReviewer #3 (Remarks to the Author):\nThis cross-disciplinary study demonstrates that stresses along convergent margin primarily\nare controlled by the down-going slab’s buoyancy and the subduction-induced flow in the\nupper mantle and are as such influenced by dynamics operating throughout a 3-D space\nencompassing plates and mantle. The study nicely encompasses three model families to\ntest different aspects systematically and individually and relate them to various\nobservables of the Sunda margin region to provide these insights. The model has been\ntested in-depth, although detailed tests (such as a resolution test) are only accessible in a\nprevious (but open-access) publication. The underlying open-source code, Underworld, is\nstate-of-the-art and well tested. Finally, the study is nicely and clearly written, and the nice\nfigures are helpful, accurate and accessible.\nThis study demonstrates that models like the ones presented here can, and here do, close\nthe gap between subduction dynamics and the long-term (multi-earthquake cycle) strength\nof megathrusts. As such, I see it as a nice fit for publication in Nature Communications. I\nonly have some minor comments outlined below.\nWe are grateful to Fabio Crameri for the encouraging review.\nSpeaking of inter-disciplinarity: I think the presentation of this study, as intended for an\ninterdisciplinary journal like Nature Communications, should get a little more attention to\nthe text. As a reviewer, I only comment on the scientifically relevant bits, as for example on\nwhat is “deep subduction”: For me this is subduction at around the upper mantle transition\nzone (670 km depth). For others this is maybe everything below the surface. And it is not\nmade clear what this study’s “deep subduction” refers to here. Could this be clarified.\nIndeed, apologies for that, we have clarified this throughout the text.\nModel setup:\nThe models apparently employ a free-slip top boundary (according to line 512)—which\nobviously caught my eye—meaning that there is no surface topography, which causes the\nsubduction zone to behave significantly different to a model with a free surface (e.g.,\nthe natural state of such models is double sided subduction, rather than a realistic single-\nsided sinking – see Crameri et al., 2012) because the subducting plate is forced to stay flat\nand at the top boundary until the very location of the plate interface. Has this simplification\nbeen tested on how it affects the stress distribution at the plate interface? If not, could this\nbe discussed shortly for a fair representation of the model?\nCrameri, F., P. J. Tackley, I. Meilick, T. V. Gerya, and B. J. P. Kaus (2012), A free plate\nsurface and weak oceanic crust produce single-sided subduction on Earth, Geophys. Res.\nLett., 39(3), L03306, doi:10.1029/2011GL050046\nIndeed, this is something that we have addressed. We have added the test in the supplementary\ninformation. The tests of the top boundary conditions and the rheology are in order and were\nperformed in 2D to reduce the cost of the computation.\nJust to clarify, in the FEM method the computational stability criterion has no dependency on the\ndimensions. Therefore, the solution in 2D and 3D is the same, for the same problem.\nWe have added the test in the supplementary information, and we are confident that the reviewer\nwill be satisfied.\nMore in general, we consider this point to be technical and will not address it in the manuscript,\nhowever we’d like to further illustrate our understanding of this problem, for the sake of scientific\ndiscussion.\nThe double-sidedness does not occur in our models because of the material tracking in the yielding\ncrust.\nFirst, this is due to our computational method’s ability to track material through Lagrangian\nparticles. This allows to differentiate the crust of the downgoing plate, which localises the shear in a\n(relatively) thin layer where conditions for plasticity are met. While this better reproduces the\nnatural prototype, it also entails a single sidedness in our models.\nThen, the crust yielding allows the effective decoupling of the two plates, as well as the decoupling\nof the subducting plate from the top boundary. These minimise the drag between the plates and the\nexcess (vertical) stress on the boundary, which otherwise weakens the whole margin, and avoid the\ndouble-sided subduction.\nThese are the major difference with the work of Fabio Crameri, where the single material modelled\ndoes not differentiate subducting and upper plates and weakening occurs broadly over the\nconvergent margins resulting in double-sided downwelling.\nIndeed, the outcomes of Fabio Crameri’s work are of great value from a general geodynamics point\nof view, aiming at the self-consistent emergence of subduction from the convection of a single-\ncomponent mantle.\nInstead, our computational approach allows differentiating the materials, which is suited to the\ninvestigation of the problem raised here. This approach creates a self-consistent shearing margin, as\nopposed to the “weak zone” with constant lowered viscosity used in other approaches, e.g., Stadler\net al., 2010, or to a “weakening zone” as described above, which will weaken plastically the whole\nmargin, from downgoing to the overriding plate.\nMost importantly, our approach allows the accurate stress distribution across the margin (Capitanio\net al., 2010) as in the “weak zone” the low viscosity accommodates all the compression, as opposed\nto a more uniform stress magnitude across the margin.\nStadler, G., Gurnis, M., Burstedde, C., Wilcox, L.C., Alisic, L. and Ghattas, O., 2010. The\ndynamics of plate tectonics and mantle flow: From local to global scales. Science, 329(5995),\npp.1033-1038.\nCapitanio, F.A., Stegman, D.R., Moresi, L.N. and Sharples, W., 2010. Upper plate controls on deep\nsubduction, trench migrations and deformations at convergent margins. Tectonophysics, 483(1-2),\npp.80-92.\nLines 54-55: I don’t understand what the authors want to communicate with this sentence.\ne.g., what does “set the stage” mean here? Should it say: “integrating [stresses] across\nthese scales”? – could you clarify this?\nWe have refined the initial sentences, here. We agree the expression was not clear.\nLines 55-58: Same here, not sure that comes across as intended. “the loading of\nmegathrusts by subducting plate motions—which increases plate interface friction” – The\n“motion” (also, why plural?) would rather reduce plate interface friction, from my point of\nview (there is more friction if an interface is static). Or you write: “the loading [of stress]\nfaces scrutiny from recent high-energy earthquakes”. Instead it should say something like:\n“Our current understanding about xxx faces scrutiny from recent insights gained from xxx”.\nI think I get what you want to say here, but it is written too confusingly, so you might want\nto consider rewriting.\nSee above\nLine 78: What does “large-scale” mean? – “Crust-scale”, “Lithosphere-scale”? “Upper-\nmantle scale”?\nChanged into “whole mantle-scale”.\nLine 148: What does “deep” refer to here?\nChanged into “whole mantle”.\nLine 197: What stability? Numerical or otherwise?\nAmended. The goal is to achieve a stable subduction, avoiding any possible artefacts due to the\ninitial conditions, such as singularities.\nLine 198 and 376: “Kyrs” – Units don’t have a plural form. Use “kyr” for a scientific notion\nand to be in line with other instances in the manuscript.\nAmended.\nLine 445: Typo in “Megathrust”\nAmended\nLines 445-456: I appreciate this modelling/study limitation paragraph. I think that is very\nuseful for any reader to understand what has been done and achieved in this study. Of\ncourse this could be extended, but I think works here as is.\nWe agree, we have expanded this part and made it a separate section.\nFigures:\nThe Scientific colour maps (Crameri 2018) are nicely put in use here. – That is great to\nsee! Would you mind also citing them, as they are under a MIT license.\nCrameri, F. (2018). Scientific colour maps. Zenodo.\nhttps://doi.org/10.5281/zenodo.1243862\nAbsolutely, apologies for this!\nWe have added it in the acknowledgment. Let us know if an in-line referencing is to be preferred.\nIf the authors want to make their figures perfect (here on in future manuscripts), avoiding\nred&green combinations at similar lightness values (as here with lines) could additionally\nbe avoided. But this is a pro-level-suggestion, in particular in comparison to other current\npublications.\nNoted, thanks for the suggestion.\nThank you for a very interesting study I enjoyed reading!\nFabio Crameri\n20.11.2024\nThanks to Fabio Crameri for the thoughtful review."}, {"round": 2, "reviewer_comments": {"reviewer_1": "The point-by-point response to my previous comments seems to have been improved, although there are some parts that I\ncannot agree with or understand. After reading the revised ms, I have many new questions about the handling of the data to\nbe compared with the model, the selection of model parameters to be compared with the observational data, and other\naspects of the description. I will comment on the points I noticed below.\nMajor points\n1. The equations for the balance of forces such as ridge push and slab pull are not written. As I commented last time, write\nthem in the form of an equation. Also, I don't understand the relationship between that equation and the equation of\nmomentum. Write it so that the relationship between the two is clear.\n2. L.263-L.279 The results near the model boundary should be strongly influenced by the boundary conditions (free slip), so\nthis influence should also be mentioned. By comparing Model 1 and Model 2, it should be possible to evaluate the influence\nof the boundary conditions on the north and east sides of Model 1. By comparing Model 1 and Model 3, it should be possible\nto evaluate the influence of the boundary conditions on the south and west sides of Model 1. By comparing Model 1, 2 and\nModel 3, it should be possible to evaluate the influence of the boundary conditions on the lower side of Model 1 and 2.\n3. L.294-307 There is a model in which the areal strain is positive in the back arc. Whether the model is acceptable or not\nshould be possible to determine by examining the focal mechanism of earthquakes in the back-arc crust. Discussion should\nbe conducted on models that can explain the observed data.\n4. L.335-356 Is trench movement (advance and retreat) an input or an output of your computation? If it is an input, the\nposition of the trench at the start of the calculation should also be specified. In that case, the calculation time is assumed to\nbe infinite, so the movement distance becomes infinite, but I cannot understand the relationship with the description in this\nparagraph. If it is an output, is it a full dynamic model? In that case, was a weak zone set at the plate boundary and the\nsimulation performed? In this case, how did you carry out the subduction along the Slab2 model ?How can it be combined\nwith data such as NNR?\n5. L.357-380 Why is τ11 used as the evaluation criterion for stress? In the Sumatra subduction zone, even though oblique\nsubduction is occurring, the P-axes of the earthquakes are oriented in a direction perpendicular to the trench axis. Can this\nmodel explain such observational facts?\n6. L. 465-520 Why do you compare earthquakes with normal stress of the plate interface? It seems that earthquakes are\nextracted from earthquakes occurring at plate boundaries, but in that case, they should be compared with shear stress.\nSince shear stress is the site of effective friction coefficient and normal stress, shear stress also changes depending on the spatial\ndistribution of effective normal stress. Many studies (especially thermal structure modeling) have shown that the effective\nfriction coefficient is very small, and therefore, the shear stress is also small, so it is premature to think that the magnitude of\nnormal stress directly corresponds to seismic activity at plate boundaries. This should be thoroughly discussed.\nMinor points\nL.282 \"buoyancy\" should be \"negative buoyancy\".\nL. 290 O(-10^-13) The maximum absolute value of the scale bar is 5×10^-14. Make sure the text and figure are consistent.\nL.291 Draw a scale of 100 km in Figure 3a.\nL.294-307 All models show a region of positive areal strain on the ocean side of the Java subduction zone. Explain why.\nL.308-315 Does this mean that the observed data do not include elastic strain? If so, how will elastic strain accumulate for\nthe next megathrust earthquake?\nL. 358, 369 Upper mantle subduction Unclear. Rewrite.\nL. 409 norther What do you mean?\nL. 403 Isn't the model a viscoplastic model? Why does the elastic parameter \"rigidity\" appear here?\nL. 417 in-plane compression Explanation required\nL. 463, 499, 504, 526 fracture energy This is a term used in the slip weakening model for dynamic fault rupture. Why is the\nterm dynamic fault rupture used here, even though it is a long-term model?\nL. 476 Since this study estimates a long-term stress field, it should be compared with steady-state seismic activity. Are\naftershocks of large earthquakes properly removed using the ETAS model, etc.? Aftershocks of large earthquakes of Mw9.0\nor more can have aftershocks with magnitude of 6 or more. If earthquakes associated with such large earthquakes have not\nbeen removed, the discussion should be based on the removed seismic activity (steady-state seismic activity).\nL.523, 537, 546, 547 Are these values ​​when pore fluid pressure is not taken into account? In reality, the effective friction\ncoefficient should be smaller because pore fluid pressure is believed to exist at the plate boundary.\nL. 545 pdfs Explanation required\nL. 553 ~10s MPa What is this unit? Explanation required\nL.643 Have you not taken into consideration the effect of the phase transition boundaries at depths of 410 km and 660 km on\nthe behavior of the slab? If not, you should write the reason why the effect can be ignored.\nL. 648 η Is the value of the viscosity a constant or a variable? If it is a constant, write the value. If it is a variable, write it in the\nform of an equation, such as a function of what.\nL.651 I and τ The definitions are incorrect.\nL.658 I don't understand the relationship between ηcreep, η in this equation and η in the equation of momentum in L.643. It\nshould be written in a way that is easy to understand.\nSupplementary Materials\nFig. S3 is divided into four figures and four captions without symbols (a-c, d, and e), which is completely different from the\nusual notation and difficult to understand. It is desirable to change Fig. S3d and after to treat them as Fig. S4, etc.\nFig. S3 Caption: What is the difference between \"slab crust\" and \"oceanic crust\"? Is it the difference before and after the\noceanic plate subducts? Why do the two give rise to different rheology?\nFig. S11 It is strange that the shear stress is larger at the plate boundary in the deep ductile region than in the shallow\nregion. Explain the cause and interpretation of this.", "reviewer_2": "I am satisfied with the authors' responses, and I am therefore in favor of publishing this manuscript as is.", "reviewer_3": "The authors have made adjustments to the manuscript following the review comments and addressed all my previous\ncomments. I think showing the differences in overall dynamics between a free and a fixed model plate surface is helpful. I\nappreciate the 2-D test results being added to the supplementary material to see and be reminded of how, and how much,\nimposing a flat plate surface changes the overall geodynamics of a subduction zone. I also found the authors' detailed\nexplanations in the response useful to read. – Thank you.\nI have no further comments and recommend it for publication."}, "author_response": "We have done our best to improve the manuscript in this and the previous version.\nConsidering the little modification to the manuscript requested, mostly clarifications, we\nfocus this rebuttal letter on the scientific debate.\nThis is, by no means, meant to solve the disagreement, on the contrary. we welcome a certain\nlevel of disagreement: disagreement is vital to science.\nWe are thankful for the opportunity to discuss, scientific progress comes from questioning,\nfrom debating, from testing ideas, as we are doing here. Different interpretations,\nmethodological approaches, or gaps in our knowledge are beneficial to science as they drive\nfurther research and a deeper understanding.\nTo keep the debate constructive, we have made our best to provide experimental and\nempirical evidence, backed by data and reproducible results, and we discussed our logical\nreasoning, grounding our arguments in scientific principles.\nNow, achieving an agreement cannot be the goal here, indeed. Instead, it would be beneficial\nif the reviewer could agree to a respectful disagreement and let the debate reach the broader\nscientific community.\nSuch disagreements have always been part of the scientific community, just as others hold\ndifferent views and may agree, much like the reviewers of this manuscript.\nMajor points\n1. The equations for the balance of forces such as ridge push and slab pull are not\nwritten. As I commented last time, write them in the form of an equation. Also, I don't\nunderstand the relationship between that equation and the equation of momentum.\nWrite it so that the relationship between the two is clear.\nApologies. This set of equations is acknowledged in geodynamics, therefore frequently\nomitted. However, we understand that this might not be so familiar to researchers in other\nfields.\nWe have added a whole section in the Supplementary Information starting from the force\nbalance, describing the forces, and deriving the equation of momentum, then conclude with\nthe Stokes flow.\nHowever, because this is textbook knowledge, we prefer to leave it in the additional material\nand keep the method section of the manuscript succinct.\n2. L.263-L.279 The results near the model boundary should be strongly influenced\nby the boundary conditions (free slip), so this influence should also be mentioned. By\ncomparing Model 1 and Model 2, it should be possible to evaluate the influence of\nthe boundary conditions on the north and east sides of Model 1. By comparing Model\n1 and Model 3, it should be possible to evaluate the influence of the boundary\nconditions on the south and west sides of Model 1. By comparing Model 1, 2 and\nModel 3, it should be possible to evaluate the influence of the boundary conditions\non the lower side of Model 1 and 2.\nYes, the boundary forces and model boundary affect the results.\nThis is what we have done and is the goal of our paper, in fact: testing the effect of different\nboundary forces on the subduction margin stresses.\nThe boundary forces either arise from plate margins or the mantle, and are directly\nimplemented or simulated through internal or boundary conditions. This is customary in\ngeodynamics, indeed, and are known to affect the velocity, under certain conditions, but not\nthe stresses.\nThis is well known in fluid dynamics and geodynamics, and follows the Stoke sinker theory,\nwhere the drag force around the slab may increase due to wall proximities or viscosity,\nwhereas the stress at the subduction zone scale with the slab buoyancy. The driving force,\nand therefore the stress, are unaffected, but the resisting forces are, determining a different\nsettling velocity.\nAs an analogy with the Stokes sinker theory, in our subduction experiment the shear stress\nmay increase along the wall, as the reviewer suggests, but the effect is mostly notable on the\nforce balance around the sinking body’, affecting drag and settling velocity. We keep this\nanalogy in our work and measured similar quantities, velocities, strain and stress at the\nsubduction margin, as the boundary forces and boundary conditions change. This is the effect\nthe boundary conditions have. Measuring it along the models’ walls does not necessarily\ntranslate into an obvious quantity that can be compared to the observations.\nBecause this is an acknowledge fluid dynamics aspect, we prefer not to discuss it in our\nmanuscript, but we have further cleared throughout the manuscript. The tracked changed\nversion will clearly show our amendments, too diffuse to report here.\nIn support of this view, common to geodynamics, we report some of the papers that have\ninvestigate the application of the fundamental Stokes sinker theory to subduction, explained,\nfor instance, in section “3.1 Stokes flow and circulation model” of the review:\nBecker, T.W. and Faccenna, C., 2009. A review of the role of subduction dynamics for\nregional and global plate motions. Subduction zone geodynamics, pp.3-34.\nwhere the scaling of the stresses is presented. This view is grounded in previous work,\nincluding ours:\nJadamec, M.A. and Billen, M.I., 2010. Reconciling surface plate motions with rapid three-\ndimensional mantle flow around a slab edge. Nature, 465(7296), pp.338-341.\nJadamec, M.A. and Billen, M.I., 2012. The role of rheology and slab shape on rapid mantle\nflow: Three‐dimensional numerical models of the Alaska slab edge. Journal of Geophysical\nResearch: Solid Earth, 117(B2).\nCapitanio, F.A., Morra, G. and Goes, S., 2007. Dynamic models of downgoing plate-\nbuoyancy driven subduction: Subduction motions and energy dissipation. Earth and\nPlanetary Science Letters, 262(1-2), pp.284-297.\nGoes, S., Capitanio, F.A. and Morra, G., 2008. Evidence of lower-mantle slab penetration\nphases in plate motions. Nature, 451(7181), pp.981-984.\nCapitanio, F.A., Morra, G., Goes, S., Weinberg, R.F. and Moresi, L., 2010. India–Asia\nconvergence driven by the subduction of the Greater Indian continent. Nature Geoscience,\n3(2), pp.136-139.\nSchellart, W.P., Stegman, D.R. and Freeman, J., 2008. Global trench migration velocities\nand slab migration induced upper mantle volume fluxes: Constraints to find an Earth\nreference frame based on minimizing viscous dissipation. Earth-Science Reviews, 88(1-2),\npp.118-144.\nStadler, G., Gurnis, M., Burstedde, C., Wilcox, L.C., Alisic, L. and Ghattas, O., 2010. The\ndynamics of plate tectonics and mantle flow: From local to global scales. science, 329(5995),\npp.1033-1038.\n3. L.294-307 There is a model in which the areal strain is positive in the back arc.\nWhether the model is acceptable or not should be possible to determine by\nexamining the focal mechanism of earthquakes in the back-arc crust. Discussion\nshould be conducted on models that can explain the observed data.\nYes, we agree, we did present the comparison in the discussion section, in fact, and have\nexpanded.\nWe have added a figure where we show the stress regime in the back arc and compare it with\nthe tectonics and the interseismic strain rates. In this sense, our models do explain the\nobserved data.\nAs noted by the reviewer in the previous review, there is no ground for a comparison between\nour averaged strain and the individual seismic events. We do not look at the rupture\nmechanisms, generating the earthquakes, therefore the strain in the upper plate, revealed by\nseismicity, are not the goal of this work.\nFor the sake of the discussion, we show here a comparison between our model’s major stress\naxes, the P-axes recovered from individual earthquakes and the CMT in the back arc, as\nrequested. Note that interpreting the P-axes of the earthquakes as indicators of stress may\nnot be warranted, as 1) these are recovered from a strain rate tensor, not the stress, and 2)\nwe must assume that the rheology is known and linearly correlates the stress and the strain\nrates. These are common simplification in the interpretation of the seismicity that we will not\ndiscuss here.\nWe focus this comparison on the back-arc regime. Part of this discussion has been folded in\nthe manuscript.\nOur model main stress axes (i.e., the largest eigenvector of the stress tensor) illustrate the\ndominant transcurrent regime in the Andaman domain (left hand side panel). This is defined\nas by the maximum and minimum stresses, s (green axes) and s, in the horizontal plane,\n1 3\nwhile s is vertical. This regime extends further south, where it becomes compressive, in west\n1\naxes) and vertical s, with minor strike-slip in the eastern side of the islands.\n(green axes) is mostly confined in the northern sector of the margin, from Andaman to\nSumatra. While analysing the orientation is beyond the scope of this letter, these are aligned\nwith the major structural features, the Andaman spreading and the Sumatran fault, which\naccommodates ~30% of the interseismic deformation, and are active since the Neogene, and\nare discussed in our manuscript, supported by published work:\nMcCaffrey, R. The tectonic framework of the Sumatran subduction zone. Ann. Rev. of Earth\nand Plan, Sci. 37, 345–366 (2009).\nSieh, K. & Natawidjaja, D. Neotectonics of the Sumatran fault, Indonesia. J. Geophys. Res. -\nS.E. 105, 28295–28326 (2000).\nFurthermore, the compression in Sumatra has seismic evidence, as suggested by the\nreviewer, and shown in the paper we have referenced:\nWidiyantoro, S., Supendi, P. & Rawlinson, N. A note on the seismicity of Sumatra, western\nSunda Arc, Indonesia, in relation to the potential for back-arc thrusting. Sci. Rep. 14,\ndoi:https://doi.org/10.1038/s41598-024-64076-7 (2024).\nThe CMT mechanisms are another way to represent the P-axes and we will not discuss it\nfurther.\n4. L.335-356 Is trench movement (advance and retreat) an input or an output of your\ncomputation? If it is an input, the position of the trench at the start of the calculation\nshould also be specified. In that case, the calculation time is assumed to be infinite,\nso the movement distance becomes infinite, but I cannot understand the relationship\nwith the description in this paragraph.\nNo, the trench motions are not an input. Trench motions, stresses and strains are the result of\nthe force balance. We have added a sentence to clear this further.\nIf it is an output, is it a full dynamic model? In that case, was a weak zone set at the\nplate boundary and the simulation performed? In this case, how did you carry out the\nsubduction along the Slab2 model ?How can it be combined with data such as NNR?\nYes, it is the output of a fully dynamic model. We have added a sentence to clear this further.\nNo “weak zone” is implemented in our approach, instead we use a more natural set up which\ninclude a plastic crust which yields self-consistently under the effect of the subduction forces.\nThe “weak zone”, meaning a layer with an imposed fixed low viscosity or lowered yield\nstrength, leads to singularities in the stress field. Instead, our setup leads to more consistent\nstress distribution across the margin, while best reproducing the natural configuration. The\n“weak zone” is usually used for initiating subduction, then the perturbation can relax. In our\nmodel it is crucial to minimise the introduction of artefact, while maximising the\nsimplification for a robust understanding.\nWe have tested this several times, among which:\nCapitanio, F.A., Stegman, D.R., Moresi, L.N. and Sharples, W., 2010. Upper plate controls\non deep subduction, trench migrations and deformations at convergent\nmargins. Tectonophysics, 483(1-2), pp.80-92.\nCapitanio, F.A., Faccenna, C., Zlotnik, S. and Stegman, D.R., 2011. Subduction dynamics\nand the origin of Andean orogeny and the Bolivian orocline. Nature, 480(7375), pp.83-86.\nThe Slab2 model is used as an initial condition, that is, the geometry. The lithospheric\nthickness is then constrained by the age of the oceanic crust and by the Slab2 model, which is\nplastic in our models, as in nature. This is explained in the manuscript and further illustrated\nthe supplementary information. Finally, we solve for the Stokes flow and obtain a fully\ndynamic subduction system with realistic geometries, motions, stress and strain rates.\nOur approach is similar to many others, where realistic geometry through these datasets are\nused, e.g.:\nFaccenna, C. and Becker, T.W., 2010. Shaping mobile belts by small-scale\nconvection. Nature, 465(7298), pp.602-605.\n5. L.357-380 Why is τ11 used as the evaluation criterion for stress? In the Sumatra\nsubduction zone, even though oblique subduction is occurring, the P-axes of the\nearthquakes are oriented in a direction perpendicular to the trench axis. Can this\nmodel explain such observational facts?\nWe have added a figure to the manuscript, as this is a point that may be misinterpreted.\nWe use here the (square root of the) second invariant of the stress tensor. This is a scalar, not\na vector, and allows estimating stress magnitude.\nWe have now added the eigenvectors of the stress tensor, and plot the longest horizontal\nvector, aka the maximum stress.\na Model3c b – Δα= 8.9˚\n10°N\n0 0\n0 0\n2 2\n0°\n9 9\n0 0\n10°S\n600 6 0\n0\n200 2 0\n0\n400 400\nThrust Azimuth misfit Δα (˚)\n40 30 20 10 0\nStrike-slip\n110°E 120°E 110°E 120°E\nThe figure above (now figure 6 in the manuscript) illustrates the models main surface stress\ns (left panel) side by side with a comparisons of megathrust stress and P-axes from the\n1\nseismicity.\nAs the reviewer correctly points out, the stress axes are perpendicular to the trench, in both\nthe observations and the models. The colour coding in (b) shows the misfit azimuth between\nmodels’ and P-axes, with a mean misfit of only 8.9˚.\nThese are comparable as the former is the modelled is the stress regime, the observed is the\nsame stress regime, yet inferred through the earthquake CMT. As indicated above, we cannot\ncompare further these as we do not model the constitutive laws correlating the stress regime\nwith the strain regime.\nYet, they are (might) close to equivalent in our models as the result of the simplified rheology\nadopted.\nWe also note that our models do reproduce the convergence obliquity along the trench (Fig.\nS8).\nHowever, the reviewer will agree that the kinematic frame for seismicity is not really viable,\ntherefore the analysis has valid structural implications, yet does not illustrate the dynamics,\nthe force balance at play, and might be misleading. The relation between plate motions and\nseismicity has been, in fact, repeatedly confuted, as we mention in the manuscript.\nThis area is, in fact, the exemplar case where this is true, where the trench-normal\nmegathrust P-axes do not easily reconcile with the obliquity convergence and oblique seismic\nslip in the fore arc and back arc. As in figure 4 from:\nMcCaffrey, R., 2009. The tectonic framework of the Sumatran subduction zone. Annual\nReview of Earth and Planetary Sciences, 37(1), pp.345-366.\n6. L. 465-520 Why do you compare earthquakes with normal stress of the plate\ninterface? It seems that earthquakes are extracted from earthquakes occurring at\nplate boundaries, but in that case, they should be compared with shear stress.\nSince shear stress is the site of effective friction coefficient and normal stress, shear stress also changes depending on the spatial\ndistribution of effective normal stress. Many studies (especially thermal structure modeling) have shown that the effective\nfriction coefficient is very small, and therefore, the shear stress is also small, so it is\npremature to think that the magnitude of normal stress directly corresponds to\nseismic activity at plate boundaries. This should be thoroughly discussed.\nWe have discussed this point in our previous letter and in the manuscript.\nThe modelling of friction in seismicity usually follows the simple assumption of lithostatic\npressure, then, given the dip of a thrust, the normal and shear components are determined.\nUnder this condition, the friction is nothing else than a trigonometric function of the dip, and\nis customary to neglect the normal stress to illustrate the shear stress.\nOur work goes beyond this limitation and illustrates the “tectonic” friction, that is the ratio\nof the shear-to-normal stresses, yet, these two are independent. The shear stress depends on\nthe constitutive laws of yielding crust, the normal stresses depends on the whole subduction\ndynamics.\nWe have cleared this further in line 528 and ff:\n“…this is the ratio of deviatoric stress due to subduction and plastic stress, therefore differs\nfrom what commonly calculated as the result of lithostatic pressure”\nThen, the normal stress and friction allows to illustrate two different points.\nThe normal stresses are compared to the inferred barriers and to patches that are locked or\nwould – ideally – yield at higher stresses, as we mention in the manuscript. The comparison\nwith this value provides insights into areas that have no major/long rupture.\nThe friction embeds the shear stress in the ratio and supports the discussion of the conditions\nin areas that did rupture. The shear stresses are shown in a supplementary figure, but for\nlimitations to the length of the manuscript we preferred to leave it there.\nWe agree that none of these stresses corresponds to seismic activity at plate boundaries.\nWe have carefully checked our work, there is no such claim.\nInstead, our major conclusion is that the tectonic (deviatoric) stresses, usually neglected, are\nin the same order of those inferred through seismological tools.\nThese stresses are to be considered complementary to the calculation suggested by the\nreviewer and may provide a future direction to refine and reconcile the different estimates.\nMinor points\nL.282 \"buoyancy\" should be \"negative buoyancy\".\nAmended accordingly\nL. 290 O(-10^-13) The maximum absolute value of the scale bar is 5×10^-14. Make\nsure the text and figure are consistent.\nApologies, that was a typo. It is amended, now\nL.291 Draw a scale of 100 km in Figure 3a.\nDone\nL.294-307 All models show a region of positive areal strain on the ocean side of the\nJava subduction zone. Explain why.\nThis is the forebulge in subduction zones, which is known to be under extension. It is larger\nin the Java segment as this is where the slab is longest, therefore the bending moment\nhighest. This decreases further north as the slab is shorter.\nWe have added a sentence and two references to work that have addressed this very topic.\nThis is a known feature of subduction zones and we believe this does not need further\nexplanation in the text.\nL.308-315 Does this mean that the observed data do not include elastic strain?\nIf so,\nhow will elastic strain accumulate for the next megathrust earthquake?\nThese questions may not be pertinent and go beyond what addressed here.\nThe interseismic strain rate is time dependent. Instead, the elastic strain is time independent,\nthey cannot show in the same observable.\nThen, it is known that elastic strain accumulates on megathrusts under the action of plate\nmotions.\nBecause on Earth the elastic moduli are isotropic, the elastic strain coincides with the elastic\nstress, therefore, the elastic stress is loaded by the motions of tectonic plates. This is also\nreferred to as stress loading or tectonic loading, and is what is modelled here: the modelled\nplate motions load the plate interface, the plate interface stress applies from the lower to the\nupper plate, deforming it.\nAt shorter timescales, transient elastic stresses are accumulated then released. This\ncomponent is removed from our modelling and from the observables we have used.\nThis is fundamental tenet of seismotectonics.\nWe are not sure this is what the reviewer questions, here.\nApologies if we have misunderstood.\nL. 358, 369 Upper mantle subduction Unclear. Rewrite.\nAmended to “slab subduction in the upper mantle”\nL. 409 norther What do you mean?\nApologies, that was a typo, now changed into “northern”\nL. 403 Isn't the model a viscoplastic model? Why does the elastic parameter \"rigidity\"\nappear here?\nAmended, changed into strength. It was meant to express the ability of plates to sustain large\nstress. Strength is more appropriate, agreed, and avoid confusion.\nL. 417 in-plane compression Explanation required\nChanged to “compression in the plane of the lithosphere”\nL. 463, 499, 504, 526 fracture energy This is a term used in the slip weakening\nmodel for dynamic fault rupture. Why is the term dynamic fault rupture used here,\neven though it is a long-term model?\nThe term is used in this speculative part of our work, where we try to extend the insights our\noutcomes provide.\nTo avoid confusion, we have reduced the discussion of this point and removed the relative\nsentences.\nL. 476 Since this study estimates a long-term stress field, it should be compared with\nsteady-state seismic activity. Are aftershocks of large earthquakes properly removed\nusing the ETAS model, etc.? Aftershocks of large earthquakes of Mw9.0 or more can\nhave aftershocks with magnitude of 6 or more. If earthquakes associated with such\nlarge earthquakes have not been removed, the discussion should be based on the\nremoved seismic activity (steady-state seismic activity).\nThis comment was already raised in the previous letter and addressed to the best of our\nabilities. While we provide here further explanation, we note that this will be the same\nprovided before. We are hoping for a more proficient discussion, this time.\nFirst, we are not sure what steady-state seismicity means. The reviewer likely refers to\nsteady-state deformation not associated with earthquakes. These are the interseimic strain\nrates, or also the long-term strain rate. As we have discussed before, this is what we have\nused, in fact.\nMost importantly, our outcomes do not rely on this dataset only.\nWe have explained in the text that our models compare to other observables, mostly the long-\nterm tectonics in the area. Additionally, we have dedicated a publication to this comparison,\nZuhair et al, 2022.\nThe dataset we use does remove stations and periods close to major earthquakes, as we have\nexplained in our previous letter.\nThe dataset used is the GEM Strain Rate Model, supported by the Global Earthquake Model\nfoundation and UNAVCO (a joint collaboration of NSF, USGS and NASA).\nThe fundamental methodology and tests have been published, and are chiefly explained here:\nKreemer, C., Holt, W.E. and Haines, A.J., 2003. An integrated global model of present-day\nplate motions and plate boundary deformation. Geophysical Journal International, 154(1),\npp.8-34.\nKreemer, C., Blewitt, G. and Klein, E.C., 2014. A geodetic plate motion and Global Strain\nRate Model. Geochemistry, Geophysics, Geosystems, 15(10), pp.3849-3889.\nL.523, 537, 546, 547 Are these values when pore fluid pressure is not taken into\naccount? In reality, the effective friction coefficient should be smaller because pore\nfluid pressure is believed to exist at the plate boundary.\nCorrect. The pore and lithostatic pressures are removed. This is explained throughout the\ntext, and we have explicitly cleared it, in this paragraph.\nWe agree with the reviewer, the pore pressure and the lithostatic pressure are indeed present,\nas well as the tectonic pressure.\nHere, we have expressly focused on the latter, we have removed the lithostatic and pore\npressure to illustrate the component of friction due to tectonic forces.\nThe conclusion of our work is that the tectonics stresses are in the same order of those\nconstrained by other methodologies. This is compatible with the idea of an arbitrary value of\nnet pressure, which may encompass the tectonic deviatoric stresses. These are commonly\nneglected, under the assumption of a quasi-lithostatic pore pressure. However, this\nassumption cannot be constrained easily.\nBy constraining one component of the pressure, that due to the force balance of subduction,\nour work helps constraining the other two, pore and lithostatic pressure.\nL. 545 pdfs Explanation required\nWe have added an explanation of the pdfs.\nL. 553 ~10s MPa What is this unit? Explanation required\nAmended, was meant to be “tens” of MPa, now reads O(10) MPa.\nL.643 Have you not taken into consideration the effect of the phase transition\nboundaries at depths of 410 km and 660 km on the behavior of the slab? If not, you\nshould write the reason why the effect can be ignored.\nThe phase changes at 410 and 660 km are included in the density scaling adopted here.\nThe averaged densities, or better the density contrasts, are constrained, together with the\nviscosities, by means of gravity, topography, plate motions, and other observables. Therefore,\nthe density contrasts constrained embeds inherently the phase changes. These are\nfundamental to the problem solved here, the Stokes flow, as constraining the mass anomaly\nallows to define uniquely the viscosity.\nThis is a customary approach in geodynamics:\nWada, I. and King, S., 2015. Dynamics of subducting slabs: Numerical modeling and\nconstraints from seismology, geoid, topography, geochemistry, and petrology. In Treatise on\nGeophysics, Vol. 7, Mantle Dynamics (pp. 339-391). Elsevier Inc.\nL. 648 η Is the value of the viscosity a constant or a variable? If it is a constant, write\nthe value. If it is a variable, write it in the form of an equation, such as a function of\nwhat.\nIt is variable and the equation is explained in detail in the following paragraph, lines 660-\n684.\nNote that this is the customary way the Stokes flow is illustrated: first we describe the\ngoverning equations, then we describe the individual functions. These are one after each\nother in our method section, therefore accessible.\nL.651 I and τ The definitions are incorrect.\nThey were inverted, apologies. We have now amended the text.\nWe also have improved consistency between manuscript and supplementary information, and\nnow use s, s’ and t the stress, the deviatoric stress and the shear stress, respectively.\nL.658 I don't understand the relationship between ηcreep, η in this equation and η in\nthe equation of momentum in L.643. It should be written in a way that is easy to\nunderstand.\nAs in the comment above, the h in the equation of momentum is explained throughout this\nparagraph. We have added a description at the beginning of the paragraph.\nWe are confident this will better lead the reader through the details of the implementation.\nSupplementary Materials\nFig. S3 is divided into four figures and four captions without symbols (a-c, d, and e),\nwhich is completely different from the usual notation and difficult to understand. It is\ndesirable to change Fig. S3d and after to treat them as Fig. S4, etc.\nWe have re-labelled all the panels in all figures, to improve consistency.\nRegarding figure S3, we prefer to keep all this material into a single figure, though, for\nreadability, they all relate to the same numerical tests and are not central to the content of\nthe paper.\nFig. S3 Caption: What is the difference between \"slab crust\" and \"oceanic crust\"?\nIs\nit the difference before and after the oceanic plate subducts? Why do the two give\nrise to different rheology?\nThis is just the labelling of the materials, all the crusts share the same rheology. This might\nnot have been clear, apologies. We have amended the caption, labelled the plot and add the\nreference to the crustal rheology, in the main text.\nThe difference between crusts is due to the geometry of the model set up: the oceanic crust\nthickness is defined by the age of the oceanic lithosphere, for age < 80 My, then is constant;\nthe thickness of the slab’s crust is equal to that of the oceanic plate at the trench, and so is\nthe upper plate’s crust is to reduce the complexity of the model.\nFig. S11 It is strange that the shear stress is larger at the plate boundary in the deep\nductile region than in the shallow region. Explain the cause and interpretation of this.\nWe are not sure we understand this comment, likely this is a misunderstanding.\nThe shear stress is is independent of the depth, being a combination of the off-diagonal terms\nof the stress tensor. The variation of shear stress, then, depends on the interface\nconfiguration, with no expected depth-dependence.\nThen, we think the reviewer is referring to what is constrained by heat flow-shear heating\nmodelling and illustrate a figure from this work:\nGao, X. and Wang, K., 2014. Strength of stick-slip and creeping subduction megathrusts from\nheat flow observations. Science, 345(6200), pp.1038-1041.\nThis figure shows that the inferred shear stress increases at depth in the Japan trench, and\nincreases down to 20 km, then decreases, in the northern Hikurangi. To our understanding,\nboth behaviours are plausible. What is the shear stress deeper than what the published\nanalyses can constrain is not known.\nWe show that the increase of shear stress in our models is in agreement with what published\nbefore in the geodynamics community. Here, we show how this varies along the trench in a\nthree-dimensional setting, which is novel.\nWe show below published two-dimensional models, from our and other groups. These\nmodels show how the deviatoric stress, that is, net of pressure, increases with depth along the\ninterface. In Capitanio et al., 2011, this is illustrated through the (square root of the) second\ninvariant of the stress tensor s , in Schellart and Moresi, 2013, and Holt et al., 2015, by the\nII\ndeviatoric horizontal stress, s , and in Beall et al., 2021, as the deviatoric shear stress, t , as\nxx\nin our manuscript.\nCapitanio, F.A., Faccenna, C., Zlotnik, S. and Stegman, D.R., 2011. Subduction dynamics and\nthe origin of Andean orogeny and the Bolivian orocline. Nature, 480(7375), pp.83-86.\nSchellart, W.P., and L. Moresi. \"A new driving mechanism for backarc extension and backarc\nshortening through slab sinking induced toroidal and poloidal mantle flow: Results from\ndynamic subduction models with an overriding plate.\" Journal of Geophysical Research:\nSolid Earth 118, no. 6 (2013): 3221-3248.\nHolt, A.F., Becker, T.W. and Buffett, B.A., 2015. Trench migration and overriding plate stress\nin dynamic subduction models. Geophysical Journal International, 201(1), pp.172-192.\nBeall, A., Fagereng, Å., Davies, J.H., Garel, F. and Davies, D.R., 2021. Influence of\nsubduction zone dynamics on interface shear stress and potential relationship with\nseismogenic behavior. Geochemistry, Geophysics, Geosystems, 22(2), p.e2020GC009267.\nReviewer #2 (Remarks to the Author):\nI am satisfied with the authors' responses, and I am therefore in favor of publishing\nthis manuscript as is.\nReviewer #3 (Remarks to the Author):\nThe authors have made adjustments to the manuscript following the review\ncomments and addressed all my previous comments. I think showing the differences\nin overall dynamics between a free and a fixed model plate surface is helpful. I\nappreciate the 2-D test results being added to the supplementary material to see and\nbe reminded of how, and how much, imposing a flat plate surface changes the\noverall geodynamics of a subduction zone. I also found the authors' detailed\nexplanations in the response useful to read. – Thank you.\nI have no further comments and recommend it for publication.\nWe are thankful to Reviewer #2 and #3\nREVIEWER COMMENTS\nReviewer #3 (Remarks to the Author):\nThis is a 3rd review regarding comments from two other reviewers and the authors'\nresponse to one of them.\nWe are grateful for the time dedicated to our manuscript\nRegarding Reviewer #1 comments:\n1) OK; has been addressed.\n2) This is a fair point raised: there is no clear testing and discussing of most model\nboundaries (only the top one has now been tested and mentioned) even though they\ncan influence the flow in the model domain. To test these, usually one can simply\nextend the model domain and rerun the same model. Here, this is done too, but by\nextending the model domain, the authors additionally add strongly varying\ngeodynamic features, instead of just extending the given physical features on the\nsides. Which is what they do to study the natural effects beyond the model 1 extend,\nrather than the artificial ones from the model limits. The suggestion of the reviewer\nto test the latter with the model setups 1, 2, and 3, would not yield very useful\ntest results, as the authors rightly point out in the rebuttal. The authors argue that they\ntherefore test it closer to the central part of the domain, where their scientific focus\nlies. That does test the impact of the added geodynamic features together with any\nartificial boundary effect, but it does not extract the latter for actual testing. The\nauthors’ response here suggest a reliance on the general theory. But citing general\nprinciples does not entirely replace demonstrating how boundary effects play out in\nspecific model results either. So, I think it would indeed be fair to, at least, mention\nthe point how the side boundaries (i.e., the spatially limited model domain) might – or\nmight not – affect the model results also from a model limitation standpoint, as was\ndone for the top boundary condition (see lines 70–73). This would clarify to the\nreader that there might be these effects, even though they are likely localised (e.g.,\nnear walls only), if I had to guess.\nAs the reviewer suggests, we have cleared this aspect in the text.\nIn fact, this test is embedded in our analysis, as expressed before. Clearing further such a\ntechnical aspect for less expert readers could go beyond the scope of the paper.\n3) I would say this is taken care of to the extent that’s sensible here.\n4) OK: this has been clarified.\n5) OK: this has been clarified.\n6) OK: this has been clarified.\nWhile I may not have the expertise to judge all of reviewer #1's minor points\nperfectly, they appear to have been adequately addressed by the authors.\n––––\nRegarding Reviewer #4 comments:\nThe authors write that they tested Non-Newtonian (/composite) rheologies and the\nresults are provided in supplementary figures S3 (of which there are a few with the\nsame tag that belong to the same test, I believe). It is not clear to me whether these\nwere simply 2-D models or not, but the authors claim that they yield “comparable”\nresults and judging based on the figures, I think I would agree. So, I would say this is\nfine as is, also because it is clarified to the reader that the model eventually used is\nNewtonian.\nWe have clarified that the tests where run in 2D.\n––––\nOverall, if I may, I think this manuscript should have presented the model setup and\nits limitation much clearer and fairer from the beginning, which would have omitted\nconfusion and saved everyone (editors, reviewers, and ultimately also the authors) a\nlot of time. From my point of view, and including my comments above, the comments\nfrom reviewers #1 and #4 have been addressed though.\nWe agree that more information is needed.\nWe have followed the logic of providing information as is customary in computational\nmodelling, then adding the information as the reviewers deemed necessary.\nAs demonstrated here, the amount and type of information needed may depend on the\nreviewer/reader expertise and interest in the field, so that anticipating what may be perceived\nas “clear” is not always possible.\nWe apologise if this came across as a lack of information and thank the reviewer for the time\nspent on our work.\nReviewer #4 (Remarks to the Author):\nI read with interest the manuscript \"Bridging the gap between subduction dynamics\nand the long-term strength of megathrusts.\" In this work the authors perform a 3D\nnumerical model of subduction along the Java-Sumatra-Andaman megathrust.\nModeling is centered on the long-term stresses on the megathrust, and is an extension\nof previous studies to include a self-consistent plate interface constitutive behavior.\nThe authors test several numerical setups and confirm that Java slab pull and mantle\nflow dominate the margin’s tectonics. I noticed that the authors highlight an interesting\nobservation: the northern Sumatra segment shows an inverse correlation: despite\nhaving the thinnest, shortest slab and weakest slab pull, it exhibits the most intense\nseismic and tectonic activity. On the contrary, the southern Java segment\ncharacterized by a thicker and older slab has a relatively reduced seismic activity.\nAlthough this looks like it is the core of this work, I do not see it reflected in the abstract.\nThe manuscript is well, written, although I find it a bit longer and sometimes difficult to\nfollow. Figures are of good quality and suitable for a top journal. Figure 2 combines a\n3D spherical view with complex mantle flow patterns, but I found it difficult to\nunderstand and see details along the thrust margin. Figure 3 shows dilatation for\nseveral models, but just looking at the figure it is difficult to see the main differences\namong the models. Regarding the main text, there are so many observations the\nmodel fit that one can easily disperse the attention of the main goal of the manuscript.\nAlso, I found that just assuming a Newtonian mantle strongly influence the solution\nalong the plastic interface. I totally understand that this is a necessary modeling\nlimitation (given the size and resolution of the model, and also the available\ncomputational resources), nevertheless linking the long-term slab dynamics with\nshort-scale strength along the megathrust would have required a more precise\napproach.\nThis statement is at odd with what demonstrated in Figure S3.\nThere we show that the mantle flow and tractions are affected by the choice of mantle\nrheology, indeed, yet the plastic interface is not.\nThe reason why this occurs is simply explained by the threshold behaviour of the crustal\nrheology, which, by definition, cannot sustain more stress than its yield strength. This\nbecomes insensitive to stress variation, once has yielded. The main driver of the yielding is\nthe pull of the heavy slab attached. This mantle rheology results in yielding of the crust, thus\nthe interface has little or no sensitivity to the additional (minor) stresses imposed by\ngradients in the flow, modulated by the non-linearities in the rheology.\nWe have cleared this further in the text.\nMy main recommendations for this work are the following:\n•⁠ ⁠Title: this is a study focused on a specific subduction region, so this should be\nmentioned in the title.\nDone\n•⁠ ⁠Text (main): it is too long, has a great diversity of information that precludes the\nfocus needed for a high-impact finding(s). It is easy to split this work in several\nworks, make them shorter and more focused, each one concentrating on a specific\ntarget, not all at once.\nWe have improved the manuscript. We have now focus and strengthened on the most relevant\nfindings by rearranging slightly the manuscript.\nWe now start showing how the models reproduce the available observations, this is what\nremains unexplained in the current literature.\nThen show how the validated models can provide relevant insights into the current\nseismotectonics.\n•⁠ ⁠Figures: several figures must be updated with 2D cuts of 3D zoom in regions, so\none can better comprehend the information (see figure 2). Other figures need to\nbetter highlight he model differences.\nWe have moved a 2D section as requested in Figure 2, while we refer the readers to the\nSupplementary Information for more cross section of the mantle flow. Adding these in the\nmanuscript would lengthen the paper unnecessarily.\nWe have improved the figures overall and we are confident that the model differences are\nnow clearer."}]

Decision letter