Mantle dynamics of continentwide Cenozoic subsidence and tilting of Australia

Dynamic topography is an important control on Earth's overall topography and the gravity field (Cazenave et al., 1989; Christensen, 1998; Colin and Fleitout, 1990; Hager et al., 1985; LeStunff and Ricard, 1995). An intriguing implication of this concept is the contribution of mantle processes toward the formation of intracontinental basins, which may be recorded by continental flooding (Artemieva, 2007; Gurnis, 1990; Mitrovica et al., 1989; Pysklywec and Mitrovica, 1998; Spasojević et al., 2008). Although mantle-driven topography on continents is a global phenomenon (Heine et al., 2008), it is generally difficult to quantitatively link continental geology to mantle processes and to resolve the dynamic component of amplitude and rate of change. This is largely due to our limited ability to separate a dynamic component of topography from other crustal and lithospheric processes. In addition, geodynamic models are limited by their input conditions and their ability to produce realistic magnitudes of dynamic topography to match surface observations (Billen et al., 2003).

While it may be difficult to isolate the mantle-driven topography from changes in lithospheric and crustal thickness (Wheeler and White, 2000), Australia is generally free from such difficulties, since it is a stable continent with distal plate boundaries (Fig.1). However, since the Eocene, rapid seafloor spreading on the Southeast Indian Ridge has moved the Australian continent away from Antarctica and toward the subduction zones in SE Asia. Indeed, despite the absence of proximal plate boundaries, it has been suggested that Australia has undergone both a broad-scale subsidence since the Late Cretaceous (Bond, 1978; DiCaprio et al., 2009; Russell and Gurnis, 1994) and a progressively increasing northeast downward tilt since the Eocene (DiCaprio et al., 2009; Sandiford, 2007; Sandiford et al., 2009). Broad-scale subsidence of the continent is inferred from the difference between relative and global sea level. As global sea level fell starting in the Cretaceous (Haq et al., 1987; Haq and Al-Qahtani, 2005; Miller et al., 2005), Australia became increasingly flooded. The preserved marine facies imply that the entire continent has experienced an ∼200 m downward shift since the Cretaceous (Bond, 1978; DiCaprio et al., 2009; Russell and Gurnis, 1994; Veevers, 2000). In addition to the broad-scale subsidence, differential vertical motion of the continent since the Cenozoic inferred from reconstructed shorelines shows that the Australian continent was progressively tilted down toward the northeast by around 300 m since the Eocene (DiCaprio et al., 2009) (Fig.1).

The downward shift and progressive tilting of Australia remain unexplained geodynamically. Here, for the first time, we attempt to interpret these first-order features in the stratigraphic record with geodynamic models. Moreover, to overcome the limitations of earlier geodynamic models (Gurnis et al., 1998; Gurnis and Müller, 2003), we include a more realistic lithosphere through the assimilation of tectonic boundary conditions. Our geodynamic models include the history of plate motions, changing plate geometries and modeled ages and thermal structures of the subducting lithosphere to investigate the time-dependent dynamically driven topography of the southwest Pacific and Australia since the Eocene. Surprisingly, we find that the subduction control on the vertical motions of Australia is insufficient to explain the continentwide subsidence observed in paleoshoreline analysis.

We used continuously changing plate velocities and plate margins at 1m.y. increments as surface velocity boundary conditions (Gurnis et al., 2011). Ocean floor paleo-age grids (Müller et al., 2008) constrain the lithospheric temperature profile. We embedded a high-resolution regional model within a low-resolution global model using nested CitcomS solvers (Tan et al., 2006) and incompressible flow. The nested solver allows us to compute the model in a high-resolution regional grid with a flow-through boundary condition, which avoids the artifact in dynamic topography caused by return flow within a confined box. The regional model has more than four times the lateral resolution of the global models, resolving features at 64 km in longitude and 40 km in latitude. The radial resolution of the regional model is 33 km to a depth of 2250 km. We used active tracers within the continental lithosphere to mimic chemical buoyancy and within the mantle wedge to lower the viscosity (Manea and Gurnis, 2007). For more information about the method and data assimilation, see GSA Data Repository material¹ and DiCaprio (2009).

Our models have a Newtonian viscosity that is dependent on temperature, depth, composition, and position. The mantle is divided into four layers: lithosphere (0–100 km), asthenosphere (100–410 km), transition zone (410–670 km), and lower mantle (670–2880 km). All models have the same radial distribution of viscosity (Table1); however, models with alternative layered viscosities were examined further in DiCaprio (2009). We varied both the strength of the asthenosphere and lower mantle. Our preferred models have a weak asthenosphere (2 × 10²⁰ Pa·s) and strong lower mantle. The weak asthenosphere allows subducted slabs to be more easily distributed laterally within the upper mantle beneath the Australian northeast margin and produces a broader topographic signal compared to those with a higher-viscosity upper mantle. Models with a high-viscosity lower mantle (1 × 10²³ Pa·s) impede the descent of slabs into the lower mantle. These models retain more cool material in the upper mantle and produce a smaller-amplitude but broad-scale length of dynamic topography. Here, we show the effect of changes in the gradient of the Clapeyron slope for the phase changes bounding the mantle transition zones, which we altered for a range of reasonable values based on mineral physics experiments (Table 2), as summarized in Billen (2008).

Our models show that since the Eocene, subducted material has accumulated beneath the northeastern margin of Australia and the Tasman Sea (Fig.2). In the model, slabs presently beneath the northeast Australian margin are the product of Melanesian subduction, while those beneath the Tasman Sea are remnants of Loyalty-Tonga-Kermadec subduction. This result is consistent with fast-velocity perturbations observed beneath the NE margin in tomography (Hall and Spakman, 2003; Ritsema et al., 1999). Between 50 and 32 Ma, slabs from Melanesian and Loyalty subduction descended through the upper mantle and were laid flat within the transition zone more than 2000 km away from the trench. At this time, dynamic topography affecting the Australian continent was negligible, while dynamic topography amplitudes of several hundred meters of subsidence were concentrated within the backarc regions (Fig.2). However, between 32 and 10 Ma, a long-wavelength component of dynamic topography developed as Australia began to drift over the slabs accumulated within the transition zone. This caused subsidence of up to 100 m along the northeastern Australian margin between 32 and 10 Ma. Between 10 and 2 Ma, the slabs for which descent was partially impeded by the phase transition at 660 km depth started to pass through the 660 km phase boundary into the lower mantle. This resulted in up to 200 m of total dynamic subsidence on the north and northeastern Australian margin, corresponding to the timing of northeastern margin reef demise (DiCaprio et al., 2010; Isern et al., 2002).

Since the continent experienced little vertical disturbance at 50 Ma (Fig.2), its topography at 50 Ma is a good proxy of isostatic topography. Differential topography is computed as the change relative to the initial topography at 50 Ma, a quantity that can be compared directly against the anomalous subsidence and tilting of the continent recovered through an analysis of paleoshorelines (DiCaprio et al., 2009). A SW-NE cross section shows that the slope of the modeled subsidence across the continent poorly matches the slope from paleoshoreline analysis since the Miocene (Fig. 3). Although, the modeled subsidence is qualitatively consistent with the trend observed by paleoshoreline analysis, the modeled subsidence is much too confined to the north. Both models and paleoshoreline analysis show that since the Cenozoic, Australia has subsided in the northeast, and the magnitude of this subsidence has increased toward the present.

The model that best approximates the subsidence estimated by paleoshoreline analysis has a steep Clapeyron slope (4 MPa/K) at the 410 km discontinuity and a shallow Clapeyron slope (2 MPa/K) at the 660 km discontinuity (Fig.3). The shallow Clapeyron slope at the 660 km phase change allows slab material to descend into the lower mantle relatively easily and produces a larger signal of surface subsidence in the northeast.

Our models produce insufficient surface subsidence since the Miocene, ∼100 m smaller than that inferred from paleoshoreline analysis. This implies that the continent requires an additional downward shift over a much longer wavelength since the Oligocene. The downward shift of the Australian continent may be related to its motion northward away from putatively hot or chemically buoyant mantle beneath Antarctica.

The mantle beneath Cenozoic Antarctica may have been relatively hot due to its connection to Gondwanaland. This is inferred from the abundance of magmatism before and during breakup (Kent, 1991; Storey et al., 1995). The presence of a large-scale upwelling beneath Antarctica during the Cenozoic and today is consistent with the distribution and geochemistry of igneous rocks (Behrendt et al., 1991; Finn et al., 2005; LeMasurier and Landis, 1996; Sutherland et al., 2010). In addition, global tomography models reveal low-velocity perturbations in the upper mantle and transition zone beneath the West Antarctic margin (Grand, 2002; Gu et al., 2001; Masters et al., 2000; Ritsema and Heijst, 2000). Recently, geodynamic models have shown that the Campbell Plateau subsided by ∼1km as it moved away from a buoyant mantle anomaly presently located beneath Antarctica (Spasojevic et al., 2010; Sutherland et al., 2010). These observations and interpretations suggest that Australia has likely also moved away from a dynamic topography high since 50 Ma.

Motivated by these observations, we modified a model named M50_1 by prescribing a 5% temperature difference (75 °C hotter) in the mantle beneath Australia and Antarctica. The buoyancy associated with this hot mantle causes a topographic high beneath both continents at 50 Ma. However, as the spreading rate at the South East Indian Ridge increases and accelerates Australia northward, the Australian continent moves away from this topographic high (Fig. DR1 [see ^{footnote 1}]). This motion causes vertical subsidence of the whole continent, in addition to a downward tilt toward the northeast, as the continent overrides the Melanesian subducted slabs. The summation of these two dynamic forces matches the predicted tilt of the continent from paleoshoreline analysis (Fig. 3). The buoyant mantle produces the additional broad-scale subsidence required to fit the paleoshoreline analysis. This can be seen by computing a residual since the Miocene (Fig. 4C) between the modeled differential topography (Fig. 4A) and the anomalous topography from paleoshorelines (Fig. 4B). However, we would like to note that the mantle buoyancy, here modeled as simply a thermal effect, is likely due to a combination of mantle hydration and other geochemical and thermal heterogeneities that may be associated with the long-lived eastern Gondwanaland slab graveyard (Spasojevic et al., 2010).

Dynamic topography from models of the evolution of subducted slabs since the Eocene is in disagreement with geologic observations of continent-scale tilting. Models with the negative thermal buoyancy associated with slabs between Melanesia and the proto–Tonga-Kermadec subduction systems only recover the tilt of the continent, while the vertical motion of the continent is underestimated by up to 200 m. By including 75 °C (5%) hotter than average mantle beneath the Australian and Antarctic continents, our models match both the tilt and vertical displacement of the Australian continent. The agreement of the models with Australia's paleoshoreline evolution suggests that the mantle anomaly was buoyant, and it may have been the result of mantle hydration or some combination of thermal and geochemical heterogeneities. During the Cenozoic, as spreading along the South East Indian Ridge accelerated Australia away from the buoyant Antarctic mantle and toward the subduction zones in Melanesia, the continent both subsided in bulk and tilted to the NE.

We are thankful to Leif Strand from Computational Infrastructure for Geodynamics for his assistance with CitcomS. The Australian Research Council Australian Postgraduate Award administered by the University of Sydney, with additional support from National Science Foundation grant EAR-0810303, financially supported Lydia DiCaprio. This work represents contribution 10022 of the Division of Geological and Planetary Sciences, California Institute of Technology.