RESEARCH ARTICLE Slow Geodynamics and Fast Morphotectonics in the Far East

10.1029/2021GC010167
Tethys
Special Section: L. Husson1 , N. Riel2, S. Aribowo1,3 , C. Authemayou4 , G. de Gelder1, B. J. P. Kaus2 ,
Tethyan dynamics: from rifting C. Mallard5, D. H. Natawidjaja3, K. Pedoja6, and A. C. Sarr7
to collision
1
Institut des Sciences de la Terre, CNRS, Université Grenoble Alpes, Grenoble, France, 2Institute of Geosciences, Johannes
Key Points: Gutenberg University Mainz, Mainz, Germany, 3Research Center for Geotechnology, Indonesian Institute of Sciences (LIPI),
• F ast and short-lived disruption in Bandung, Indonesia, 4LGO, IUEM, CNRS, Université de Brest, Plouzané, France, 5School of Geosciences, University of
surface morphotectonics occur at all Sydney, Sydney, NSW, Australia, 6M2C, Université de Caen, Caen, France, 7CEREGE, Aix-Marseille Université, Aix-en-
spatial scales during transient phases
Provence, France
of mantle flow
• Uncommon SE Asia seismotectonic
activity, dynamic topography, and
vertical land motion all triggered by Abstract  How can the sluggish, long-wavelength mantle convection be expressed by so many time
the transience of mantle flow and space scales of morphotectonic activity? To investigate these relationships, we explore the Java-Banda
• Depiction of the intricate relationships
subduction zone, where geodynamic records cluster. In the far-East Tethys, the exceptionally arcuate Banda
between mantle dynamics and surface
morphotectonics subduction zone circumscribes the deepest oceanic basin on Earth, seismotectonic activity slices the upper plate
more efficiently than anywhere else, and uncommonly vast expanses of continents are flooded in Sundaland
and Northern Australia. By comparing numerical simulations of subduction dynamics to a set of independent
Supporting Information:
Supporting Information may be found in
observations, we reveal the many facets of tectonic and physiographic changes that the sole docking of the
the online version of this article. Australian continent onto the subduction zone triggered. While mantle flow remains slow and long-wavelength
at depth, intense tectonic activity, and dynamic uplift and subsidence profoundly rework the physiography at
Correspondence to: many spatial scales. These results demonstrate that a modest disruption in the slow geodynamic tempo may
L. Husson, trigger manifold morphotectonic disturbances.
[email protected]

1. Introduction
Citation:
Husson, L., Riel, N., Aribowo, S., “Plate tectonics is the surface expression of mantle convection” is commonplace, but ever since Holmes (1931)
Authemayou, C., de Gelder, G., Kaus, seeded this concept, geological observations seem to defy this convenient rationalization of plate tectonics and
B. J. P., et al. (2022). Slow geodynamics
and fast morphotectonics in the far East mantle flow. In the vicinity of plate boundaries, the morphotectonic evolution of active plate boundaries occurs
Tethys. Geochemistry, Geophysics, at length and timescales that are seemingly at odds with the long-wavelength and slow pace of mantle convection.
Geosystems, 23, e2021GC010167. https:// Therefore, the links between mantle convection and surface tectonics remain challenging to disentangle (Davies
doi.org/10.1029/2021GC010167
et al., 2019; Jolivet et al., 2018; Mallard et al., 2016). At the far East end of the Tethys subduction zone, the
Received 14 SEP 2021 Java-Banda subduction zone displays a unique collection of morphotectonic events that occur at various time
Accepted 9 JAN 2022 and space scales that epitomizes this paradox. Whether these events are caused by the congestion of profuse
subductions and collisions is plausible, but only few studies attempt to comprehensively analyze this possibility
(Hall, 2017).

Perhaps the most puzzling observation is the pervasive flooding of vast continental units: Sahul in North Aus-
tralia, and Sundaland in West Indonesia, are the largest expanses of continents to lie below sea level on Earth
and are currently subsiding (Sarr, Husson, et al., 2019; Solihuddin et al., 2015). Another salient characteristic
of the region is the Weber Deep, which lies at depths of more than 7 km, making it the deepest region on Earth,
besides oceanic trenches. After the Banda slab rolled back to the eastern dead-end of the Banda embayment
(Milsom, 2001; Spakman & Hall, 2010), the Banda basin is now nested and circumscribed by the extremely tight,
horseshoed Banda arc (Figure 1). Intense seismotectonic activity is yet another remarkable facet of the region, the
back-arc domain being more commonly hit by devastating earthquakes (like the Flores, Alor, and Lombok earth-
quakes, Figure 1) than in any other subduction zone. Each of the above geodynamic oddities has been scrutinized
© 2022 The Authors. in one way or another but seldom as components, or surface expressions, of a unique, coupled mantle-lithosphere
This is an open access article under system. Here instead, we hypothesize that a single event, the inception of the Australian continent into the Ja-
the terms of the Creative Commons va-Banda trench, is in fact sufficient to explain the development of multiple morphotectonic events that altogether
Attribution-NonCommercial License,
which permits use, distribution and profoundly and swiftly remodeled the physiographic layout of the entire region.
reproduction in any medium, provided the
original work is properly cited and is not The Indian-Australian oceanic lithosphere subducted underneath SE Asia during the Cenozoic closure of the
used for commercial purposes. Neo-Tethys. This regime prevailed until the inception of the continent in the subduction zone in the East during

HUSSON ET AL. 1 of 17
 15252027, 2022, 1, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2021GC010167 by Nat Prov Indonesia, Wiley Online Library on [11/04/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
Geochemistry, Geophysics, Geosystems 10.1029/2021GC010167

Figure 1.  Morphotectonic map of SE Asia, southeastward view. The map depicts tectonic structures (faults: red curves; subduction trenches: red bold curves), envelope
of the subducting plate (depth isocontours 20 m), uplift and subsidence rates (triangles, see Section 2), selected focal mechanisms (Mw > 6.5, Wetar-Flores back-arc
thrust) and geodetic velocities (Koulali et al., 2016). Geodynamic extremes include [1] Papua orogen, [2] Weber superdeep basin, [3] subsiding Sahul platform, [4]
Timor slab hole, [5] Savu Sea, [6] Flores back-arc thrust, [7] Lombok 2018, Alor 2004, Flores 1992 earthquakes, [8] Java back-arc thrusts (Baribis and Kendeng belts),
and [9] subsiding Sunda shelf.

the Miocene (Hall,  2017), while oceanic subduction continued to the West. At the onset of collision of the
Australian continent with the island arc, the lithospheres on both sides of the subduction zone were paved with
alternating oceanic and continental lithospheres. The subducting plate was composed of the Indian oceanic lith-
osphere in the West and the Australian continent in the East, respectively converging with continental Sundaland
and a mostly oceanic domain. At present, the tessellated tectonic layout mostly remains, except that the northeast-
ern quadrant has been intensively reworked to form the complex assemblage of continental and oceanic units of
Wallacea (Figure 1). Here, by means of three-dimensional numerical models of mantle and lithosphere dynamics
(Kaus et al., 2016) that we match to a set of surface and subsurface tectonic, physiographic, and kinematic obser-
vations, we analyze the multiple consequences of the transient regime that followed the arrival of the Australian
continent into the subduction zone. We focus on the remodeling of the large-scale mantle and lithosphere dynam-
ics and how this disruption triggered a collection of surface morphotectonic events that are commonly analyzed
in isolation from one another. Besides reproducing the first-order regional features of dynamic topography and
regional plate tectonics, the simulations also outline the many time and space scales of the surface responses to
the sluggish dynamics of the Earth interior.

2.  Materials and Methods

2.1.  Numerical Model

We perform 3-D geodynamic simulations using LaMEM (Kaus et  al.,  2016). This finite difference staggered
grid discretization code uses particle-in-cell methods to solve the energy, momentum, and mass conservation
equations. The rheologies of the rocks are assumed to be visco-elasto-plastic and the total deviatoric strain rate
is given by:
1 1 𝜕𝜕𝜕𝜕𝑖𝑖𝑖𝑖 𝜕𝜕𝜕𝜕
 𝜖𝜖̇ 𝑖𝑖𝑖𝑖 = 𝜖𝜖̇ 𝑖𝑖𝑖𝑖𝑣𝑣𝑣𝑣𝑣𝑣 + 𝜖𝜖̇ 𝑖𝑖𝑖𝑖𝑒𝑒𝑒𝑒 + 𝜖𝜖̇ 𝑖𝑖𝑖𝑖𝑝𝑝𝑝𝑝 = 𝜏𝜏𝑖𝑖𝑖𝑖 + + 𝜆𝜆̇ ,
2𝜂𝜂𝑒𝑒𝑒𝑒 𝑒𝑒 2𝐺𝐺 𝜕𝜕𝜕𝜕 𝜕𝜕𝜕𝜕𝑖𝑖𝑖𝑖

where 𝐴𝐴𝐴 𝑖𝑖𝑖𝑖𝑣𝑣𝑣𝑣𝑣𝑣 , ̇𝜖𝜖𝑖𝑖𝑖𝑖𝑒𝑒𝑒𝑒 , and ̇𝜖𝜖𝑖𝑖𝑖𝑖𝑝𝑝𝑝𝑝 are the viscous, elastic, and plastic strain rates, respectively. ηeff is the effective viscosity, G
𝐴𝐴
the elastic shear modulus, τij the deviatoric stress tensor, t the𝐴𝐴time, 𝜆𝜆̇ is the plastic multiplier, Q the plastic flow
potential and σij = −P + τij the total stress. The effective viscosity ηeff is given by:
( ) 1 −1
1 −1𝑛 𝐸 + 𝑃𝑉
 𝜂𝑒𝑓 𝑓 = 𝐴 exp 𝑛
𝜖̇ 𝐼𝐼 ,
2 𝑛𝑅𝑇

HUSSON ET AL. 2 of 17
 15252027, 2022, 1, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2021GC010167 by Nat Prov Indonesia, Wiley Online Library on [11/04/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
Geochemistry, Geophysics, Geosystems 10.1029/2021GC010167

Table 1
Material Parameters Used in This Study
Unit ρ (kg m−3) Rheology A (Pa−n s−1) E (J mol−1) n V (m−3 mol−1) ϕ (°) C (MPa)
Upper mantle 3,300 Dry olivine (Hirth & Kohlstedt, 2004) 1.1 × 10 5
530 × 10 3
3.5 12 × 10 −6
20 30
Lower mantle 3,600 Dry olivine (Hirth & Kohlstedt, 2004) 1.1 × 105 530 × 103 3.5 12.6 × 10−6 20 30
Oceanic crust 3,300 Diabase (Mackwell et al., 1998) 8 485 × 10 3
4.7 0 0 5
Continental crust 2,800 Quartzite (Ranalli, 1995) 6.7 × 10−6 156 × 103 2.4 0 20 30
Note. Additional parameters, uniform across all materials are the elastic shear modulus G  =  5  ×  1010  Pa, thermal expansivity α  =  3  ×  10−5  K−1, heat capacity
Cp = 1,050 J K−1, and thermal conductivity κ = 3.0 W m1 K−1. The minimum and maximum viscosities are capped to 1 × 1019 and 1 × 1023 Pa s1, respectively.

creep, 𝐴𝐴𝐴 𝐼𝐼𝐼𝐼 is the second invariant of

where A is the exponential prefactor, n the stress exponent of the dislocation 𝐴𝐴
the viscous strain rate tensor, E, V are the activation energy and volume, respectively, P is the pressure, R is the
gas constant, and T is the temperature. Plasticity is modeled using the Drucker-Prager yield criterion given by:

 𝜏𝜏𝑌𝑌 = sin(𝜙𝜙)𝑃𝑃 + cos(𝜙𝜙)𝐶𝐶𝐶

where τY is the yield stress, ϕ the friction angle, and C the cohesion. Strain softening is taken into account by
linearly reducing both the friction angle and the cohesion of the material by a factor of 100 between 10% and 60%
of accumulated strain. Minimum cohesion is set to 0.01 MPa and maximum yielding stress to 900 MPa.

The modeled region is a 3-D box of dimensions 5,000 × 4,000 × 1,050 km with a resolution of 320 × 256 × 128 cells.
It comprises a 50 km thick air layer, lithospheric plates (oceanic and continental) embedded in a 660 km thick upper
mantle, and a 340 km thick lower mantle. Boundary conditions are no slip on the bottom, free slip on all sides and
free surface with open boundary on top. The temperature is set to 20°C at the top of the model and within the air lay-
er, and 1602°C at the bottom boundary. In order to keep the air layer at 20°C the conductivity in the air is artificially
set to κ = 100 W m−1 K−1 and the heat capacity is set to Cp = 106 J K−1. Initial temperature profiles within the plates
are prescribed according to the half-space cooling model (Turcotte & Schubert, 1982). For continental plates (Sun-
daland and Australia) the thickness of the boundary layer is set to 120 km and the half-space cooling age is 110 Ma.
For the lithosphere of the volcanic arc the thickness of the boundary layer is set to 100 km and the cooling age to
70 Ma. For oceanic plates, the thermal profile is set according to the half space cooling model with an 80 km thick
boundary layer and the cooling age is initialized as a function of the distance to the oceanic ridge, using a spreading
rate of 1.5 cm yr−1. In the upper mantle and uppermost lower mantle, the initial temperature gradient is set to 0.3 K
km−1. Continental and oceanic crusts are respectively, 35 and 15 km thick. Material properties are given in Table 1.
In practice, we performed a set of simulations that gradually incorporated complexity, starting from box models with
few rheological components, toward the final model that displays a complex arrangement of oceanic and continental
units as well as an island arc. In order to continuously have a weak subduction interface, the maximum strength of
the oceanic crust is capped by setting up the plastic cohesion to a low value (C = 5 MPa). Convergence between the
Indian-Australian and the Sundaland/Wallacea system is triggered by incorporating a pre-subducted segment of the
oceanic plate. The subsequent geodynamic evolution is self-consistent, entirely buoyancy driven.

2.2.  Bulk, Isostatic, and Dynamic Uplift Rates

Model simulations output the bulk uplift rates, accounting for both isostatic and dynamic components. At some
locations, both can be significant when strain rates and mantle flow efficiently distort the surface. We disentangle
the two components by approximating that the bulk uplift 𝐴𝐴 rate 𝐴𝐴𝑣𝑣bulk results from the joint effects of instantaneous
crustal isostatic𝐴𝐴uplift 𝐴𝐴𝑣𝑣iso and dynamic𝐴𝐴uplift 𝐴𝐴𝑣𝑣dyn :

 𝑢𝑢𝑣𝑣dyn = 𝑢𝑢𝑣𝑣bulk − 𝑢𝑢𝑣𝑣iso

zontal velocity at the surface ∇uh and crustal thickness S such that
We first approximate isostatic uplift rates by computing crustal thickening rates from the divergence of the hori-

1
 𝑢𝑢𝑣𝑣iso = 𝜌𝜌𝑚𝑚 −𝜌𝜌𝑐𝑐 ∇𝑢𝑢ℎ 𝑆𝑆
1+ 𝜌𝜌∗ −𝜌𝜌𝜌𝜌

HUSSON ET AL. 3 of 17
 15252027, 2022, 1, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2021GC010167 by Nat Prov Indonesia, Wiley Online Library on [11/04/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
Geochemistry, Geophysics, Geosystems 10.1029/2021GC010167

where ρm and ρc are the reference densities of the mantle and crust, while ρ∗ is the density of the loading material
(air, water, and sediments).

The model output dynamic topography and dynamic uplift rates are air-compensated. Water-compensated and
sediment-compensated deflections are respectively increased by a factor Δρa/Δρw and Δρa/Δρs, where Δρa, Δρw,
and Δρs are the density contrasts between the mantle and air, water, and sediments. Reference mantle density is
set in our model to 3,300 kg m−3, seawater has a density of 1,030 kg m−3 and we assume a density of 2,000 kg
m−3 for recent sediments. The “maritime continent” that we model includes air-, water- and sediment-loaded
topographies, and appropriate compensations need to be carefully chosen based on the geological record, when
comparing model outputs to observations. In the deep Banda and Weber basins, we thus account for water com-
pensation, while in the Sahul shelf, which is only flooded in the late stage, water- or sediment-compensations
have negligible impacts. Conversely in the Sunda shelf, the platform is largely flooded and covered by layers of
sediments and seawater of variable, yet limited thicknesses, indicating a combined compensation by air, water,
and sediments, which does not permit to precisely a priori elicit one of the three compensation mechanisms.
Nevertheless, we retain the intermediate situation (water-compensation) as a plausible situation between the two
end-members (air- and sediment-compensation) for the time-integrated dynamic deflection. Conversely, uplift
rates at present are the time derivative of the current elevation, and attest for the instantaneous compensation
process. In that case, sediment-compensated subsidence below sea level most appropriately reflects the current
sedimentation pattern in the Sunda shelf.

2.3.  Observed Uplift and Subsidence Rates

Mean Pleistocene uplift and subsidence rates are derived from unpublished and published estimates (see Table S1
and Figure S1 in Supporting Information S1). Most rates are estimated from the geomorphology of uplifting or
subsiding sequences of coral reefs. Some uplift rates are derived from the inversion of river profiles. In addition
in the Sunda shelf, the shallow stratigraphy of the platform (as seen on high-resolution seismic profiles) reveals a
series of incised sedimentary layers that are interpreted as fingerprints of the quickly oscillating sea level during
the Quaternary. Alternating periods of marine transgressions and regressions imprint the sedimentary cover by
episodes of sedimentary deposition followed by episodes of drainage. The number of these episodes indicates
the number of glacial cycles during which the platform was reached by marine transgressions. This number is
contingent on sea-level oscillations and bathymetry but also on subsidence rates (Sarr, Husson, et al., 2019). We
use this relationship to derive new estimates of subsidence rates from a survey in North East Java (Susilohadi
& Soeprapto, 2015), where 4–5 units (or 4 to 8 sub-units) have been identified. Using sea level reconstructions
(Waelbroeck et al., 2002), this allows bracketing subsidence rates in the area between 0.1 and 0.25 mm/yr (Figure
S2 in Supporting Information S1).

3. Results
3.1.  Australian Continent Disrupts Subduction Dynamics and Mantle Flow

We explore the regional subduction dynamics and surface morphotectonics using Cartesian three-dimensional
thermo-mechanical models. These simulations are “Earth-like,” in the sense that the setup attempts to objectively
integrate the specific geometries of the region without any a priori filtering of selected features. This choice
allows us to diagnose the impact of the arrival of the Australian continent into the subduction zone on regional
morphotectonics in both their temporal and spatial frameworks since the early Oligocene. More specifically, the
close resemblance to real Earth facilitates the comparison of the model outcomes to our set of observations (see
initial setup Figure 2). Our scenario starts prior to the collision, at a time when the Indo-Australian lithosphere
subducted North of Papua during the early Oligocene. At that stage, the existence of the Banda embayment
prior to collision, as well as an island arc and alternating landmasses and oceanic basins in the upper plate best
account for the pre-collisional conditions (Harris, 2011). We only present our final simulation, extracted from a
set of experiments on its ability to reproduce kinematic and physiographic observations (Hall, 2017; Spakman &
Hall, 2010). Complete modeling approach, alternative simulations, and animations are presented in Supplemen-
tary Materials (Figures S3–S8 in Supporting Information S1).

Our simulation unravels the joint evolution of the regional physiography and mantle flow (Figure 3). The first
event to disrupt the northward journey of the Australian continent is the inception of the Papuan promontory

HUSSON ET AL. 4 of 17
 15252027, 2022, 1, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2021GC010167 by Nat Prov Indonesia, Wiley Online Library on [11/04/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
Geochemistry, Geophysics, Geosystems 10.1029/2021GC010167

5LGJH )5((6
8 5)$&(
5LGJH
,QGLDQ
SODWH
D
UF
6XQGDODQG 1RUWK
QLF
FD
$XVWUDOLD RO

)5((6/,3
Y
6XODZHVL NP
3DSXD
8SS
HUP 3KLOLSSLQH
DQWO SODWH
H
/RZ
HUP NP
DQW ,QLWLDOO\VXEGXFWHG 
OH VHJPHQW ƒ&

  & 

)5((6/,3
& & 
 
NP /
/
ƒ&
 /
ƒ
&   

80 80 80

'HSWK
NP
&RQWLQHQWDO 9ROFDQLF 2FHDQLF
DUF

Figure 2.  Initial setup of the model. The model is entirely buoyancy driven, forces arising from mantle thermal convection,
from which ridge push and slab pull arise. Initial configuration is based on reconstructions (Ely & Sandiford, 2010;
Hall, 2012, 2017; Hall & Spakman, 2015; Harris, 2011; Spakman and Hall., 2010). The model accounts for 3 different
vertical stratifications (as depicted on the columns). The thermo-mechanical structure of the Indian and Philippine oceanic
plates is set by half-space cooling away from the ridge. The shapes of continental units are defined according to their current
and reconstructed geometries prior to collision. The small continental units North of Australia represent the envelope of the
loosely constrained assemblage of small-size continental units prior to collision with Papua; they play a negligible role in
the modeled dynamics. The volcanic arc has a more ductile behavior. See Section 2 for technical description of the thermal/
mechanical aspects and Table 1 for material properties.

into the subduction zone. Until that stage, the mantle flow is essentially poloidal (in the vertical plane, mostly
normal to the trench), accompanying the descent of the oceanic lithosphere in the upper mantle. Soon after, the
slab breaks off and collision proceeds in Papua. In the upper mantle, a toroidal (in the horizontal plane, parallel
to the surface) component of mantle flow expands during the breakoff of the slab segment underneath Papua.
Mantle flow accelerates and conveys the mantle from the sub-Australian (or Indian) reservoir into the Eurasian
reservoir. This revised scheme of mantle circulation facilitates rollback of the Banda slab as the trench winds
around the Bird's Head, West Papua and eventually retreats and nests in the Banda embayment (Figures 3 and 4).
This scenario conforms to previous kinematic reconstructions based on surface geological observations and seis-
mic tomography (Milsom,  2001; Spakman and Hall.,  2010), which depict a continuous slab that increases in
convexity as it molds in the Banda embayment. At first sight, our model predictions validate the mechanics of
this conceptual model. However, the situation is more complex at depth, as such trench runaway requires extreme
distortions of the subducting plate. The once more linear Miocene trench running from Java to Papua is now
bent enough to bear a unique slab that subducts with opposite vergences within only a few hundred kilometers
along the Banda arc (Figure 1 and Figure S1 in Supporting Information S1). While breakoff underneath Papua
in the North eased the retreat of the trench and the subsequent opening of the Banda back-arc basin, the situa-
tion differs in the South. After 26.8 My (forward model time, 0 My being the beginning of the simulation), the
setting resembles that of present day (0 Ma actual time, Figures 3e and 3f). At that stage, not only the Australian
continent has already docked onto the trench, which causes vertical tension in the subducted lithosphere (Ely &
Sandiford, 2010; Royden & Husson, 2009), but eastward slab rollback also requires some longitudinal stretch-
ing of the slab, in an East-West direction, mirroring the trench elongation along the margin of North Australia
by at least a factor of two (see for instance reconstructions by Spakman and Hall (2010), and Figure 4). Both
geometrical constraints impose considerable distortions of the lithosphere. At the surface, this is consistent with

HUSSON ET AL. 5 of 17
 15252027, 2022, 1, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2021GC010167 by Nat Prov Indonesia, Wiley Online Library on [11/04/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
Geochemistry, Geophysics, Geosystems 10.1029/2021GC010167

Figure 3.  Model simulations of the Indo-Australian subduction zone. Time advances forward from top to bottom (model time is given in My since the beginning of
the experiment, approximate equivalent Earth time is given in Ma). (a, c, e, and g) Westward views of the continental topography and Indian oceanic lithosphere (blue
surface and thermal structure, eastern side view). (b, d, f, and h) Close ups of the Banda subduction zone, southward views. Indian oceanic lithosphere (blue surface)
and continental units in brown (Australia) and white contours (Wallacea islands). Streamlines show mantle flow. (i, j, k, and l) Power spectrum of the strain rate in
horizontal slices (normalized to the highest power). See animations in Supporting Information S1.

HUSSON ET AL. 6 of 17
 15252027, 2022, 1, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2021GC010167 by Nat Prov Indonesia, Wiley Online Library on [11/04/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
Geochemistry, Geophysics, Geosystems 10.1029/2021GC010167

m
−17 −16 −15 −14 −13 −6 −3 0 3 6
strain rate (log10 s-1) elevation (km)
A) [2] B)
1000
[1] [1]

−1000
N 4.4 My [24.3 Ma]

C) D)
1000
[5] [2] [5]
[3] [4] [4]
[3]
0

[6]
−1000
15.2 My [13.5 Ma]

E) F)
1000

Sundaland
[2]
[3] [8] [9] [8] [9]
[3]
0
[10]

Weber
Sahul
[6]
−1000
26.9 My [0 Ma]

G) H)
1000
[11]
[2] [2]

[6]
−1000
31.2 My [+4.3 Ma]
−2000 −1000 0 1000 2000−2000 −1000 0 1000 2000

Figure 4.  Modeled strain rates and topography. Time advances forward from top to bottom (model time in My, approximate
equivalent Earth time in Ma). (a, c, e, and g) Surface strain rates, continental units in solid black contours. (b, d, f, and h)
Elevation and interpretative tectonic framework. [1] arc-continent collision, N Papua [2] Southward tilt, Sundaland [3]
incipient compression, Java [4] back-arc thrust, Wetar-Flores [5] Central range orogeny, Papua [6] Northward tilt, Australia
[8] back-arc propagation to Sundaland [9] slab rollback and back-arc opening, Banda Sea [10] Savu Sea opening, and
compression in Sumba [11] subduction polarity reversal.

geological observations of lithospheric delamination and with the back-arc basins that spread in the Banda Sea as
the trench retreats and stretch at high strain rates (Spakman & Hall, 2010 and Figure 4). At depth, slab pull tends
to neck the slab, while its eastward rollback extensively stretches it. Whether such a scenario is mechanically vi-
able without disrupting the slab underneath Timor or Wetar is disputed. Seismic models yield contrasted results:
while seismicity catalogs (Das, 2004; Ely & Sandiford, 2010) as well as some seismic tomography models (Widi-
yantoro & van der Hilst, 1997; Widiyantoro et al., 2011) tend to confirm slab detachment in this region, other
tomographic models are less clear (Hall & Spakman, 2015) or even indicate no detachment (Harris et al., 2020;
Porritt et al., 2016; Zenonos et al., 2019). Within the range of our parametric search, our simulations show little

HUSSON ET AL. 7 of 17
 15252027, 2022, 1, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2021GC010167 by Nat Prov Indonesia, Wiley Online Library on [11/04/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
Geochemistry, Geophysics, Geosystems 10.1029/2021GC010167

alternative to slab tearing and breakoff. At depth, the main consequence of this event is to gradually reinforce the
vigor of the toroidal cell, which in turn boosts the poloidal flow (Figure 3f). Likewise, our model indicates that
the slab enters the lower mantle only in the easternmost side, when the Banda slab becomes isolated from the
rest of the Indo-Australian slab, and ponds at the base of the upper mantle elsewhere (Figure S9 in Supporting
Information S1). This is at odds with most seismic tomography models that indicate that the slab has penetrated
the lower mantle all along its strike, indicating that our model dynamics are late with respect to actual penetration
of the slab in the lower mantle.

This model reveals the links between complex surface tectonics (slab rollback, trench rotation by more than 90°
and encroachment in the Banda embayment) and a simple pattern of mantle flow, only composed of a coun-
terclockwise toroidal flow cell and a poloidal cell. This behavior is captured by the power spectrum of the
strain rate (Figures 3i–3l). During early simulation stages, the spectrum is roughly depth-independent, the low
frequencies being present at all depths mirroring the concordance between surface flow and surface tectonics,
in a quasi-steady state (Figure 3i). After the collision in Papua and subsequent slab breakoff, high frequencies
are excited near the surface, but gradually vanish with depth (Figure 3j). That depth dependence of the power
spectrum remains from that stage onward, and is further reinforced once the subduction wraps around the Banda
embayment and the slab detaches underneath Timor-Wetar (Figures 3k and 3l). This temporal evolution of the
power spectrum is conditioned by the structure of the mantle; in the early stages, the Tethyan subduction dom-
inates the poloidal flow; complexity complexity arises once smaller-size structures appear. Overall during tran-
sient subduction dynamics, the deep pattern of deformation remains restricted to long-wavelengths, irrespective
of the high degree of tectonic complexity that arises in the shallowest layers. The short wavelengths of tectonic
activity instead originate from the disruption of the complex structure of the uppermost layers. Prior to the colli-
sion of the Australian continent, the surface lithosphere was heterogeneously structured: continental Sundaland
in the West and the complex assemblage of continental and oceanic units of Wallacea in the East. The power
spectrum during and after the collision of Australia indicates that this small-scale tessellation does not prevent
the long-wavelength pattern of surface deformation to resemble that of the underlying mantle flow, but instead
superimposes shorter wavelengths that are the cause of the observed short-lived and fast morphotectonic activity
that dominate the region.

3.2.  Back-Arc Thrust Propagates and Builds Up Toward Subduction Reversal

Intense tectonic activity in and around the subduction arc is another fundamental expression of regional geody-
namics. Mega-earthquakes occur on the subduction megathrust as in other subduction zones, but the seismotec-
tonic activity of the back-arc domain, underscored by large magnitude earthquakes like the 1992 Mw 7.8 Flores
and 2004 Mw 7.5 Alor earthquakes, is more unusual. The epicenters of these destructive earthquakes punctuate
the Flores and Wetar back-arc thrusts that developed during the Australian continent impingement into the sub-
duction zone (Figure 1) (Koulali et al., 2016; Silver et al., 1983). The 2018 Mw 6.9 Lombok earthquake sequence
highlighted the westward propagation of the thrust from the Flores and Wetar back-arc thrusts. The Lombok
events occurred on blind thrusts (Yang et al., 2020), indicative of an early stage of faulting, which invigorates
the question of the propagation of the fault into continental Sundaland. Mounting tectonic and geodetic evidence
suggest that the Baribis and Kendeng thrusts of Java island (Figure 1) stem from the Flores and Wetar back-arc
thrusts (Koulali et al., 2016) and form a continuous Java back-arc thrust. Are these discontinuous faults merging
into a continuous plate scale, an even more forceful thrust that would mark the inception of a subduction reversal?

This scenario is supported by the pattern of surface deformation in our simulation (Figure 4). Strain rates peak
in the arc until collision occurs in Papua (Figures 4a and 4b), after which contraction initiates in the back-arc,
approximately in the Wetar Flores thrusts (Figures 4c and 4d) and eventually gently propagates westward to enter
continental Sundaland at a stage that best compares to the present (Figures 4e and 4f). On the eastern side, the

drop more than tenfold when entering Sundaland (from ∼5 × 10−15 s−1 to ∼5 × 10−16 s−1, Figure 4e). At present,
modeled tectonic activity reproduces the severe deformation of the Flores and Wetar back-arc thrusts, strain rates

observed and modeled surface deformation fields (Figure 5) are in good agreement. Model velocities after the
collision of the Australian continent decrease to rates that are approximately twice smaller than current geodetic
rates (see also Figure S1 in Supporting Information S1), but show a remarkably consistent pattern of deformation.
Both geodetic and modeled velocities strongly decrease westward along the island arc (with respect to stable
Sundaland) with a reduction in the northward velocity at the edge of continental Sundaland. This result is further

HUSSON ET AL. 8 of 17
 15252027, 2022, 1, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2021GC010167 by Nat Prov Indonesia, Wiley Online Library on [11/04/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
Geochemistry, Geophysics, Geosystems 10.1029/2021GC010167

−2000 −1000 0 1000 2000 mm/yr

25
a) absolute model velocities

1000

20

15

−1000

10
500
b) model velocities relative to stable Sundaland

0
−1000 −500 0 500

−5˚ c) geodetic velocities relative to stable Sundaland

50 mm/yr
−10˚

110˚ 115˚ 120˚ 125˚

Figure 5.  Observed and modeled surface velocities. (a) Global velocities in the model reference frame; (b) zoomed in
velocities, with respect to stable Sundaland (at coordinate −500/500 km). (c) Geodetic velocities relative to Sundaland
(Koulali et al., 2016) at an approximately comparable geographical scale than in panel (b).

reinforced by the similarity between strain rates (Figure 4e) and geodetically derived strain rates (Figure S10 in
Supporting Information S1).

Structural observations, indicative of the bulk deformation, are also compatible: while shortening in the eastern
part of the arc (including Sumba, Flores, and Timor islands) is characterized by important nappe stacking, back-
arc thrusting and fast uplift, Java island at the southern rim of Sundaland is affected by the modest Plio-Quater-
nary fold and thrust belt in the Kendeng hills of East Java, and the even more modest thrusts of the Baribis fault
zone in the West (Clements et al., 2009; Simandjuntak & Barber, 1996). The model outcome, when compared to
these tectonic and kinematic observations, suggests that the Wetar-Flores back-arc thrusts are currently merging
and prograding westward, entering Java island where the Java back-arc thrust would be the active propagator
of the fault. As it is common practice in atmospheric science, numerical models permit to be predictive by

HUSSON ET AL. 9 of 17
 15252027, 2022, 1, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2021GC010167 by Nat Prov Indonesia, Wiley Online Library on [11/04/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
Geochemistry, Geophysics, Geosystems 10.1029/2021GC010167

deduction. Here, predictions for future times from our simulations (Figures 4g and 4h) suggest that the deforma-
tion in the back-arc will eventually evolve into a major and continuous south verging thrust that will ultimately
become the most active structure of the region. During the process, subduction will initiate after small oceanic
basins in Wallacea begin to subduct along that fault, leading to subduction polarity reversal. These small oceanic
basins, which spread southward during Banda slab rollback, will disappear while the subduction reverses, with a
northward-retreating trench.

3.3.  Dynamic Topography Drowns Down the Weber Basin to 7,000 m

Banda subduction dynamics quickly evolved throughout the course of the Australian convergence. In particular,
as the trench hides away in the Banda embayment, its curvature increases, and it now tightly confines the Weber
and Banda basins. As a young oceanic basin (∼3 Ma or less; Hinschberger et al., 2005), the Weber basin is ex-
pected to stand at shallower depths. Likewise, the adjacent North and South Banda Basin, which opened earlier
during the Miocene, lie at considerable depths (4,000–5,000 m). The anomalous depth of the Weber Deep has
been tentatively linked to tectonics (Pownall et al., 2016), but it is unclear how even an overly thinned lithosphere
would reach much greater depths than that of back-arc basins or oceanic ridges, for instance. If the comparison
with ridges holds, the residual topography (i.e., the component of the topography that departs from isostatic equi-
librium) of such a young basin is undoubtedly in excess of several kilometers, for which explanations are lacking.

During the retreat and curvature of subduction trenches, the convexity of slabs also increases as they sweep into
the mantle (Loiselet et al., 2009). The mass anomalies of dense slabs get concurrently more and more focused
underneath their back-arc basins, similarly to mantle downwellings that accompany slab descent. The vertical
traction exerted by the viscous flow underneath the Earth's surface accordingly increases and warps the surface
downward above the mass anomaly (Husson, 2006). As the slab breaks off, the mass excess gets deeper, vertical
tractions at the surface decrease, and the topography returns to isostatic equilibrium. This scenario applies very
well to the Weber and Banda basins for which our model permits the evaluation of the magnitude of this dynamic
topography (Figure 6). As the slab enters the Banda embayment (25.5 My model time), it gets tightly squeezed
between the northern arc (equivalent to Seram) and the northward migrating Australian margin; the predicted
dynamic deflection reaches its maximal value (Figures 6b and 6c). Retreat proceeds and dynamic topography
decreases because the slab flattens as it lies at the bottom of the upper mantle, and because the slab tears to the
North of Australia. The situation that most resembles the present in the model (when the trench meets the eastern

geodynamic evolution (by ∼1.8 My, 28.7 My instead of 26.9 My, see above). At that stage, the predicted dynamic
end of the Banda embayment) is slightly delayed with respect to the reference stage that regionally matches the

∼4 km in the West to ∼7 km in the East (Figure 6c). Our model suggests that this ∼3 km deep deflection could
deflection above the subduction (from the trench to the Western part of the basin) warps down the basin from

the topography peaked to ∼4 km and drowned the basin to ∼8 km before declining. As the slab breaks off, the
have been even larger during the earlier stages (25.5 and 26.9 My, Figure 6c) when the dynamic deflection of

amplitude decreases and topography will be almost back to isostasy at 31.2 My (model time).

These simulations suggest that the Weber Deep and more generally the Banda basin largely owe their exceptional
depths to dynamic topography. This deflection is likely the most important one on Earth. Other back-arc basins
lie at considerable depths in similar settings, but none is as deep as the Weber Deep. In continental domains, the
basement below the Focsani (Carpathians) and Ganga (Himalaya) foreland basins is subsided down to 13 and
8 km depths by dynamic subsidence (Husson et al., 2014; Şengül-Uluocak et al., 2019), but they only surpass the
Weber Deep because the sediment load reinforces the deflections. Our model results also indicate that the Banda
arc and affiliated basins are currently uplifting as dynamic topography decreases. This result partly explains the
ubiquitous uplifted Pleistocene coral reef sequences that festoon the islands of the Banda arc (Figure  1). Al-
though some are uplifted due to tectonic contraction (for instance, Sumba island, Authemayou et al., 2018, Miller
et al., 2021, or SE Sulawesi, Pedoja et al., 2018), the background uplift of the region is induced by the relaxation
of the dynamic topography. This motion is well localized around the basin, while the Australian margin converse-
ly continues to subside at present.

HUSSON ET AL. 10 of 17
 15252027, 2022, 1, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2021GC010167 by Nat Prov Indonesia, Wiley Online Library on [11/04/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
Geochemistry, Geophysics, Geosystems 10.1029/2021GC010167

0˚
B
Se
31.2 My
-4
Sul Bu
−4˚ nB P

z (km)
W
sB -6
−8˚
F T 25.5 My
28.7 My
Su Ar

1000 km
26.9 My -8
−12˚ A) A C)
120˚ 124˚ 128˚ 132˚ 136˚ 140˚
-600 -400 -200 0
B) distance (km)
500 500

0 0

25.5 My 28.7 My
−500 −500
500 500

0 0

26.9 My 31.2 My
−500 −500
0 1000 2000 0 1000 2000
m
−6 −3 0 3 6
elevation

Figure 6.  Dynamic topography in the Banda arc and Weber Deep. (a) Bathymetry, Sul, Sulawesi; Bu, Buru; Se, Seram; B,
Bird's Head; P, Papua; A, Australia; T, Timor; F, Flores; Su, Sumba; Ar, Arafura sea; nB, sB, North and South Banda basins;
W, Weber Deep. (b) Modeled topography during Banda arc formation, (c) Profiles of predicted dynamic topography (water-
compensated—Section 2), truncated in the East at the location of the subduction trench to show the upper plate only. Location
Figure 3f.

3.4.  Dynamic Subsidence in the Sunda and Sahul Shelves

The Banda basin and neighboring islands broadly define the mostly uplifting region of Wallacea, while at its
South and West, it is surrounded by the flooded continental platforms of Sahul and Sundaland (Figure 1). They
stand above a “subduction graveyard,” where considerable amounts of oceanic lithosphere have been drawn
down in the mantle within multiple subduction zones. These flooded continents have long been assigned to the
dynamic topography associated with whole-mantle circulation at present-day (Conrad & Husson,  2009; Gur-
nis, 1993; Ricard et al., 1993). The vertical traction exerted by the downwelling mantle warps the surface of the
Earth downward enough to expose continents to marine transgressions. More specifically, subsidence in Sahul is
thought to have increased over time as the continent approaches the subduction and overrides the Indo-Australian
subducting plate (DiCaprio et al., 2011); long term subsidence of Sundaland is proposed to result from slab av-
alanching in the lower mantle (T. Yang et al., 2016). At present-day, both continental platforms subside at a few
tens of mm/yr (Figure 1), about 10 times faster than the slow subsidence predicted from subduction dynamics
models (DiCaprio et al., 2011; T. Yang et al., 2016). In addition, hectic movements of uplift and subsidence are
reported in Sahul since Late Miocene (Gurnis et al., 2020).

by ∼0.4‰ on average. Our models predict extremely comparable values of dynamic tilt (Figure 7b). In addition,
The topography of North Australia, the “tilting continent” (Sandiford, 2007) dips toward the trench (Figure 7a)

HUSSON ET AL. 11 of 17
 15252027, 2022, 1, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2021GC010167 by Nat Prov Indonesia, Wiley Online Library on [11/04/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
Geochemistry, Geophysics, Geosystems 10.1029/2021GC010167

Figure 7.  Australia, “the tilting continent” (Sandiford, 2007). (a) Swath topographic profile across the northern half of
Australia (solid curve: average, location in inset map). (b) Predicted topography (water-compensated—gray dashed line—and
sediment-compensated—gray dotted line—platform bathymetries are shown for completeness). (c) Uplift rates at three model
time steps. Diamonds show observations of Late Pleistocene subsidence rates (S, Scott reef; A, Adele reef; C, Cockatoo
Island, see Supporting Information S1).

the topographic profiles are very similar, showing maximal slopes between 500 and 1,000 km from the trench
(Figures 7a and 7b, particularly after 31.2 My model time; again, slightly delayed with respect to the regional
“reference stage” at 26.9 My). Subsidence rates also vary over time, and at present-day, they only slightly ex-
ceed estimates of Pleistocene subsidence based on observations (Figure 7c). Our simulation confirms that the
northward tilt of Australia is caused by the dynamic deflection of the descending slab. Revived subsidence at
present-day happens in our simulation during slab tearing. The toroidal flow strengthens from this slab hole
(online animation S2) which facilitates the descent of the slab and its retreat underneath NE Australia, relocating
dynamic subsidence from the Banda basin toward Australia. This process makes subsidence rates quickly vary in
time, but also spatially in North Australia (Figure 8). This quickly evolving pattern in turn explains the fast and
“reversible subsidence” observed in the NE margin of Australia (Gurnis et al., 2020).

ure 9a) reveals that the flat continental platform is tilted southward by ∼0.6% and that sedimentary layers level the
The Sunda shelf, in the North of the Indonesian arc, lies at shallower depths than 100 m. Seismic reflection (Fig-

surface. Our model predicts comparable tilting due to dynamic topography (Figure 9b), which permits to assign
a dynamic origin to the tilt of Sundaland. The amplitude of the deflection is larger than that of Sahul because the
Sunda shelf lies on the upper plate and is therefore located above the dense, northward dipping subducting mass
anomaly of the slab. The sedimentary infill is unconformably disposed and shows a northward progradation of

HUSSON ET AL. 12 of 17
 15252027, 2022, 1, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2021GC010167 by Nat Prov Indonesia, Wiley Online Library on [11/04/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
Geochemistry, Geophysics, Geosystems 10.1029/2021GC010167

mm/yr onlapping sediments; the sedimentary layers are thickening southward in a

−0.4 −0.2 0.0 0.2 0.4 fan shape. While the unconformity indicates that the tilt precedes to a large
uplift rate extent the deposition of the sediments since Miocene (Smyth et al., 2008), the
deposition pattern is that of a broad trenchward-increasing subsiding basin.
At present-day, subsidence rates amount to 0.1–0.3 mm/yr (Figure 1). Our
1000
simulation shows that most of the deflection was acquired during early stages
of subduction (Figure  4d), but that subsidence resumed after the northern
0 margin of Australia collided with the arc (Figure  4e), yielding subsidence
rates in the range of the geological estimates (Figures 7c, 9b, and 9c). During
these events, dynamic subsidence rates in Sundaland varied through time:
−1000 during the initial stages of oceanic subduction, poloidal flow over the large
A) bulk
Indo-Australian subducting plate deflected the upper plate. Following the
collision of Australia, both slab detachment underneath Papua and slab tear-
1000 ing underneath Timor fostered the development of a toroidal cell over a now
narrower Indian oceanic slab (Figures 3d–3f). Subduction and mantle flow
rates increased, reinforcing the dynamic warping of the upper plate originally
0
acquired during pre-collisional stages, thereby causing the subsidence that
currently prevails.
−1000
B) isostatic
4. Discussion
1000 The setup of our model was designed to resemble Early to Mid-Miocene
reconstructions at the starting stage, in a need to facilitate a close match to an
0 extensive set of morphotectonic observations. One of the main outcomes of
this simulation is to reveal how surface complexity can be induced by a much
simpler mantle circulation scheme at depth. While the Banda-Java subduc-
−1000 tion zone displays an apparent tumult of various morphotectonic expressions
C) dynamic, air comp. at the surface, the underlying mantle flow is the mere adjunction of a coun-
terclockwise toroidal cell (starting when the slab initially detached under-
1000 neath Papua) to the dominantly poloidal flow that exists as long as subduc-
tion proceeds. The slowly evolving pattern of mantle flow in the aftermath
of the collision of the Australian continent with Papua and with the Banda
0
arc remains simply structured. It nevertheless remodels the physiography at
various scales: jerks of dynamic subsidence are excited by slower subduction
−1000 dynamics. If morphotectonic time and space scales do not seem to match
D) dynamic, sedim. comp. between the surface and the mantle, it is because mantle flow enters a tran-
−2000 −1000 0 1000 2000 sient phase that relocates stresses underneath the already structured surface
lithosphere. Surface tectonics indicate that horizontal stresses quickly vary,
Figure 8.  Modeled uplift rates after 26.9 My (model time). (a) Bulk, (b) and subsidence rates reveal that changes in dynamic topography can be very
isostatic, (c) dynamic, assuming air-loading, and (d) dynamic, assuming fast. The “surface expression of mantle convection” may also occurs at faster
sediment-loading (see Section 2).
rates and shorter space scales than the very mechanism it originates from
(Figure 3, Coltice et al., 2019; Mallard et al., 2016). While the long-wave-
length pattern of deformation adjusts to that of the underlying mantle flow,
the short-wavelength deformation is the transient response of the tessellated, heterogeneous lithosphere. Such
relationships are likely valid at other plate boundaries that the SE Asian region illustrates more than elsewhere.

From a regional standpoint, our model mechanically validates the earlier conceptual scenario proposing that the
Indo-Australian subducting slab rolled back and curled up around Papua before molding the Banda embayment
(Milsom, 2001; Spakman & Hall, 2010). However, in our simulation, the fate of the Banda slab differs. It is inex-
orably torn and eventually detached underneath Timor-Wetar before it reaches the far East Banda arc. Although
seismic models are inconsistent regarding the existence of a slab window in this region (Ely & Sandiford, 2010;
Hall & Spakman, 2015; Harris et al., 2020; Porritt et al., 2016; Widiyantoro & van der Hilst, 1997; Widiyantoro
et al., 2011; Zenonos et al., 2019), this result is difficult to reconcile with the most recent seismic studies that state
that slab integrity is preserved. If the rollback model holds, geometrical consistency (or, put differently, mass

HUSSON ET AL. 13 of 17
 15252027, 2022, 1, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2021GC010167 by Nat Prov Indonesia, Wiley Online Library on [11/04/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
Geochemistry, Geophysics, Geosystems 10.1029/2021GC010167

Plio - IV
1 km
V

Miocene
dip ~0.6%

Paleogene
1 s twt
20 km
km A)
Java island Sunda shelf

500
0
elevation (m)

-1000

-2000 .6%
~0
dip
B)
-3000
0.1
uplift rate (mm/yr)

-0.1
S

Ba
-0.2
B
C)
S 0 500
distance from tectonic front (km) N
Figure 9.  “Subsiding Sundaland” (Sarr, Husson, et al., 2019). (a) Seismic profile across the Sunda shelf (courtesy Pusdatin
ESDM Indonesia, location as a red bar in inset map; conversion from two-way times to meters assuming a uniform wave
speed of 1,800 m/s in sediments); V, Volcanic rocks (undated, but still active). (b) Predicted topography (model-equivalent
location in inset map), for three surface loading hypotheses: air-compensated (orange), water-compensated (green), sediment-
compensated (yellow). Scales are different for panels (a and b). (c) Predicted uplift rates. Bulk (gray dashed line) and dynamic
subsidence (solid curve, see Methods). Diamonds show observations of Late Pleistocene subsidence rates (B, Belitung; Ba,
Bawean; S, Singapore, see Section 2).

conservation) invites reappraising seismic tomography models. One limitation of the model could conversely
be turned into a means to lift the paradox between the dynamic and seismological standpoints: while we used
state-of-the art rheologies in our models, we acknowledge that the mechanical behavior is only known to a certain
extent. Less viscous oceanic lithospheres or a more dynamic rheology at continental margins could possibly hold
the Banda slab longer while allowing for its efficient stretching without tearing. Likewise, decoupling of the man-
tle lithosphere from the crust could facilitate the subduction of the mantle lithosphere (Capitanio et al., 2010).
However, given the current knowledge, such fine-tuning of the models would be fraught with unsupported spec-
ulations. Conversely, the successful and scrupulous model predictions of several morphotectonic features that
we deem diagnostic of the processes at play come in support of the dynamic scenario. The excellent agreement
between observations and models indicates that, following slab rollback in the Banda embayment, subsequent
rerouting of mantle flow during slab detachment underneath Timor-Wetar subsided the Weber basin to abyssal
depths and reinvigorated warping and flooding in the Sundaland and Sahul continental platforms until present.
Concurrently, compression increases in the arc, causing the Wetar-Flores back-arc thrust to quickly propagate
westward into continental Sundaland while subduction polarity is on the brink to reverse. Other data, such as the

HUSSON ET AL. 14 of 17
 15252027, 2022, 1, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2021GC010167 by Nat Prov Indonesia, Wiley Online Library on [11/04/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
Geochemistry, Geophysics, Geosystems 10.1029/2021GC010167

nature of the volcanism or a comprehensive analysis of seismic anisotropy could help settle the long-lived debate
regarding existence of a detached slab underneath eastern Indonesia.

A natural follow-up of this study would be to address the interaction of the solid Earth with the soft Earth. The
quickly evolving Late Cenozoic physiography of Southeast Asia interacts with the hydro-, atmo-, and bio-spheres
of the Earth likely more dynamically than anywhere else (Molnar & Cronin, 2015; Salles et al., 2021; Sarr, Sep-
ulchre, & Husson, 2019). High-resolution geodynamic models like the current one permit to envision a consistent
integration of all Earth's spheres at once.

5. Conclusion
At the eastern end of the Tethys, the subduction of the Indo-Australian plate underneath Southeast Asia ideally
permits to investigate the expression of a single, long-wavelength geodynamic event on the morphotectonic activ-
ity at various time and space scales. The current “Earth-like” numerical simulation of the subduction zone shows
that: (a) Docking of the Australian continent onto the subduction arc disrupted the mostly cylindrical mantle flow
by adding a long-wavelength toroidal flow cell to the existing long-wavelength poloidal flow cell. (b) Mantle
coupling to the previously structured lithosphere causes many short-lived and small-scale tectonic perturbations,
as exemplified by slab rollback and curling into the Banda embayment, back-arc thrusting and future subduc-
tion reversal to the North of Timor-Wetar, or lateral propagation of the back-arc thrust onto Java island. (c) The
Weber superdeep basin owes its exceptional bathymetry to the concentration of mass anomalies related to the
subducting slab around the Banda arc. This is the largest amplitude of dynamic topography on Earth. (d) Dynamic
subsidence in the North of Australia increases as it gets closer and overrides its own oceanic slab. (e) Dynamic
topography over the Indian slab floods the Sunda shelf. The revised subduction dynamics following the collision
of the Australian continent fosters dynamic subsidence at present-day. The magnitudes, rates, and velocities of
all these morphotectonic events are not only high, compared to other subduction zones, but are also caused by a
unique geodynamic event.

Data Availability Statement

Seismic data were obtained courtesy of Pusdatin ESDM, Indonesia (MoU LIPI, Indonesia). Numerical code LaM-
EM (9) is available at https://bitbucket.org/bkaus/lamem. All maps were done with GMT (Wessel et al., 2019).
All other data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary
Materials.

Acknowledgments References
This research was partially funded by
National Geographic Society, Lembaga Authemayou, C., Brocard, G., Delcaillau, B., Molliex, S., Pedoja, K., Husson, L., et al. (2018). Unraveling the roles of asymmetric uplift, normal
Pengelola Dana Pendidikan, Nusantara faulting and groundwater flow to drainage rearrangement in an emerging karstic landscape. Earth Surface Processes and Landforms, 43,
Hubert Curien partnership, and Euro- 1885–1898. https://doi.org/10.1002/esp.4363
pean Research Council. The authors are Capitanio, F. A., Morra, G., Goes, S., Weinberg, R., & Moresi, L. (2010). India-Asia convergence driven by the subduction of the Greater Indian
grateful to O. H. Göğüş and to our editors continent. Nature Geoscience, 3, 136–139. https://doi.org/10.1038/ngeo725
for their constructive reviews. No other Clements, B., Hall, R., Smyth, H. R., & Cottam, M. A. (2009). Thrusting of a volcanic arc: A new structural model for Java. Petroleum Geosci-
unpublished data were used or created for ence, 15, 159–174. https://doi.org/10.1144/1354-079309-831
this research. Coltice, N., Husson, L., Faccenna, C., & Arnould, M. (2019). What drives tectonic plates? Science Advances, 5, eaax4295. https://doi.org/10.1126/
sciadv.aax4295
Conrad, C. P., & Husson, L. (2009). Influence of dynamic topography on sea level and its rate of change. Lithosphere, 1, 110–120. https://doi.
org/10.1130/l32.1
Das, S. (2004). Seismicity gaps and the shape of the seismic zone in the Banda Sea region from relocated hypocentres. Journal of Geophysical
Research, 109, B12303. https://doi.org/10.1029/2004jb003192
Davies, D. R., Valentine, A. P., Kramer, S. C., Rawlinson, N., Hoggard, M. J., Eakin, C. M., & Wilson, C. R. (2019). Earth's multi-scale topo-
graphic response to global mantle flow. Nature Geoscience, 12, 845–850. https://doi.org/10.1038/s41561-019-0441-4
DiCaprio, L., Gurnis, M., Müller, R. D., & Tan, E. (2011). Mantle dynamics of continentwide Cenozoic subsidence and tilting of Australia.
Lithosphere, 3, 311–316. https://doi.org/10.1130/l140.1
Ely, K. S., & Sandiford, M. (2010). Seismic response to slab rupture and variation in lithospheric structure beneath the Savu Sea, Indonesia.
Tectonophysics, 483, 112–124. https://doi.org/10.1016/j.tecto.2009.08.027
Gurnis, M. (1993). Phanerozoic marine inundation of continents driven by dynamic topography above subducting slabs. Nature, 364, 589–593.
https://doi.org/10.1038/364589a0
Gurnis, M., Kominz, M., & Gallagher, S. J. (2020). Reversible subsidence on the North West Shelf of Australia. Earth and Planetary Science
Letters, 534, 116070. https://doi.org/10.1016/j.epsl.2020.116070
Hall, R. (2012). Late Jurassic-Cenozoic reconstructions of the Indonesian region and the Indian Ocean. Tectonophysics, 570–571, 1–41.
https://doi.org/10.1016/j.tecto.2012.04.021

HUSSON ET AL. 15 of 17
 15252027, 2022, 1, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2021GC010167 by Nat Prov Indonesia, Wiley Online Library on [11/04/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
Geochemistry, Geophysics, Geosystems 10.1029/2021GC010167

Hall, R. (2017). Southeast Asia: New views of the geology of the Malay Archipelago. Annual Review of Earth and Planetary Sciences, 45,
331–358. https://doi.org/10.1146/annurev-earth-063016-020633
Hall, R., & Spakman, W. (2015). Mantle structure and tectonic history of SE Asia. Tectonophysics, 658, 14–45. https://doi.org/10.1016/j.
tecto.2015.07.003
Harris, C. W., Miller, M. S., Supendi, P., & Widiyantoro, S. (2020). Subducted lithospheric boundary tomographically imaged beneath arc-conti-
nent collision in eastern Indonesia. Journal of Geophysical Research: Solid Earth, 125. https://doi.org/10.1029/2019JB018854
Harris, R. (2011). The nature of the Banda Arc–continent collision in the Timor region. In D. Brown & P. D. Ryan (Eds.), Arc continent collision.
Frontiers in Earth Sciences (pp. 163–211). Springer Berlin Heidelberg. https://doi.org/10.1007/978-3-540-88558-0_7
Hinschberger, F., Malod, J.-A., Réhault, J. P., Villeneuve, M., Royer, J. Y., & Burhanuddin, S. (2005). Late Cenozoic geodynamic evolution of
eastern Indonesia. Tectonophysics, 404, 91–118. https://doi.org/10.1016/j.tecto.2005.05.005
Hirth, G., & Kohlstedt, D. (2004). Rheology of the upper mantle and the mantle wedge: A view from the experimentalists. In J. Eiler (Ed.), Inside
the subduction factory (Vol. 138, pp. 83–106). American Geophysical Union. https://doi.org/10.1029/138GM06
Holmes, A. (1931). Radioactivity and Earth movements. Nature, 128, 419. https://doi.org/10.1038/128496e0
Husson, L. (2006). Dynamic topography above retreating subduction zones. Geology, 34, 741–744. https://doi.org/10.1130/g22436.1
Husson, L., Bernet, M., Guillot, S., Huyghe, P., Mugnier, J. L., Replumaz, A., et al. (2014). Dynamic ups and downs of the Himalaya. Geology,
42, 839–842. https://doi.org/10.1130/g36049.1
Husson, L., Boucher, F. C., Sarr, A. C., Sepulchre, P., & Cahyarini, S. Y. (2019). Evidence of Sundaland's subsidence requires revisiting its bio-
geography. Journal of Biogeography, 47, 843–853. https://doi.org/10.1111/jbi.13762
Jolivet, L., Faccenna, C., Becker, T., Tesauro, M., Sternai, P., & Bouilhol, P. (2018). Mantle flow and deforming continents: From India-Asia
convergence to Pacific subduction. Tectonics, 37, 2887–2914. https://doi.org/10.1029/2018tc005036
Kaus, B. J. P., Popov, A. A., Baumann, T. S., Pusok, A. E., Bauville, A., Fernandez, N., & Collignon, M. (2016). Forward and inverse modelling
of lithospheric deformation on geological timescales. Proceedings of NIC Symposium (Vol. 48, pp. 978–983).
Koulali, A., Susilo, S., McClusky, S., Meilano, I., Cummins, P., Tregoning, P., et al. (2016). Crustal strain partitioning and the associated earth-
quake hazard in the eastern Sunda-Banda Arc. Geophysical Research Letters, 43, 1943–1949. https://doi.org/10.1002/2016gl067941
Loiselet, C., Husson, L., & Braun, J. (2009). From longitudinal slab curvature to slab rheology. Geology, 37, 747–750. https://doi.org/10.1130/
g30052a.1
Mackwell, S. J., Zimmerman, M. E., & Kohlstedt, D. L. (1998). High-temperature deformation of dry diabase with application to tectonics on
Venus. Journal of Geophysical Research, 103(B1), 975–984. https://doi.org/10.1029/97JB02671
Mallard, C., Coltice, M., Seton, M., Müller, R. D., & Tackley, P. J. (2016). Subduction controls the distribution and fragmentation of Earth's
tectonic plates. Nature, 535, 140–143. https://doi.org/10.1038/nature17992
Miller, M. S., Zhang, P., Dahlquist, M. P., West, A. J., Becker, T. W., & Harris, C. W. (2021). Inherited lithospheric structures control arc-conti-
nent collisional heterogeneity. Geology, 49, 652–656.
Milsom, J. (2001). Subduction in eastern Indonesia: How many slabs? Tectonophysics, 338, 167–178. https://doi.org/10.1016/s0040-1951(01)
00137-8
Molnar, P., & Cronin, T. W. (2015). Growth of the maritime continent and its possible contribution to recurring ice ages. Paleoceanography, 30,
196–225. https://doi.org/10.1002/2014pa002752
Pedoja, K., Husson, L., Bezos, A., Pastier, A. M., Imran, A. M., Arias-Ruiz, C., et al. (2018). On the long-lasting sequences of coral reef ter-
races from SE Sulawesi (Indonesia): Distribution, formation, and global significance. Quaternary Science Reviews, 188, 37–57. https://doi.
org/10.1016/j.quascirev.2018.03.033
Porritt, R. W., Miller, M. S., O’Driscoll, L. J., Harris, C. W., Roosmawati, N., & Teofilo da Costa, L. (2016). Continent–arc collision in the Banda
Arc imaged by ambient noise tomography. Earth and Planetary Science Letters, 449, 246–258. https://doi.org/10.1016/j.epsl.2016.06.011
Pownall, J. M., Hall, R., & Lister, G. S. (2016). Rolling open Earth's deepest forearc basin. Geology, 44, 947–950. https://doi.org/10.1130/
g38051.1
Ranalli, G. (1995). Rheology of the earth (2nd ed.). Chapman & Hall.
Ricard, Y., Richards, M., Lithgow-Bertelloni, C., & Le Stunff, Y. (1993). A geodynamic model of mantle density heterogeneity. Journal of Geo-
physical Research, 98, 21895–21909. https://doi.org/10.1029/93jb02216
Royden, L. H., & Husson, L. (2009). Subduction with variations in slab buoyancy: Models and application to the Banda and Apennine sys-
tems. In S. Lallemand & F. Funiciello (Eds.), Subduction zone geodynamics (pp. 35–45). Springer Berlin, Heidelberg. https://doi.
org/10.1007/978-3-540-87974-9_2
Salles, T., Mallard, C., Husson, L., Zahirovic, S., Sarr, A.-C., & Sepulchre, P. (2021). Quaternary landscape dynamics boosted species dispersal
across Southeast Asia. Communications Earth & Environment, 2, 240. https://doi.org/10.1038/s43247-021-00311-7
Sandiford, M. (2007). The tilting continent: A new constraint on the dynamic topographic field from Australia. Earth and Planetary Science
Letters, 261, 152–163. https://doi.org/10.1016/j.epsl.2007.06.023
Sarr, A. C., Husson, L., Sepulchre, P., Pastier, A.-M., Pedoja, K., Elliot, M., et al. (2019). Subsiding Sundaland. Geology, 47, 119–122. https://
doi.org/10.1130/G45629.1
Sarr, A. C., Sepulchre, P., & Husson, L. (2019). Impact of the Sunda Shelf on the climate of the Maritime Continent. Journal of Geophysical
Research: Atmospheres, 124, 2574–2588. https://doi.org/10.1029/2018jd029971
Şengül-Uluocak, E., Pysklywec, R. N., Göğüş, O. H., & Ulugergerli, E. U. (2019). Multidimensional geodynamic modeling in the Southeast
Carpathians: Upper mantle flow-induced surface topography anomalies. Geochemistry, Geophysics, Geosystems, 20, 3134–3149.
Silver, E. A., Reed, D., McCaffrey, R., & Joyodiwiryo, Y. (1983). Back arc thrusting in the eastern Sunda arc, Indonesia: A consequence of
arc-continent collision. Journal of Geophysical Research, 88, 7429–7448. https://doi.org/10.1029/jb088ib09p07429
Simandjuntak, T. O., & Barber, A. J. (1996). Contrasting tectonic styles in the Neogene orogenic belts of Indonesia. Geological Society, London,
Special Publications, 106, 185–201. https://doi.org/10.1144/gsl.sp.1996.106.01.12
Smyth, H., Hall, R., Hamilton, J., & Kinny, P. (2008). East Java: Cenozoic basins, volcanoes and ancient basement. In Proceedings Indonesia
Petroleum Association, 30th Annual Convention (pp. 251–266).
Solihuddin, T., Collins, L. B., Blakeway, D., & O’Leary, M. J. (2015). Holocene coral reef growth and sea level in a macrotidal, high turbidity set-
ting: Cockatoo Island, Kimberley Bioregion, northwest Australia. Marine Geology, 359, 50–60. https://doi.org/10.1016/j.margeo.2014.11.011
Spakman, W., & Hall, R. (2010). Surface deformation and slab–mantle interaction during Banda arc subduction rollback. Nature Geoscience, 3,
562–566. https://doi.org/10.1038/ngeo917
Susilohadi, S., & Soeprapto, T. A. (2015). Plio-Pleistocene seismic stratigraphy of the Java Sea between Bawean Island and East Java. Marine
Geology of Indonesia, 32, 5–16.
Turcotte, D., & Schubert, G. (1982). Geodynamics. John Wiley and Sons.

HUSSON ET AL. 16 of 17
 15252027, 2022, 1, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2021GC010167 by Nat Prov Indonesia, Wiley Online Library on [11/04/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
Geochemistry, Geophysics, Geosystems 10.1029/2021GC010167

Waelbroeck, C., Labeyrie, L., Michel, E., Duplessy, J. C., McManus, J. F., Lambeck, K., et  al. (2002). Sea-level and deep water tempera-
ture changes derived from benthic foraminifera isotopic records. Quaternary Science Reviews, 21, 295–305. https://doi.org/10.1016/
s0277-3791(01)00101-9
Wessel, P., Luis, J. F., Uieda, L., Scharroo, R., Wobbe, F., Smith, W. H. F., & Tian, D. (2019). The generic mapping tools version 6. Geochemistry,
Geophysics, Geosystems, 20, 5556–5564. https://doi.org/10.1029/2019gc008515
Widiyantoro, S., Pesicek, J. D., & Thurber, C. H. (2011). Subducting slab structure below the eastern Sunda arc inferred from non-linear seismic
tomographic imaging. Geological Society, London, Special Publications, 355, 139–155. https://doi.org/10.1144/sp355.7
Widiyantoro, S., & van der Hilst, R. (1997). Mantle structure beneath Indonesia inferred from high-resolution tomographic imaging. Geophysical
Journal International, 130, 167–182. https://doi.org/10.1111/j.1365-246x.1997.tb00996.x
Yang, T., Gurnis, M., & Zahirovic, S. (2016). Mantle-induced subsidence and compression in SE Asia since the early Miocene. Geophysical
Research Letters, 43, 1901–1909. https://doi.org/10.1002/2016gl068050
Yang, X., Singh, S. C., & Tripathi, A. (2020). Did the Flores backarc thrust rupture offshore during the 2018 Lombok earthquake sequence in
Indonesia? Geophysical Journal International, 221, 758–768. https://doi.org/10.1093/gji/ggaa018
Zenonos, A., De Siena, L., Rawlinson, N., & Rawlinson, N. (2019). P and S wave travel time tomography of the SE Asia-Australia collision zone.
Physics of the Earth and Planetary Interiors, 293, 106267. https://doi.org/10.1016/j.pepi.2019.05.010

Reference From the Supporting Information

Kreemer, C., Blewitt, G., & Klein, E. C. (2014). A geodetic plate motion and global strain rate model. Geochemistry, Geophysics, Geosystems,
15, 3849–3889. https://doi.org/10.1002/2014GC005407

HUSSON ET AL. 17 of 17