During the latest Cretaceous–early Cenozoic, the northern margin of the Australian plate was characterized by a large (4000 km wide) north- to northeast-dipping subduction zone (New Guinea–Pocklington subduction zone) consuming a marginal basin. Geological and geophysical data imply that the subduction zone was active ca. 71–50 Ma, and suggest that it was responsible for plate acceleration from ∼1.0 to ∼7.3 cm/yr ca. 64–59 Ma, and plate deceleration from ∼7.3 to ∼0.3 cm/yr at 52–49 Ma. This paper presents a numerical model of buoyancy-driven subduction to test if the rates of Australian plate acceleration and deceleration can be ascribed to the progressive evolution of a subducting slab. The geodynamic model reproduces the first-order plate velocity evolution of the Australian plate, with a transient ∼5 m.y. period of acceleration from 2 to 8 cm/yr during upper mantle slab lengthening, an ∼5 m.y. period of rapid plate motion (∼5–8 cm/yr), and a short, 3.9 m.y., period of plate deceleration, starting with a 2 cm/yr velocity drop during 3.1 m.y. of continental subduction and followed by ∼0.8 m.y. of rapid deceleration (4 cm/yr velocity drop) during slab detachment. The geodynamic model demonstrates that plate velocity increases or decreases of ∼4–6 cm/yr can occur over a period lasting <1 m.y. to a few million years, comparable to what is observed for the latest Cretaceous–early Cenozoic evolution of the Australian plate. Such rates of plate acceleration and deceleration could be tested against plate kinematic data for other subduction settings on Earth.


The Earth’s lithosphere is segmented into numerous larger and smaller plates (Bird, 2003) that move at a wide variety of plate tectonic speeds, to as much as ∼10 cm/yr. The velocities of these plates have traditionally been calculated in hotspot reference frames, such as the Pacific (Minster and Jordan, 1978; Gripp and Gordon, 2002; Wessel and Kroenke, 2008), Indo-Atlantic (e.g., O’Neill et al., 2005), or global (e.g., Gordon and Jurdy, 1986; Doubrovine et al., 2012) hotspot reference frames, as well as no-net-rotation reference frames (e.g., Argus and Gordon, 1991; Kreemer et al., 2003; Argus et al., 2011). More recently, other geological and geophysical features have been used to constrain the absolute velocities of the tectonics plates, such as subducted slabs in the upper and lower mantle (e.g., van der Meer et al., 2010; Schellart, 2011; Butterworth et al., 2014), upper mantle seismic anisotropy (e.g., Long and Silver, 2009; Kreemer, 2009; Conrad and Behn, 2010), and spreading ridges at the Earth’s surface (e.g., Becker et al., 2015; Wessel and Müller, 2016). A number of conceptual and geodynamic models have been put forward to explain the variations in plate velocities. Although such models vary as to what might be the main physical mechanism that explains variation in absolute plate velocities, such as subducting plate age (e.g., Carlson et al., 1983; Goes et al., 2008), slab width (Schellart et al., 2010), plate circumference attached to a slab (Forsyth and Uyeda, 1975; Gripp and Gordon, 2002), double subduction (Jagoutz et al., 2015), and the transition from normal (i.e., oceanic) subduction to collision and/or continental subduction (e.g., Patriat and Achache, 1984; Klootwijk et al., 1992; Molnar and Stock, 2009; van Hinsbergen et al., 2011; White and Lister, 2012; Zahirovic et al., 2012; Capitanio and Replumaz, 2013), the models generally have the common element of subduction to explain rapid plate motion (e.g., Conrad and Lithgow-Bertelloni, 2002). An important exception is the so-called plume push force, which has been proposed to explain the extreme plate velocity of the Indian plate in the latest Cretaceous–earliest Cenozoic (Cande and Stegman, 2011; van Hinsbergen et al., 2011). The ridge push force is generally thought to be an order of magnitude smaller than the slab pull force (Forsyth and Uyeda, 1975).

The Australian plate is home to the Australian continent, which is currently the fastest moving continent on Earth (Keep and Schellart, 2012) with north- to northeast-directed velocities as high as 6.5 cm/yr (Fig. 1). These high velocities can be ascribed to the wide subduction zones that are present along the northern boundary of the Australian plate, namely the north- to northeast-dipping Sunda subduction zone in the northwest and the north- to northeast-dipping Melanesian subduction zone in the northeast. Over the past ∼75 m.y., the Australian plate has generally undergone rapid plate velocities, except for 2 periods, at 73–63 Ma and 49–40 Ma, when plate velocities reached minima of 0.8–2.9 cm/yr and 0.3–1.6 cm/yr, respectively, in the Indo-Atlantic moving hotspot reference frame (Fig. 2). During the intervening period the velocities were much higher with a maximum of 7.3 cm/yr at 59–52 Ma. During the latest Cretaceous (Maastrichtian) to early Eocene, the Australian plate was bordered to the east, south, and west by spreading ridges in the Tasman Sea, Southern Ocean, and Indian Ocean (e.g., Royer and Sandwell, 1989; Yan and Kroenke, 1993; Gaina et al., 1998; Hall, 2012; Seton et al., 2012); this cannot explain the period of rapid plate motion and the preceding and following phases of plate acceleration and deceleration. In Schellart and Spakman (2015), a 4000-km-wide (trench-parallel extent) subduction zone was identified in the New Guinea region active ca. 71–50 Ma, and the phases of plate acceleration, rapid plate motion, and plate deceleration of the Australian plate were ascribed to the geodynamic evolution of this subduction zone. The subduction zone consumed a small oceanic backarc basin, the Emo backarc basin, with a maximum north-south extent of ∼1100 km. The phase of plate acceleration was ascribed (in Schellart and Spakman, 2015) to progressive upper mantle slab lengthening during the initial transient subduction phase and the phase of plate deceleration was ascribed to slab detachment during the final stage of subduction.

I present here a geodynamic numerical model of buoyancy-driven subduction to test if the rates of plate acceleration and deceleration, and the high rate of Australian plate motion during the intervening period, can be ascribed to the progressive evolution of the subduction zone in the New Guinea region, from transient slab lengthening to the final stage of slab detachment. The geodynamic model provides new insight into the rates of plate acceleration and deceleration during subduction of the edge of a plate that contains a large continent.


Numerical models are presented that have specifically been designed to investigate the subduction evolution of a small oceanic backarc basin (Emo backarc basin), with a particular focus on the latest stage of subduction when all oceanic lithosphere is consumed and continental passive margin lithosphere approaches the trench. The models use the code Underworld (Stegman et al., 2006; Moresi et al., 2007, 2014), in which plate motion, subduction and mantle flow are modeled in a Cartesian domain using compositional buoyancy contrasts in an incompressible Boussinesq fluid at a very low Reynolds number. Distinct volumes (e.g., continental upper crust, lower crust and lithospheric mantle, sublithospheric upper mantle, lower mantle) are represented by sets of Lagrangian particles that are embedded within a standard Eulerian finite element mesh, which discretizes the problem to solve the governing equations. For additional information on the numerical technique and the nondimensional equations, see Moresi et al. (2007, 2014) and Stegman et al. (2006). Velocities in the models are scaled following the scaling formulations presented in Schellart and Moresi (2013).

The modeling domain is 1100 km deep and 6200 km long and has free-slip boundary conditions along the top surface, bottom surface, and lateral side walls (Fig. 3). Mesh resolution in the 6200 × 1100 km numerical domain is 1024 (length) by 256 (depth) elements. An adaptive spatial mesh refinement has been implemented to increase resolution in the region of the subduction zone and close to the surface. The smallest cell dimensions in this central region (rectangle with dotted outline in Fig. 3), 220 km thick and 1860 km wide, are 2.79 km (length) by 2.15 km (depth). Initial particle distribution is 20 particles per cell (total of 5,242,880 particles).

The reference model A involves a layered mantle volume incorporating a free, four-layer subducting plate, a fixed, neutrally buoyant, non-stratified overriding plate, a low-viscosity sublithospheric upper mantle, and a high-viscosity lower mantle. The top 100 km of the model domain consists of a low-viscosity, low-density air layer to provide a free surface to the tectonic plates (sticky air layer following Schmeling et al., 2008; Crameri et al., 2012). The subduction zone interface is implemented in the same way as in the numerical models in Schellart and Moresi (2013) and the laboratory models of Duarte et al. (2013, 2015), in which a weak zone develops that consists of the weak materials from the top layer of the subducting plate, which allows for progressive one-sided subduction. The subducting plate has a negatively buoyant, 1120-km-long oceanic lithosphere segment representing the Emo oceanic backarc basin and a trailing, positively buoyant 3600-km-long continental lithosphere segment representing the Australian continent (which includes a 100 km passive margin). The length of the oceanic lithosphere segment is based on the reconstruction from Schellart and Spakman (2015), while the length of the continental part is based on the north-south extent of the Australian continent in the latest Cretaceous–earliest Cenozoic.

The oceanic part of the subducting plate is 80 km thick and consists of 4 layers with their own thickness and rheology following earlier modeling approaches (e.g., Schellart et al., 2010; Stegman et al., 2010), except that the models presented here use viscoplastic rheologies with a von Mises yield stress throughout the lithosphere to facilitate possible detachment. The physical properties of the different layers are listed in Table 1. The 80 km thickness of the oceanic lithosphere is based on the age of the Emo backarc basin, which is derived from the age of the protoliths of the Emo metamorphics and Owen Stanly metamorphic rocks that are thought to originate from this basin. Early research on the metamorphics proposed a Mesozoic age (Davies and Warren, 1992) or Late Cretaceous age (Worthing and Crawford, 1996). In more recent work, the plate tectonic reconstruction from Hall (2012) implies that an ocean basin or backarc basin north of the New Guinea passive margin already existed in the Late Jurassic–Early Cretaceous. In other recent work, the age of the protolith of the metamorphics has been described as Cretaceous (Baldwin et al., 2012), while Davies (2012) recorded a mid-Cretaceous age with radiometric dates of 120–107 Ma. This would imply an oceanic backarc basin with an age range of ∼36–70 million years at the time of subduction ca. 71–50 Ma. Considering this age range, an 80-km-thick oceanic lithosphere (implying an average age of ∼50 million years at the time of subduction) is justified. An average density contrast of 80 kg/m3 between the oceanic slab and ambient sublithospheric mantle is adopted in the models, as it is assumed that there is pervasive eclogitization of the basaltic crust into eclogite facies rocks (Cloos, 1993).

The continental part of the subducting plate, which includes the 100-km-long passive margin, has the same rheological layering with the same thickness as the oceanic part (Fig. 3, bottom panel). The only difference between the oceanic and continental domains is the density structure, with a lighter, positively buoyant, continental lithosphere due to a 30 km continental crust, a denser, negatively buoyant, oceanic lithosphere (Table 1), and a passive margin across which the density changes linearly due to a linear change in continental crustal thickness from 30 km at the continental side to 10 km at the oceanic side. This set-up was specifically chosen to investigate where strain localization and slab detachment would occur in a scenario without any lateral rheological contrasts, but with only lateral density contrasts in the subducting plate.

The model overriding plate is neutrally buoyant, implying that it is oceanic in nature with a density that is comparable to that of the sublithospheric mantle. This is most likely if one considers the tectonic setting to the north of Australia in the latest Cretaceous–early Eocene, where reconstructions indicate that there was an ocean to the north (Hall, 2012). Another reason for implementing a neutral buoyancy was to minimize the effect of the overriding plate on the subduction process, such that the change from oceanic to continental subduction and the process of slab detachment could be studied without other complicating factors. For the rheology the simplest set-up was chosen (following Duarte et al., 2013; Schellart and Moresi, 2013) that allows for continuous and progressive subduction, because little is known about the properties and rheology of the overriding plate. The model overriding plate has a forearc region with a higher viscosity than the backarc region (Table 1) (following Schellart and Moresi, 2013), as it is generally thought that the forearc is relatively cool and strong with depressed isotherms due to the cold subducting plate diving below it, and therefore has a higher effective viscosity than the arc and backarc regions, where isotherms are elevated.

Apart from reference model A, one other model is presented, model B, with a higher yield stress compared to the reference model (Table 1). A three-dimensional geometrical set-up would be ideal for the geodynamic models, as presented in other modeling work of subduction and slab detachment (e.g., van Hunen and Allen, 2011; Capitanio and Replumaz, 2013; Chertova et al., 2014; Capitanio et al., 2015), because such models can take into account three-dimensional (3D) spatial effects, including lateral migration of slab detachment, tearing, and trench-parallel flow of mantle material. However, this was found to be unfeasible at this time, considering the very large size of the subduction zone (∼4000 km), the required resolution, and the required number of time steps (several thousand). These earlier works generally used lower spatial resolution and modeled smaller plates and subduction zones. Furthermore, considering that the natural prototype is a wide subduction zone, its geodynamic behavior, in particular that of its central portion, can be approximated with a model using a 2D spatial set-up (Schellart and Moresi, 2013).


Reference model A shows a general style of subduction with progressive slab lengthening and an increase of the slab dip angle (averaged along the entire slab length) during the free sinking phase, from 21° at the start (Fig. 3) to 40° at an intermediate stage (Fig. 4A), to 46° at an advanced stage (Fig. 4B). As the slab approaches the 660 discontinuity its tip is slightly deflected to a lower slab dip angle. During the free sinking phase the trench-normal subducting plate velocity vSP⊥ increases rapidly to a maximum of ∼8 cm/yr, when the slab tip is located at ∼578 km depth (82 km above the 660 km discontinuity) (Fig. 5A). The trench-normal trench velocity (vT⊥) is significantly less than vSP⊥ and increases to a maximum of ∼2.5 cm/yr, which is reached at a time before the maximum of vSP⊥. During the slab tip–660 km discontinuity interaction phase (∼8–12 m.y.), vSP⊥ decreases from ∼8 to 5–6 cm/yr, while vT⊥ remains relatively unchanged. From ∼12.2 to 15.3 m.y. continental subduction takes place (Figs. 6A, 6B), during which vSP⊥ decreases from ∼5 to 3 cm/yr. This is followed by a short phase of slab detachment from 15.3 to 16.1 m.y. (Figs. 6C–6E), during which vSP⊥ decreases from ∼3 to –1 cm/yr and then increases again to ∼0 cm/yr (Fig. 5A). During continental subduction and slab detachment, vT⊥ decreases from ∼2 cm/yr to 0 cm/yr. After slab detachment, subduction stops, vSP⊥ and vT⊥ are close to zero, while the detached slab continues to sink into the mantle (Figs. 4H, 4I, and 6F).

During the phase of continental subduction, the continental lithosphere is pulled into the mantle by the oceanic slab segment, and the maximum subduction depth of the continental crust (its base, originally at 30 km depth), just before slab detachment is complete, is 96 km. This indicates a depth increase of 66 km. The slab detachment process starts with the formation of two conjugate shear zones (Fig. 6C) and occurs mostly by symmetrical necking (Figs. 6C–6E). Once slab detachment is complete, the trailing slab segment slightly rebounds.

The other model (model B with a stronger subducting plate than model A) generally shows the same evolution as model A until ∼12 m.y.; model B starts continental subduction just after ∼12 m.y., and this continues until ∼20 m.y. when vSP⊥ has decreased to zero, vT⊥ is slightly negative (∼–0.5 cm/yr), continental subduction stops and the only remaining activity is the very slow progressive steepening of the slab (Figs. 5B and 7). The maximum subduction depth of the continental crust (its base, originally at 30 km depth) at the end of the model run is 144 km.


Plate Velocity Changes in Model and Nature

The subduction models A and B show a general style of subduction behavior during the oceanic subduction period that is comparable to previous buoyancy-driven subduction modeling studies, with progressive slab steepening (e.g., Jacoby, 1973; Funiciello et al., 2004; Schellart, 2004a; Guillaume et al., 2009) and an increasing trenchward subducting plate velocity during the free sinking phase to a maximum velocity when the slab tip is just above (<∼100 km from) the 660 km viscosity discontinuity (e.g., Schellart, 2008; Capitanio et al., 2010; Schellart et al., 2011; Meyer and Schellart, 2013; Chen et al., 2015, 2016) (Fig. 5). In the current models it takes only ∼5 m.y. to increase from vSP⊥ = ∼2 cm/yr to ∼8 cm/yr. Such a velocity increase over such a time period is comparable to that implied by the velocity curve for the Australian plate in the Indo-Atlantic moving hotspot reference frame (Fig. 2), showing an increase from ∼1 to ∼7.4 cm/yr between ca. 64 and ca. 59 Ma and an increase in the northward velocity component, which is approximately perpendicular to the strike of the fossil subduction zone, from ∼0.2 to ∼5.0 cm/yr. The increase in subducting plate velocity in the models and for the Australian plate can be ascribed to the progressive lengthening of the slab and deepening of the slab tip (Figs. 4A, 4B, and 5), and thereby the increase in slab negative buoyancy and net slab pull (Schellart, 2004b). These high subducting plate velocities can be reached, despite the fact that the plate is relatively long in the latest Cretaceous–early Cenozoic (∼4000–5000 km long Australian plate in the north-south direction).

The decrease in subducting plate velocity during the slab tip–660 km discontinuity interaction phase is generally also observed in earlier works (e.g., Schellart, 2008; Capitanio et al., 2010; van Hunen and Allen, 2011; Schellart et al., 2011; Meyer and Schellart, 2013; Chen et al., 2015, 2016). Reference model A further shows a period of decrease in vSP⊥ from ∼4.5 cm/yr to ∼3 cm/yr during ∼3.1 m.y. of continental subduction, followed by an ∼0.8 m.y. phase of slab detachment, during which vSP⊥ decreases from ∼3 cm/yr to –1 cm/yr (Fig. 5A). The slowdown during continental subduction can be explained by the decrease in buoyancy force of the slab (i.e., becoming less negative) due to addition of positively buoyant continental slab material (Figs. 6A–6C, images on left). The rapid slowdown during slab detachment can be explained by the rapid decrease in coupling between the detaching slab segment and the trailing slab segment due to the rapid decrease in effective viscosity in the detachment zone (Figs. 6B–6E, images on right). During the ∼3.9 m.y. period of continental subduction and slab detachment the velocity decreases ∼5.5 cm/yr, which is comparable to the ∼3 m.y. period from ca. 52 Ma to ca. 49 Ma during which the Australian plate velocity decreased ∼7 cm/yr (Fig. 2A) and the northward velocity component of the Australian plate decreased ∼5 cm/yr (Fig. 2B) in the Indo-Atlantic moving hotspot reference frame. Subduction model A thereby demonstrates the physical viability of the conceptual model in which slowdown of the Australian plate in the latter part of the early Eocene is ascribed to continental subduction and subsequent slab detachment, as illustrated in Figures 5A and 6.

In the preceding two paragraphs the subducting plate velocities from the numerical models have been compared to the subducting Australian plate velocities in the Indo-Atlantic moving hotspot reference frame from O’Neill et al. (2005). As explained in detail elsewhere (e.g., Schellart et al., 2008; Schellart, 2011; Schellart and Spakman, 2012, 2015), this Indo-Atlantic hotspot reference frame is preferred over other reference frames for calculating plate velocities and plate boundary velocities, which is particularly evident for the western Pacific domain. The reasons are many, including minimization of global trench migration velocities and rollback-induced toroidal volume fluxes in this reference frame (Schellart et al., 2008), optimal agreement between predicted and observed fossil slab locations in the southwest Pacific (Schellart and Spakman, 2012) and slab structures in the western Pacific (Schellart, 2011), and highest coincidence with fossil slab sinking-induced dynamic topography in the southwest Pacific (Schellart and Spakman, 2015). For comparison, the velocities have also been plotted in the global moving hotspot reference frame from Doubrovine et al. (2012). The velocity profile for the northward velocity component is roughly comparable in shape to the Indo-Atlantic velocity profile in the period 65–49 Ma (Fig. 2B), but the total acceleration, peak velocity, and total deceleration are reduced.

In Schellart and Spakman (2015) geological and geophysical evidence was presented to propose the existence and the timing of activity (ca. 71–50 Ma) of the fossil New Guinea–Pocklington subduction zone; it was argued that this subduction zone was largely responsible for velocity changes of the Australian plate, with an increase in speed from ∼64 Ma to ∼59 Ma due to slab lengthening, rapid plate motion ca. 59–52 Ma during mature subduction, and slowing ca. 52–49 Ma due to slab detachment. The reference numerical model A is largely consistent with this conceptual model, but also shows that part of the slowing (∼27%) can be ascribed to the phase of continental subduction prior to slab detachment.

The numerical geodynamic model can explain the subducting plate velocity changes and velocity magnitudes in the period 65–45 Ma, but does not explain the significant north-directed plate acceleration that started ca. 44 Ma. This new phase of plate acceleration has been ascribed to the formation of new subduction zones along the northwestern, northern, and northeastern boundaries of the Australian plate between ca. 50 and ca. 40 Ma (Schellart and Spakman, 2015), and this phase of renewed subduction has not been included in the numerical model. Formation of the north- to northeast-dipping Sunda subduction zone ca. 45 Ma (Hall, 2012), North Sulawesi–Halmahera subduction zone ca. 45 Ma (Hall, 2012), and northeast- to east-dipping New Caledonia–Northland subduction zone at 50–40 Ma (Schellart and Spakman, 2012) would have resulted in increasing slab pull forces, much like the free slab sinking phase and subducting plate acceleration shown in Figure 5, causing fast plate acceleration at 45–40 Ma and fast north- to northeast-directed plate velocities since ca. 40 Ma.

Continental Subduction and Slab Detachment

Velocities and Duration of Continental Subduction and Slab Detachment

Reference model A with slab detachment shows a phase of continental subduction and slab detachment lasting only ∼3.9 m.y., with a decrease in vSP⊥ of ∼5.5 cm/yr. Earlier works on slab detachment generally show comparable reductions in subducting plate velocity (e.g., Burkett and Billen, 2010; van Hunen and Allen, 2011; Capitanio and Replumaz, 2013; Capitanio, 2014). Normal (oceanic) subduction models from Burkett and Billen (2010) show comparable reductions in vSP⊥ of ∼5–9 cm/yr during slab detachment, as do the continental subduction break-off models of Capitanio (2014) and Capitanio and Replumaz (2013) with a reduction in vSP⊥ of ∼5 cm/yr and a reduction in convergence velocity of 2–4 cm/yr, respectively. The continental subduction break-off model of van Hunen and Allen (2011) shows a higher reduction of ∼14 cm/yr, but their model generally shows much faster subducting plate velocities, to ∼23 cm/yr.

The duration of slab detachment for model A is relatively short (∼0.8–1.0 m.y.), but similar short durations (∼0.3–3.0 m.y.) have been reported in earlier works on 2D numerical modeling of subduction and slab detachment (Duretz et al., 2012, 2014), as well as 2D and 3D numerical modeling of normal (oceanic) subduction (∼1–3 m.y.) by Burkett and Billen (2010). The 2D numerical models of subduction followed by continental subduction of van Hunen and Allen (2011) report a long delay time between the start of continental subduction and slab detachment (in the range ∼9–23 m.y.), much longer than reported here (∼3.9 m.y.), which is possibly due to the higher yield stress in their models (400 MPa) compared to the one used here for reference model A (352.8 MPa). Model B, with a yield stress of 382.2 MPa, did not show any yielding for the duration of the simulation, which was allowed to continue for ∼9.5 m.y. since the start of continental subduction.

The models presented here have a 2D spatial set-up, and thus represent a subduction scenario in which slab detachment occurs contemporaneously along the entire lateral extent of the subduction zone. Earlier work, however, proposed that slab detachment and slab tears migrate subhorizontally (e.g., Wortel and Spakman, 2000). Numerical subduction models with a 3D spatial set-up have reproduced this process, demonstrating that lateral tear migration can be very rapid, as much as ∼80 cm/yr (van Hunen and Allen, 2011). For a 4000-km-wide subduction zone this would imply that slab detachment can take place in 2.5 m.y. if it initiated at one location in the center and migrated outward. It is plausible that along such a wide subduction zone, slab detachment starts at multiple locations, and migrates laterally from these places, thereby reducing the time of detachment along the entire extent of the subduction zone.

Earlier work has found that temperature effects on slab detachment duration are relatively small (decrease of ∼8%–10%) (Gerya et al., 2004; Duretz et al., 2012), indicating that the duration of slab detachment for the isothermal model A presented here is of the right order of magnitude and would only be marginally affected in case thermal gradients were incorporated. Instead, the variability in duration of slab detachment has been ascribed mainly to differences in rheology and yield stress of the slab material (e.g., Andrews and Billen, 2009; Duretz et al., 2012). The yield stress of 352.8 MPa used here for model A with slab detachment is comparable to that reported by Andrews and Billen (2009; 300 MPa) to allow for slab detachment.

Model B without slab detachment shows a decrease in vSP⊥ from ∼4.5 cm/yr to ∼0 cm/yr during continental subduction over an ∼8 m.y. period. Laboratory models of subduction and late-stage continental subduction from Edwards et al. (2015) also lack slab detachment, and show periods of continental subduction that generally last longer, ∼10–30 m.y. This can be partly ascribed to the neutral buoyancy of the continental lithospheric mantle in the current models compared to the negative buoyancy for the models from Edwards et al. (2015), giving the continental lithosphere in the current work a larger net positive buoyancy force, thereby reducing continental subduction. Most of the discrepancy, however, can be ascribed to the fact that 14 out of 15 models in Edwards et al. (2015) have a larger passive margin (250–500 km). The one model with the same passive margin size as here (100 km) has a very comparable continental subduction period of ∼10 m.y.

Forces During Continental Subduction and Slab Detachment

In reference model A, the slab detaches in a late stage of continental subduction, while there is no detachment in model B. This can be explained by the buoyancy forces that operate during continental subduction and the different yield stress in models A and B. Figure 8A shows the evolution of the buoyancy force of the oceanic slab segment FBu-Oc (negative) (per meter length of the trench) and the buoyancy force of the continental slab segment FBu-Co (positive) (per meter length of the trench). These forces have been calculated with the following equations: 
where ΔρOc is the density contrast between ambient mantle and oceanic slab (ΔρOc = –80 kg/m3), LOc is the length of the oceanic slab segment that contributes to the buoyancy force (so excluding the segment subhorizontally overlying the 660 km discontinuity; this parameter changes with time), TOc is the thickness of the oceanic slab segment (changes with time and space), g is gravitational acceleration (9.8 m/s2), ΔρCo is the density contrast between ambient mantle and continental crust (ΔρCo = 400 kg/m3), LCo is the length of the continental slab segment (changes with time), and TCo is the thickness of the continental crustal slab segment (changes with time and space; note that the continental lithospheric mantle is neutrally buoyant). As can be observed, subduction of oceanic lithosphere creates a negative (downward) buoyancy force with a minimum of ∼–5.2 × 1013 N/m that is reached between the start of continental subduction and the start of slab detachment. Subduction of continental lithosphere creates a positive (upward) buoyancy force with a maximum of ∼6.4 × 1012 N/m that is reached during slab detachment. The differential tensional buoyancy force (per meter trench length) between the continental slab segment and oceanic slab segment, with FBu-Diff = FBu-CoFBu-Oc, is focused in the boundary region between the oceanic and continental slab segments and reaches a peak at the onset of slab detachment with FBu-Diff = ∼5.6 × 1013 N/m (Fig. 8A).

The strength of the slab is effectively determined by the yield stress of the strong core layer of the subducting plate (layer 3 and Top LithM in Table 1), which is 352.8 MPa. With such a yield stress and a 15–20-km-thick core layer in the vicinity of the ocean-continent transition during the continental subduction phase, this gives a yield force, the maximum in-plane slab pull force (per meter trench length) in the slab, FS-Max, of 5.3–7.1 × 1012 N/m. The evolution of this yield force for model A is plotted in Figure 8B. We can see that FS-Max is about an order of magnitude smaller than FBu-Diff, which might suggest at first that the slab is much too weak to sustain such a large buoyancy force and that slab detachment should have occurred in a much earlier phase of subduction. A significant part of the buoyancy force, however, is sustained elsewhere in the system, in particular by viscous forces in the ambient mantle (Conrad and Lithgow-Bertelloni, 2002; Schellart, 2004b; Krien and Fleitout, 2008; Leng and Zhong, 2010). Earlier work showed that only ∼10%–12% of the buoyancy forces is transmitted in the plane of the subducted slab toward the trailing plate (Schellart, 2004b; Sandiford et al., 2005). This would imply that the effective slab pull force (per meter trench length) FSP-Eff = 0.10–0.12 × FBu-Diff. Figure 8B shows the evolution of FSP-Eff, demonstrating that a maximum effective slab pull force is reached at the start of slab detachment, with FSP-Eff = 5.6–6.7 × 1012 N/m. Figures 6C and 8 further show that slab detachment starts at a time when a significant section of continental lithosphere has been subducted (LCo ≈128 km) and when FS-Max ≈0.10FBu-Diff, indicating that at this time the effective slab pull force is large enough to start rapid shearing, necking, and yielding in the slab. Before slab detachment, the amount of subducted continental lithosphere was insufficient to create a FBu-Diff of sufficient magnitude, and the slab was not sufficiently stretched yet to decrease FS-Max, giving FSP-Eff < FS-Max, and so slab detachment was not possible yet.

Applying calculations as above to model B, for a late stage as depicted in Figure 7C, FBu-Oc = –3.40 × 1013 N/m (±9.5%), FBu-Co = 1.07 × 1013 N/m (±17.2%), and FBu-Diff = 4.47 × 1013 N/m (±26.6%). This would imply an effective slab pull force (per meter trench length) FSP-Eff = 4.47–5.37 × 1012 N/m. This is less than the estimated yield force (per meter trench length) in the slab of FS-Max = 7.64 × 1012 N/m (±5%), thereby explaining why slab detachment did not occur in model B.

Continental Subduction Depth

Models A and B show different maximum depths of continental subduction, 96 km and 144 km, respectively, which can simply be explained by the weaker subducting plate in A compared to B. The weaker plate in A resulted in slab detachment prior to the continental slab segment reaching its maximum depth that could have been reached with the available negative buoyancy force of the oceanic slab segment. This resulted in a shorter duration of slab pull during continental subduction and thus a shallower depth of continental subduction.

The observed continental subduction depths reported here are comparable to, although in the lower range of, depths reported in previous modeling studies of buoyancy-driven subduction (e.g., van Hunen and Allen, 2011; Sizova et al., 2012; Edwards et al., 2015). For example, Edwards et al. (2015) reported depths in the range 70–560 km for experiments with different passive margin width, continental crustal thickness, and oceanic lithospheric thickness. Their experiment with a comparable passive margin width (100 km), continental crustal thickness (30 km), density contrast (366 kg/m3), and oceanic lithospheric thickness (100 km) showed a low continental subduction depth (70 km), comparable to the work reported here.


In this work a numerical geodynamic model of subduction, continental subduction, and slab detachment has been presented to simulate the subduction evolution at the northern margin of the Australian plate during the latest Cretaceous and early Cenozoic, and to quantify the velocity of the Australian plate. The numerical model produces a comparable subducting plate velocity evolution as documented for the Australian plate between ca. 66 and 49 Ma, with a first phase of plate acceleration due to progressive upper mantle slab lengthening, which increases the slab pull force, a period of relatively rapid plate motion during mature subduction, and a phase of plate deceleration due to continental subduction and subsequent slab detachment. Most (∼73%) of the plate deceleration is caused during the short (∼0.8 m.y.) phase of slab detachment, while a smaller component of plate deceleration (∼27%) is caused by the preceding phase of continental subduction, which lasted ∼3.1 m.y. A geodynamic subduction model that shows continental subduction but lacks a final phase of slab detachment is characterized by a slower and more gradual plate deceleration that lasts ∼8 m.y., which is in conflict with observations of fast plate deceleration for the Australian plate between ca. 52 and 49 Ma. It is thus concluded that slab detachment is an essential element to explain the rapid plate deceleration of the Australian plate. This conclusion is consistent with the evidence presented in Schellart and Spakman (2015) that the fossil New Guinea–Pocklington slab detached in the early Eocene and is currently located in the upper part of the lower mantle below central and southeastern Australia.

More generally, I propose here that rapid plate deceleration (∼3–5 cm/yr) over periods lasting not more than ∼1 m.y. can be ascribed to the process of slab detachment following continental subduction. It is thus conceivable that rapid decelerations that have occurred in the geological past for other tectonic plates have a comparable underlying cause, namely slab detachment.


I thank Anne Replumaz and two anonymous reviewers for their constructive and helpful comments. This research was supported through a Vici Fellowship from the Dutch Science Foundation (NWO) and has been supported by computational resources from the NCI National Facility in Australia through the National Computational Merit Allocation Scheme (projects ei8 and qk0).