Global patterns in Earth’s dynamic topography since the Jurassic

. We evaluate the spatial and temporal evolution of Earth’s long-wavelength surface dynamic topography since the Jurassic using a series of high-resolution global mantle convection models. These models are Earth-like in terms of convective vigour, thermal structure, surface heat-ﬂux and the geographic distribution of heterogeneity. The models generate a degree-2-dominated spectrum of dynamic topography with negative amplitudes above subducted slabs (i.e. circum-Paciﬁc regions and southern Eurasia) and positive amplitudes elsewhere (i.e. Africa, north-western Eurasia and the central Paciﬁc). Model predictions are compared with published observations and subsidence patterns from well data, both globally and for the Australian and southern African regions. We ﬁnd that our models reproduce the long-wavelength component of these observations, although observed smaller-scale variations are not reproduced. We subsequently deﬁne “geodynamic rules” for how different surface tectonic settings are affected by mantle processes: (i) locations in the vicinity of a subduction zone show large negative dynamic topography amplitudes; (ii) regions far away from convergent margins feature long-term positive dynamic topography; and (iii) rapid variations in dynamic support occur along the margins of overriding plates (e.g. the western US) and at points located on a plate that rapidly approaches a subduction zone (e.g. India and the Arabia Peninsula). Our models provide a predictive quantitative framework linking mantle convection with plate tectonics and sedimentary basin evolution, thus improving our understanding of how subduction and mantle convection affect the spatio-temporal evolution of basin architecture.


Introduction
At short spatial scales, Earth's topography depends on interactions between crustal structures, tectonic deformation and surface processes.However, at length scales exceeding the flexural strength of the lithosphere, topography is an expression of the force balance acting at Earth's free surface.This balance has (i) an isostatic component, which depends on the composition and thickness of the crust and lithospheric mantle, and (ii) a dynamic component, which is the deflection of Earth's surface resulting from vertical stresses generated by convection in the underlying mantle.Their relative importance varies according to the tectonic setting.Most orogens show patterns of thickened continental crust and hence isostatically compensated topography (e.g.Frisch et al., 2011), while the anomalous elevation of the stable southern African craton has been attributed to dynamic topography generated by deep-mantle upwellings (e.g.Lithgow-Bertelloni and Silver, 1998;Gurnis et al., 2000).
Dynamic topography generates significant surface deflections in both continental and oceanic regions (e.g.Ricard et al., 2006).The connection to sea level and continental flooding is direct: as continents migrate over mantle upwellings (regions of positive dynamic topography) and downwellings (negative dynamic topography), vertical motions may be sufficiently large to control the emergence and submergence of continents (e.g.Mitrovica et al., 1989), thus respectively generating surface depressions in which sediments accumulate or erosional surfaces.Dynamic topography is identified in Earth's present-day relief via deviations from isostasy, yielding the so-called residual topography field (e.g.Kaban et al., 1999;Doin and Fleitout, 2000;Crosby et al., 2006;Winterbourne et al., 2009;Flament et al., 2013;Hoggard et al., 2016).It is recorded via anomalous vertical motion events in the basin infill (e.g.Mitrovica et al., 1989;Heine et al., 2008) or long-wavelength uplift in non-basinal areas, as seen in plateau exhumation and geomorphic features (e.g.Roberts and White, 2010).It is, therefore, a fundamental component of Earth's topographic expression and an improved understanding of the history of Earth's surface topography can only be achieved by successfully accounting for its spatial and temporal variations.
There are two principal ways to retrieve the component of dynamic support: (i) the "top-down" inverse approach, which is based on geophysical and geological observations and subtracts isostatic topography from observed surface elevations, and (ii) the "bottom-up" strategy, which infers dynamic support by determining mantle flow from Earth's internal buoyancy and viscosity distribution.In this study, we follow the bottom-up strategy and discuss our findings in relation to results from the top-down approach (e.g.Panasyuk and Hager, 2000;Kaban et al., 2003Kaban et al., , 2004;;Steinberger, 2007;Heine et al., 2010;Flament et al., 2013;Hoggard et al., 2016).Discrepancies between these two basic strategies lead to much discussion in the literature (e.g.Flament et al., 2013;Molnar et al., 2015;Colli et al., 2016;Hoggard et al., 2016).Molnar et al. (2015) argue that only small free-air gravity anomalies (< 30-50 mGal) are caused by long-wavelength topography, interpreting most of the signal as isostatic.They suggest that a maximum of 300 m of dynamic support is indicated by long-wavelength gravity anomalies, which is in contrast to many published geodynamic models (e.g.Flament et al., 2013).However, Colli et al. (2016) demonstrate that mantle convection with a depthdependent viscosity generates small perturbations (geoid 37 m; gravity 8 mGal), considerable dynamic topography (1.1 km) and a gravitational admittance of mostly less than 10 mGal km −1 .This supports the idea that relatively small long-wavelength free-air gravity anomalies do not preclude dynamic topography at amplitudes substantially larger than 300 m.It also helps to underpin the results of Hoggard et al. (2016) on variations of 1 km in dynamically induced topography inferred from observations.Numerous bottom-up studies have been undertaken in recent years in an attempt to constrain the amplitude, spatial pattern and time dependence of dynamic topography using numerical modelling.These studies examined the temporal evolution of dynamic support by forward modelling (e.g.Ricard et al., 1993;Zhang et al., 2012;Flament et al., 2013), backward advection (e.g.Conrad and Gurnis, 2003;Moucha et al., 2008;Heine et al., 2010), adjoint schemes (e.g. Bunge et al., 2003;Liu et al., 2008;Spasojevic et al., 2009;Colli et al., 2017) or a back-and-forth iterative method (Glišović andForte, 2016, 2017).Present-day dynamic topography is discussed on a global scale by Ricard et al. (1993), Steinberger (2007), Zhang et al. (2012) and Flament et al. (2013), but relatively few studies (e.g.Zhang et al., 2012;Flament et al., 2013) which use the forward modelling approach have examined the temporal evolution of dynamic support globally.
In this paper, we also use numerical forward modelling to compute the transient evolution of dynamic topography.The distribution of heterogeneity in our models is generated by solving the equations relevant to mantle convection inside a global spherical shell, with plate velocities derived from a global plate kinematic model (Seton et al., 2012, with modifications by Shephard et al., 2013) and featuring topological plate boundaries extending back to 200 Ma imposed as a surface boundary condition.We compare our predicted present-day dynamic topography to previous numerical studies before discussing the temporal evolution of dynamic topography within selected regions.Finally, we propose a categorisation of regions on Earth's surface in terms of their exposure to the effects of dynamic support over the past 200 Myr, which is corroborated by an automated classification technique.

Model set-up
We use the global mantle convection code TERRA (e.g.Baumgardner, 1985;Bunge et al., 1997;Davies and Davies, 2009;Wolstencroft et al., 2009;Davies et al., 2013;Wolstencroft and Davies, 2017), which solves the Stokes and energy equations relevant to mantle convection inside a spherical shell, with appropriate material properties and boundary conditions (Fig. 1).Calculations are performed on a mesh with ∼ 80 million nodal points (lateral nodal spacing of ∼ 28 km at the surface and ∼ 14 km at the core-mantle boundary (CMB); radial spacing of ∼ 22 km across the entire domain), thus providing the resolution necessary to explore mantle flow at Earth-like convective vigour.Our models incorporate compressibility via the anelastic liquid approximation (ALA), with radial reference values represented through a Murnaghan equation of state (Murnaghan, 1944).Isothermal boundary conditions are specified at the surface (273 K) and CMB (models M1 and M3 4000 K, model M2 3500 K), with the mantle also heated internally at roughly chondritic rates.The presentday surface heat flux of model M1 is 34 TW, which is consistent with current estimates for Earth's surface heat flux (∼ 47 TW) minus the contributions from crustal radioactivity (Davies and Davies, 2010).Of this 34 TW, ∼ 30 % is transferred across the CMB with the remainder introduced through internal heating.
A free-slip boundary condition is specified at the CMB, whilst surface velocities are assimilated via 200 Myr of plate motion histories (Seton et al., 2012;Shephard et al., 2013) generated by the plate tectonic reconstruction software GPlates (Boyden et al., 2011;Fig. 1) at discrete 1 Myr intervals.In our reference model (M1), we prescribe a radial viscosity profile with four layers, namely the lithosphere (10 23 Pa s), the sub-lithospheric upper mantle (10 20 Pa s), the transition zone (10 21 Pa s) and the lower mantle (10 23 Pa s).Viscosity also varies with temperature T , following the relation: where T * represents the non-dimensionalised temperature.
For initial conditions, we approximate the unknown Jurassic mantle heterogeneity structure by applying a 70 Myr preconditioning stage during which global plate configurations are fixed to the oldest available reconstruction at 200 Ma.
From here, models run forward for 200 Myr toward the present day.Note that for simplicity our models assume no chemical heterogeneity.The distribution of heterogeneity is principally dictated by the imposed plate motion histories, the material properties of the model and the heating mode.Purely thermal models, as in the present study, have previously been shown to provide a good match to seismic observations (Schuberth et al., 2009;Davies et al., 2012Davies et al., , 2015) ) and dynamic topography (Colli et al., 2017).Key model parameters are listed in Table 1 with reference values also illustrated in Fig. 2. Three separate cases are examined with the variations between these cases summarised in Table 2. Dynamic topography is computed from normal stresses at the model surface using the following formula: where σ rr is the normal stress, ρ is the density difference between surface and air and g is the acceleration due to gravity.This is likely to yield elevations that are too high, given that we are restricted in our choice of the asthenospheric viscosity, which we limit to a minimum of 10 20 Pa s.Estimates for the viscosity of Earth's asthenosphere are between 10 19 and 10 21 Pa s (e.g.Mitrovica and Forte, 2004;Steinberger and Calderwood, 2006;Petersen et al., 2010).Whilst we cannot compensate for the weaker coupling of tectonic plates to the mantle that this entails, we can account for the increased elevations predicted in our models through the introduction of a scaling.Thus, we assume simple linear scaling and introduce a modified formula for dynamic topography calculations: where η 0 = 3 × 10 19 Pa s is an estimate of the asthenospheric viscosity, while η ref is the model reference viscosity (see Table 2).This is less than the lowest reference viscosity used in our model (10 20 Pa s) and subsequently acts to scale our predicted dynamic topography to the observed amplitudes.
The scaling factor thus encapsulates our assumption that the stresses scale linearly with the ratio of asthenosphere to model viscosity found at the isostatic compensation depth of the lithosphere.Note that to calculate dynamic topography, we focus on long-wavelength subduction-driven global flow below the lithosphere, as this is the contribution to dynamic topography that does not depend on the specific structure of the lithosphere and the depth of the lithosphereasthenosphere boundary.Following these considerations, we replace the model thermal structure above 225 km of depth with the constant median mantle temperature of 2150 K.We subsequently re-run the model for one time step employing a rigid (zero-slip) boundary condition at the model surface and use the resulting surface normal stresses to evaluate dynamic topography.The depth above which we applied the blanking procedure has been chosen such that all lithospheric structures are encompassed.However, we varied this limiting depth between 300 and 150 km and the impact on the results was negligible.We also varied the blanking temperature of 2150 K, but since the chosen temperature is applied to the entire spherical shell, it is always neutral with respect to the mantle buoyancy structure and the amplitude of dynamic topography does not depend on the specific temperature choice.Similar approaches to discard the lithospheric buoyancy structure from the dynamic topography signal have been applied by Lithgow-Bertelloni and Silver (1998), Steinberger (2007), Conrad and Husson (2009) and Flament et al. (2013) for structures above 325, 220, 300 and 350 km of depth, respectively.

Model robustness
To assess the robustness of our results, we vary two key parameters in our reference model, M1, and compare the temporal evolution of dynamic topography at selected locations.These synthetic sample points have been picked to illustrate the effect on regions where either the influence of dynamic support at present day or through time is debated, or the geological setting of which prompts that question.Model M2 employs the same viscosity profile as M1 (Fig. 2 and Table 2) but differs in CMB temperature; a temperature of 3500 K is prescribed instead of 4000 K. Model M3, on the other hand, utilises the same CMB temperature of 4000 K, but its viscosity profile differs from model M1: the lower mantle viscosity is reduced from 1 × 10 23 to 3 × 10 22 Pa s.
The computed dynamic topography is not strongly affected by these variations, as demonstrated by the elevation histories of chosen sample locations in Fig. 3. Decreasing the CMB temperature has a negligible effect, with models M1 and M2 predicting very similar dynamic topography amplitudes and trends.Decreasing the lower mantle viscosity has a more pronounced influence, as slab penetration into the lower mantle is facilitated and slab-sinking rates are increased.This has two principal effects: (i) due to viscous coupling, faster-sinking slabs increase the amplitude of negative dynamic topography, and (ii) the suction effect of a faster-sinking detached slab fades quicker after subduction slows down or ceases.Nevertheless, the firstorder characteristics of all three models (Fig. 3) are in good agreement, indicating that long-wavelength dynamic topography is principally controlled by the plate tectonic history and the location of subduction zones.Accordingly, for the remainder of this paper, we focus solely on our reference model, M1.
3 Model results and implications 3.1 Predicted mantle structure Snapshots of the mantle structure predicted in model M1 are presented in Fig. 1, illustrating that the upper mantle planform is dominated by strong downwellings in regions of present-day plate convergence.In the mid-mantle, cold downwellings are prominent beneath North America (the subducted Farallon slab) and southern Eurasia (the former Neo-Tethys ocean), whilst remnants of older subduction (those encircling the former supercontinent Pangaea) are visible above the CMB.In our model M1, slabs sink at rates of 1.5 cm year −1 , which is comparable to those www.solid-earth.net/8/899/2017/Solid Earth, 8, 899-919, 2017 inferred from seismic tomography (e.g.Li. et al., 2008;Simmons et al., 2010).These downwellings modulate the location of upwelling flow, such that it becomes concentrated beneath Africa and the Pacific (e.g.McNamara and Zhong, 2005;Davies et al., 2012Davies et al., , 2015)).The Pacific anomaly is approximately circular, whilst the African anomaly is a NW-SE trending structure, which to the north curves eastward under Europe and to the south extends into the Indian Ocean.Note that due to (i) the short evolution times of our models and (ii) the limited temperature dependence of viscosity, plumes play a negligible role in the convective planform.Accordingly, the regions of upwelling flow beneath Africa and the Pacific discussed above principally represent passive return flow from subduction.Nevertheless, we also examined models with different initial conditions, under which plumes were able to form.We found that dynamic topography is modified in regions of active upwelling.However, these dynamic topography highs constitute shorter wavelength features (Moucha and Forte, 2011;Hoggard et al., 2016) and, on a global scale, are less significant than the role of subducting slabs; this is why we do not examine plumerelated dynamic topography in this study.

Global and regional dynamic topography
The spatio-temporal evolution of dynamic topography from case M1 is illustrated in Fig. 4a in the mantle reference frame  and in Fig. 4b in a plate frame of reference.Note the strong correspondence of the upper mantle thermal heterogeneity distribution with negative dynamic topography in regions of present-day and former subduction, and positive dynamic topography elsewhere.Before we further interpret our results we compare them with previous geodynamic modelling studies and with selected observations on a regional scale.We evaluate our model against well data at Australian locations, which have experienced two marked periods of anomalous subsidence and uplift since the Jurassic.Finally we consider observations that constrain the uplift history of southern Africa.Due to its particular geological setting, southern Africa has been the target of many studies on constraining its uplift history and assessing to what extent and when it experienced dynamic support.

Global models of dynamic topography
To place our results in the context of previous studies, we compare our predictions of present-day dynamic topography ( www.solid-earth.net/8/899/2017/Solid Earth, 8, 899-919, 2017  different tomography models after obtaining mantle density through the thermodynamic modelling of various chemical compositions in the mineralogical assemblage.This yields a peak-to-peak dynamic topography amplitude exceeding 3 km for all models.Other studies suggest that amplitudes for long-wavelength (10 3 km) dynamic topography should be smaller than ±1 km (Wheeler and White, 2002;Molnar et al., 2015;Forte et al., 2015;Hoggard et al., 2016).In contrast to some previous models, our model M1 agrees well with this suggestion, mapping deep mantle return flow above the Pacific and African LLSVPs at amplitudes not exceeding ±1 km.

Australia
The tectonic stability of the Australian continent and its relatively large displacement relative to the underlying mantle since the Cretaceous makes it an ideal location for investigating mantle-induced topography changes.It is the only continent to have undergone two distinct episodes of dynamic topography during the last 200 Myr: (i) a pronounced west-east-dominated warping during the middle and late Jurassic and early Cretaceous (Gurnis et al., 1998) caused by subduction along its eastern margin as the Phoenix plate subducted westward under East Gondwana, and (ii) a south-north tilt due to subduction of the Australian plate underneath the Eurasian and Pacific plates (Fig. 6) that started in the Miocene and continued through to the present day (Sandiford, 2007).The effect of transient mantleconvection-induced vertical Australian motions has been the subject of many observation-and model-based research efforts (Veevers, 1984;DiCaprio et al., 2009;Heine et al., 2010;Czarnota et al., 2014).
To ground-truth our predicted dynamic topography we focus on tectonic uplift and subsidence histories from four selected Australian exploration wells located in opposite corners of the Australian continent.The well data were retrieved from the Geoscience Australia PIMS system online (http:// dbforms.ga.gov.au/pls/www/npm.pims_web.search;Fig. 7).
For the well on the North West Shelf we find good agreement of the long-term Cenozoic elevation history (Fig. 7b; Müller et al., 2000).The modelled vertical motion history shows that the area was subjected to gradually increasing subsidence during the Cenozoic, which is in good agreement with observations (Czarnota et al., 2013).Also for Figure 10.Volume of subducted oceanic lithosphere (in km 3 ) for the western Tethyan region for the past 200 Myr.The amount of material was calculated based on paleo-age grids with a 6 arcmin grid cell size for the subducted oceanic lithosphere (see Müller et al., 2008, for methodology) and convergence velocities derived from the input plate kinematic model.We assume a simple halfspace cooling model to compute the thickness of the oceanic lithosphere for ages < 69.9 Myr as z = 11.36 √ age, where age is the age of the oceanic lithosphere.Older lithosphere is set to a default thickness of 95 km (Kanamori and Press, 1970).Volumes are derived by multiplying thickness with computed trench-normal convergence in 1 Myr time steps.Grid cell size for the subduction volume is 30 arcmin.Map background is ETOPO1 topography (Amante and Eakins, 2009).
the remaining wells (Fig. 7a, c and d), our modelled vertical motion histories follow the backstripped subsidence patterns generally quite well.
The Gandara 1 well in the north-western part of Australia (Fig. 7b) shows a distinctly different temporal subsidence behaviour from the most eastern one in the Surat Basin, Cecil Plains West 1 (Fig. 7d).We explain this through the evolution of dynamic topography shown in Fig. 6.The former well is not affected by subducted, cold, dense slab material from the Proto-Western Pacific/East Gondwana subduction zone.When, from around 70 Ma onward, Australia starts progressively moving north, this area experiences the increasing influence of downward mantle flow related to the subduction of the Australian lithosphere along the Indonesian subduction zone system.The Cecil Plains well, on the other hand, mostly displays the Jurassic and early Cretaceous East Gondwana subduction setting, while Australia's Cenozoic tilting only has a moderate effect: this location indicates phases of accelerating subsidence until about the latest Jurassic, followed by a phase of non-deposition and erosion and a recent phase of gentle subsidence.Our models can reproduce the initial Jurassic phase of subsidence.We attribute the time of the depositional hiatus (∼ 155 to 25 Ma) to the rebounding of the west-east tilt, while the renewed phase of subsidence seems to be related to the Cenozoic warping of the Australian continent.
For the remaining bore holes, the modelled Duyken 1 and Mallabie 1 dynamic topography histories can be explained along analogous lines of argument (Fig. 6).Both reflect Cretaceous subsidence, with the former being more pronounced as the subducted material progressively extends farther west in the north of Australia.A decrease in subsidence can be seen when that influence fades, just to be affected even more profoundly later when the Australian north is drawn down by the negative dynamic topography due to sinking slabs in the Indonesian subduction zone.Keeping in mind that dynamic topography varies only slowly over time, we observe a good correlation between modelpredicted vertical motions and reconstructed vertical motions from well data.Hence, these data corroborate the dynamic topography amplitudes from our study and are in line with other studies (e.g.DiCaprio et al., 2011;Heine et al, 2010;Czarnota et al., 2013).

Southern Africa
Southern Africa's unusually high elevation relative to other continents, despite being surrounded by passive margins, makes it an area of particular interest.Braun et al. (2014) analysed sedimentary flux data, which confirmed a major phase of erosion on the southern African plateau in the late Cretaceous.By examining paleoprecipitation indicators across the plateau, they ruled out climate variability as the main driver of increased erosional activity.Thus, Braun et al. (2014) concluded that the continent migrated over a source of positive dynamic topography at that time.Tinker et al. (2008) quantified exhumation using apatite fission track dating.While denudation rates during the Cenozoic indicate that mantle processes did not primarily cause the uplift of the region, the authors suggest that mantle buoyancy during the Cretaceous was responsible for the enhanced erosion across southern Africa, which in turn generated up to 3 km of exhumation (Tinker et al., 2008).Stanley et al. (2013) performed a detailed thermochronometry study across the southern African plateau by employing (U-Th)/He dating.They suggest that a mid-Cretaceous buoyancy increase in the lithospheric mantle caused coeval denudation.Our model predicts dynamic support equivalent to ∼ 200 m underneath southern Africa gradually developing in the Cretaceous period (Fig. 8), which is in line with previous modelling studies (Zhang et al., 2012;Flament et al., 2014).However, it does not reproduce a significant post-30 Ma uplift phase as revealed by the inversion of river profiles (Paul et al., 2014), which may reflect the lack of active upwellings in our convection model or convective processes at (sub-) lithospheric levels.

Modes of long-term dynamic surface topography evolution
In this section we use our model results to establish a set of "geodynamic patterns" for how different continental regions are affected by underlying mantle processes.We define three categories that are exemplified by selected locations using two alternative approaches.First, we plot time series of their dynamic topography and visually classify them into these categories based on the overall temporal trends in elevation.Secondly, we apply a k-means-based cluster analysis to the time series of dynamic topography as an independent means to validate our interpretation and explore the consequences of classifying the time series into different numbers of groups.

Category I -Topographic stable areas
The first category is comprised of locations in topographically stable regions.These areas experience only positive dynamic support over the 200 Myr time period examined, generating stable positive topography of a few hundred metres.As can be seen from Fig. 9, these regions lie on a curved band that stretches from the Eurasian to the African plate and extends into the Somali plate, encompassing locations in Africa, Madagascar, Western and Northern Europe and Western Siberia.This category generally consists of regions with no adjacent subduction zone.Where points are located adjacent to a subduction zone (for example the Rhenish Massif next to the western Tethyan subduction zones; Fig. 9h), there is insufficient cold, dense slab material in the mantle to cause negative dynamic www.solid-earth.net/8/899/2017/Solid Earth, 8, 899-919, 2017 topography, as evidenced by a computation of subducted oceanic lithosphere volumes (Fig. 10) for the past 200 Myr.The total volume of subducted material in the western Tethyan and present-day Mediterranean region is on the order of 10 000-15 000 km 3 per 30 arcmin grid cell, yielding very low average volume fluxes of around 50-75 km 3 Myr −1 .

Category II -Dynamically subsiding areas
Locations that fall into the second category experience a long-term continuous decrease in dynamic support.
Examples include East and Western India, the Arabian Peninsula, north-eastern Brazil, the north-eastern United States, north-eastern Canada and Western Australia (Fig. 11a  and b).These regions subside as they approach sites of active or recent subduction.Note that Australia is a special case in which two subduction zones exert a sequential influence (Sect.3.2.2).First the East Gondwana and subsequently the Indonesian subduction zone give rise to a decrease in elevation, as described above.

Category III -Fluctuating areas
Locations in this category (Fig. 12) exhibit a nonmonotonic temporal evolution with mostly negative dynamic topography.This is likely caused by two factors: (i) the influence of multiple subduction zones or (ii) the influence of a single subduction zone that exhibits strongly timedependent characteristics.
Multiple subduction zones affected the dynamic topography evolution of Australia, as discussed above: two lows in dynamic subsidence can be identified when the particular location is influenced by one of these.The Chukchi Sea region (Fig. 12a) is a second example where three local minima in dynamic topography correlate with sequential activity in Arctic subduction zones: the Koni-Taigonos subduction zone at ∼ 200 Ma, the Koyukuk and Nutesyn subduction between 140 and ∼ 125 Ma and the Alaskan subduction zone since ∼ 110 Ma (Shephard et al., 2013; see Fig. 4a and the corresponding Supplement).A third example is the Central Siberian Plateau (Fig. 12g), which was first affected by the East Gondwana subduction A key time-dependent parameter of a subduction zone is the subduction velocity.Rapid convergence tends to introduce more cold, dense material into the mantle, which leads to a more pronounced negative dynamic topography.To isolate the role of subduction velocity, we compare dynamic topography time series at three different subduction zones with the orthogonal component of the relative velocity between the subducting and overriding plate at the closest subduction zone (Fig. 13).For instance, the Chukchi Sea location exhibits a significant increase in negative dynamic topography about 0 m at 145 Ma to about −500 m at 120 Ma.This modelled dynamic subsidence phase is followed by a short period of rebound (+150 m over 15 Myr) before the area is again modelled as experiencing dynamic subsidence of more than 500 m between 110 and 80 Ma.Here, we can correlate the modelled trench-normal convergence velocities from our input plate kinematic model (Seton et al., 2012;Shephard et al., 2013) directly with predicted dynamic topography evolution over time.Similarly, accelerating convergence from 200 Ma until ca.150 Ma at the Western Interior Seaway sample point corresponds to decreasing dynamic topography.This is followed by a short uplift period due to slow convergence before elevation decreases again, driven by the increased convergence velocity.An analogous situation occurs in the South American Amazonian region, where the general trend of orthogonal convergence correlates well with the evolution of dynamic topography.It should be noted that the slabs in our models are generally weaker than slabs on Earth due to our simplistic representation of a primarily depth-and temperaturedependent rheology (see Garel et al., 2014, for comparison).As a consequence, upon slab descent, slabs often tear and break, which leads to rapid, localised changes in regional dynamic topography (e.g.Fig. 12c-e; h-l).Similar tears or breaks occur due to abrupt changes in plate motions.The closer a point lies to a subduction zone where this occurs, the larger the associated amplitude of dynamic topography change.

Cluster analysis
In the previous section, we isolated three characteristic groups which represent distinct trends in vertical motion.This classification was based on a simultaneous consideration of absolute elevations, changes in elevation and the geodynamic context.Here we use cluster analysis, based on the k-means algorithm (e.g.Lekic and Romanowicz, 2011), as an independent means to validate our classifications (see the Supplement for a detailed description).To this end, we search for clusters of similarity within our dynamic topography time series at individual surface points using a grid of nearly 2 × 10 5 equidistantly spaced nodes on the sphere.This provides a resolution of ∼ 28 km on the surface sampled at 5 Myr intervals in time.The cluster analysis was applied to the topography time series from 150 Ma to present; this excludes the first 50 Myr of the model run in which our dynamic topography predictions will be strongly sensitive to the initial condition.The k-means algorithm partitions a given dataset into k clusters such that the sum of the squares of the deviation from the cluster mean is minimised.The number of clusters k is an input parameter for the algorithm.Note that this approach only evaluates elevations at the particular grid point at a certain time and does not account for relative variations (i.e.uplift or subsidence), nor does it take into account the history of elevations.For a detailed description of the kmeans algorithm and the results for all clusters (k), see the Supplement.
We vary k between 2 and 6, which allows us to highlight a variety of spatio-temporal patterns in dynamic topography.
Here we discuss k = 5 because this is the minimal number of clusters that reflects the underlying geodynamics to a first order (Fig. 14).Cluster 5 can be directly associated with our Category I ("topographic stable regions"), while Cluster 4 is an approximate representation of our Category II regions ("dynamically subsiding areas").To first order, the map in Fig. 14 mirrors our dynamic topography maps in that the topographically high regions (dark blue areas) are concentrated in a region that stretches from Russia over Europe to Africa, mostly far away from subduction zones.Cluster 4 (light blue) corresponds to subsiding areas which approach subduction zones.Also, North India falls into this category (as opposed to k < 5; see the Supplement), which describes the history of the subcontinent starting out as a part of the highly elevated Pangaea supercontinent to the plunging of its northern edge as it collides with Asia (Fig. 11a, b, g).The Arabian Peninsula, however, has been determined not to reside in this group, but in the "topographic high" group instead; on average for most of the observed time period, this location displays positive elevations (Fig. 11a, l).The same holds true for the borderline case of the Cuvier Abyssal Plain (Fig. 11a, k).Category III ("fluctuating areas") disintegrates into three distinct clusters (red, orange and green): Cluster 1 with points that remain at negative dynamic topography throughout the entire model run, and Clusters 2 and 3 with locations that experience an initial phase of long-term dynamic subsidence followed by a subsequent uplift.Higher cluster numbers k serve to further partition our categories, yet do not reveal any additional information (see the Supplement).To summarise, the k-means scheme provides an independent instrument for the quantitative analysis of our model data and we find good correspondence to our classification.

Conclusions
We have demonstrated that the long-wavelength influence of dynamic topography in different continental regions can be classified via their spatio-temporal vicinity to a subduction zone.(i) Locations that are far away from subduction zones are generally stable and dynamically supported; these lie on a curved band stretching from the Eurasian plate to the African plate and extending into the Somali plate.(ii) Locations that move toward subduction zones dynamically subside over time, and (iii) locations that are being influenced by numerous subduction zones or a subduction zone that shows strong time dependence (i.e.abrupt changes in plate motions of both the overriding and subducting plate) will display an oscillatory history of dynamic topography.While this classification is based on a geodynamic understanding, it is confirmed through an independent approach based on a k-means cluster analysis.As geosciences moves toward more integrated approaches that couple deep Earth and surface processes, our pattern analysis of dynamic topography provides a useful means to evaluate the dynamic relationships between basin infill and sediment-sourcing regions to reconstruct regional, spatiotemporal vertical motion histories and support resource exploration.
Competing interests.The authors declare that they have no conflict of interest.
Acknowledgements.This research was supported by resources provided by the Pawsey Supercomputing Centre with funding from the Australian Government and the Government of Western Australia and with the assistance of resources from the National Computational Infrastructure (NCI), which is supported by the Australian Government.Sascha Brune was funded by the Marie Curie International Outgoing Fellowship 326115 and the Helmholtz Young Investigators Group CRYSTALS.Christian Heine was supported by ARC Linkage Project LP0989312 with Shell E & P and TOTAL.D. Rhodri Davies is funded by an ARC Future Fellowship (FT140101262) and Simon Williams and R. Dietmar Müller are supported by ARC grants DP130101946 and IH130200012.Leonardo Quevedo is acknowledged for the numerical routines to compute subduction volumes and the implementation of the GPlates TERRA output routines, along with John Cannon.The authors thank the employees of both supercomputing centres for their generous support and Geoscience Australia for their vast, open and easily accessible database.Some figures were generated using the Generic Mapping Tools (GMT; Wessel et al., 2013).
Edited by: J. Huw Davies Reviewed by: two anonymous referees

Figure 1 .
Figure 1.Schematic depiction of model run showing the evolution of the model from 200 Ma to present day.The velocity field as denoted by the coloured arrows is generated by the plate tectonic reconstruction software GPlates based on the Earthbyte 2013 plate model (see text).This velocity field is used as a time-dependent velocity boundary condition for the TERRA model runs.The model temperature field within the lower mantle is shown in terms of the 2000 and the 2800 K isotherm depicted as blue and red isosurfaces, respectively.

Figure 3 .
Figure 3. Dynamic topography time series of selected locations for all three models (M1-M3).Vertical motions are dominated by the plate tectonic history, while variations in core-mantle boundary temperature (M2) and lower mantle viscosity (M3) have only marginal effects on dynamic topography.For locations of sample points, see Figs. 11 and 12.

Figure 4 .
Figure 4. Evolution of global dynamic topographic elevation from 200 Ma to present.Blue colours denote negative dynamic topography, and red colours mark positive dynamic topography.(a) Dynamic topography in the absolute reference frame.The indicated points show synthetic sample locations as discussed in the text.Note that sample points move along with their associated tectonic plates.The evolution of dynamic topography in 10 Myr intervals can be found in the Supplement.(b) Dynamic topography in the plate frame of reference.The evolution of dynamic topography in 10 Myr intervals can be found in the Supplement.

Figure 6 .
Figure 6.Dynamic topography evolution for Australia, analogous to Fig. 4a.Between 200 and 120 Ma, Australia's eastern margin experiences westward-directed subduction of the Phoenix plate.From 50 Ma onward Australia's dynamic topography is dominated by the northward-directed subduction of the Australian plate under the Indonesian archipelago.See the Supplement for dynamic topography evolution in 10 Myr time steps.

Figure 7 .
Figure 7. Tectonic subsidence curves for selected wells in Australia (red) compared to our model output (coloured).Well data curves (accordingly coloured) are offset to model elevation at present day.

Figure 8 .
Figure 8.Southern African (sample point southern Africa east) uplift resulting from our models.While it successfully reproduces the Cretaceous uplift as described in previous work, uplift in the Miocene period remains comparatively small (see text).

Figure 9 .
Figure 9. Overview map (top) and time series (bottom) of dynamic topography for Category I, "topographic stable areas".

Figure 11 .
Figure 11.(a) Overview (top) and time series (bottom) plot of dynamic topography for Category II, "dynamically subsiding areas".(b) Overview (top) and time series (bottom) plot of dynamic topography for Category II, "dynamically subsiding areas", continued.

Figure 12 .
Figure 12.Overview (top) and time series (bottom) plot of dynamic topography for Category III, "fluctuating areas".

Figure 13 .
Figure 13.Comparison of dynamic topography time series at three sample locations with the relative orthogonal velocity of the down-going and overriding plate at the nearest point on the subduction zone.For locations of sample points, see Fig. 12.

Figure 14 .
Figure 14.The cluster analysis underpins the suggested categorisation; k = 5.(a) The dark blue cluster consists of locations in stable areas, light blue corresponds to (long-term) dynamically subsiding areas and the fluctuating category is encompassed by three clusters (green, orange and red).(b) Different clusters of elevation time series.The average trend is marked in black.See the Supplement for alternative choices of k.

Table 1 .
Physical parameters common to all model runs.See Table2for the parameters varied between models.