A review of observations and models of dynamic topography

Knowledge of the effect of mantle flow on surface topography has considerably increased over the past 30 yr. The rapid improvement in computational algorithms and computing resources has facilitated the modeling of global mantle flow at increasing resolution, which now achieves Earth-like convective vigor. The dramatic expansion of global and regional seismic data sets and numerical methods has led to commensurate improvements in the resolution and global coverage of mantle tomographic images (e.g., Romanowicz, 2008). In turn, these have led to an improved understanding of the thermo-chemical structure of the mantle (Grand, 2002; Masters et al., 2000; Simmons et al., 2007). Independently, global plate reconstructions that are used to constrain the evolution of mantle flow and to interpret seismic images of the mantle have likewise improved, and such reconstructions are now available back to the Triassic in million-year increments (Seton et al., 2012). As a consequence, plate kinematics and seismic tomography can now be assimilated into both forward and backward mantle flow models. Exploiting such improvements, various known long-wavelength vertical motions of continental plates have been attributed to mantle flow and hence used as constraints on the spatial character, amplitude, and time dependence of mantle flow. Examples include the Cenozoic uplift of southern Africa (Gurnis et al., 2000), the late Cenozoic uplift of the Colorado Plateau (Moucha et al., 2009), the Cretaceous subsidence and subsequent uplift of the interior of North America (Liu et al., 2008; Mitrovica et al., 1989), the tilt of northern South America to the east during the Miocene (Shephard et al., 2010), the tilt of Australia since the Late Cretaceous (DiCaprio et al., 2009), and the vertical motions of the Slave and Kaapvaal Cratons since the Paleozoic (Zhang et al., 2012). An important concept that derives from such studies is that there is likely “no such thing as a stable continental platform” (Moucha et al., 2008a), which in turn implies that long-term global sea-level (eustasy) change cannot be defined based on the detailed analysis of the stratigraphic record in a single area, as has been argued previously (e.g., Miller et al., 2005).

Nevertheless, multiple challenges still limit our understanding of global dynamic topography, including different definitions of the meaning of the term “dynamic topography,” inaccurate estimates of present-day dynamic topography in the absence of a detailed global model of the structure of the lithosphere, and contradicting model predictions of vertical motions for the same area (Moucha et al., 2008a; Müller et al., 2008b; Spasojevic et al., 2008).

In this review, we first attempt to clarify the definition of dynamic topography. We then explore observational constraints on present-day and time-dependent dynamic topography. For the present day, we compare three published residual topography fields (Kaban et al., 2003; Panasyuk and Hager, 2000; Steinberger, 2007) and a residual topography field that we derive, and we show that they are in broad agreement. We then carry out a non-exhaustive review of analytical and numerical models of dynamic topography, both instantaneous and time dependent, with an aim of not only showing the progressive improvements of these models, but also to expose the strengths and limits of each. As part of this process, we compare the present-day dynamic topography predicted by four published models (Conrad and Husson, 2009; Ricard et al., 1993; Spasojevic and Gurnis, 2012; Steinberger, 2007) and by the model that we introduce in this study. This comparison shows that the predicted dynamic topography fields are in good agreement at long wavelength, although the predicted amplitudes differ from one model to another. The residual topographies and modeled dynamic topographies are in broad agreement, although there are important regional mismatches. For example, East and Southeast Asia are highs in most residual topography fields, but they are also dynamic topography lows. We attribute this discrepancy to the difficulty in constraining residual topography with limited extents of dated oceanic lithosphere. In addition, the northeast Pacific present-day residual topography disagrees with modeled changes in dynamic topography, which shows the limitations in constraining a process that occurs over tens of millions of years using a snapshot of residual topography. We next discuss the consequences of dynamic topography for long-term sea-level change. Although the global effect of dynamic topography on long-term sea-level change, dominated by changes in the volume of ocean basins, is secondary (<100 m; Conrad and Husson, 2009; Spasojevic and Gurnis, 2012), its regional effects are important. To illustrate this, we carry out a detailed analysis of the anomalous subsidence of COST-B2 well offshore New Jersey, which illustrates that this margin has subsided by up to 400 ± 45 m since 85 Ma. We show that this dynamic subsidence is consistent with the model presented in this study and with two previously published models (Müller et al., 2008b; Spasojevic et al., 2008). However, given the large uncertainty associated with the dynamic models, the detailed sea-level curve derived from the stratigraphy of the New Jersey margin cannot be corrected for dynamic effects, which implies that it should not be used as long-term global sea-level curve. From an observational point of view, detailed knowledge of the thermal and chemical structure of the lithosphere is a key factor in constraining present-day dynamic topography. From a modeling point of view, explicitly including the lithosphere in future global mantle flow models should not only facilitate the comparison between modeling results and the geological record, but also make dynamic topography an additional constraint on the nature of the coupling between the mantle and the lithosphere.

Dynamic topography is that topography due to flow within the mantle, an idea that was first proposed by Pekeris (1935). It is “dynamic” because the mass anomalies driving the density are moving, and this contrasts with isostatic topography, in which the mass anomalies are in a state of quasi-equilibrium. Due to plate tectonics, the topography of Earth consists of continents and oceans with an average elevation difference of ∼4.5 km. This first-order difference in elevation is of isostatic origin and results from the different average densities and thicknesses of continental and oceanic crust and lithosphere. The long-known poor correlation between the ocean-continent function and the long-wavelength geoid (Kaula, 1972) supports this view. If the topography of Earth were not affected by any dynamic processes, flat continents would stand ∼4.5 km above flat abyssal plains. This is clearly not the case, and the topography of Earth results from an array of dynamic processes operating at different scales in time and space. In the continents, tectonic convergence results in mountain belts such as the Himalayas and the Andes, and tectonic divergence results in rifts, such as the East African Rift or the Salton Sea in California. The process leading to these topographic features is clearly dynamic, but even in active and geologically recent tectonic zones, the topography is largely compensated isostatically: The continental crust is thicker than average in mountain belts, and thinner than average in rifts (Mooney et al., 1998).

In the oceans, the bathymetry is primarily determined by cooling and thickening of oceanic lithosphere as it moves away from spreading centers (Langseth et al., 1966; McKenzie, 1967). While it is ultimately dynamic in origin, the bathymetry of the young ocean floor (<80 Ma) is well explained by combining the principles of isostasy and of the cooling of a semi-infinite half space (Davis and Lister, 1974). Ocean floor older than ca. 80 Ma appears to flatten compared to this simple and elegant physical model, although whether “normal” old seafloor flattens is subject to some debate (Hillier, 2010; Korenaga and Korenaga, 2008). The reason for the apparent shallowing of old ocean floor is the subject of ongoing debate, and proposed explanations include thermal rejuvenation by hotspots (Smith and Sandwell, 1997), small-scale convection (e.g., Afonso et al., 2008), dynamic support by mantle flow (Kido and Seno, 1994; Zhang et al., 2012), and departure from the topography predicted by thermal boundary layer theory due to radioactive heat production within the mantle (Jarvis and Peltier, 1982). In the absence of consensus regarding the physical process explaining this flattening, the age-depth relationship of the ocean floor can still be accommodated in empirical “plate models” (e.g., Crosby and McKenzie, 2009; Parsons and Sclater, 1977). As in the continents, the topography of the ocean floor is affected by variations in crustal thickness, notably in areas affected by magmatism such as oceanic plateaus. Similarly, variations in lithospheric temperature and magmatic processes also affect continental topography; however, there is no simple tectonic model to describe the topography of continents due to their long and complex tectonic history.

Given the number of dynamic processes that affect topography, it is not surprising that the term “dynamic topography,” which generally describes the effect of mantle convection on surface topography, does not have a single, agreed-upon definition. As discussed already, the bathymetry of the ocean floor is compatible with that of the thermal boundary layer of the convecting mantle, and, in this respect, the subsidence of oceanic lithosphere with age should be included in the definition of dynamic topography. However, the dynamics of the continental lithosphere are more complex and cannot be explained by boundary layer theory. For this reason, defining dynamic topography globally as the difference between observed topography and isostatic crustal topography (e.g., Forte et al., 1993) results in continents dynamically depressed by up to 3 km, which is inconsistent with observations (Gurnis, 1993a). This discrepancy is in part due to the non-negligible isostatic contribution of the continental lithosphere (in the past referred to as tectosphere; Jordan, 1975) to topography. Given the relatively poor global constraints on the complex structure and density of the continental lithosphere (Zoback and Mooney, 2003), dynamic topography is usually defined as that topography originating from sources beneath the upper thermal boundary layer of mantle convection (e.g., Le Stunff and Ricard, 1995; Panasyuk and Hager, 2000). We use this latter definition in the present paper, keeping in mind that it ignores the many lithospheric-scale dynamic processes that affect topography, including those described earlier herein.

Constraining dynamic topography requires removing the isostatic contribution of the sediments, ice, crust, and lithosphere from the observed topography. The remaining topography is called residual topography (Crough, 1983). In the oceans, a choice must be made between correcting either for the physical half-space cooling model (as in Davies and Pribac, 1993; Gurnis et al., 2000) or for a plate model accounting for the apparent flattening of the ocean floor with increasing age (Kaban et al., 2003; Le Stunff and Ricard, 1995; Steinberger, 2007). The former approach implies that mantle support may contribute to the observed flattening of old ocean floor (Kido and Seno, 1994), whereas the latter requires some other physical process to operate. In the continents, a common approach is to remove the isostatic topography due to the continental crust (Forte et al., 1993; Steinberger, 2007). However, correcting for the isostatic contribution of the continental crust alone can result in continents being either systematically lowered or elevated by 1 to 2 km, which is inconsistent with observation (Gurnis, 1993a). Zoback and Mooney (2003) showed that there is no simple relationship between continental crustal thickness and surface elevation, and that isostatic models based on crustal thickness alone systematically overpredict continental elevations, mainly because Archean continental lithosphere is depleted and of intrinsically lower density (Jordan, 1975). Several authors have proposed residual topography maps that remove the isostatic contribution of both the continental crust and continental lithosphere (Kaban et al., 2003; Le Stunff and Ricard, 1995; Panasyuk and Hager, 2000). Such calculations require assumptions relative to the thermo-chemical state of the continental lithosphere, which result in relatively large uncertainties (Kaban et al., 2003; Le Stunff and Ricard, 1995; Panasyuk and Hager, 2000) and underline the need for a better knowledge of the thermal and chemical structure of the lithosphere in order to constrain mantle flow models.

The long-wavelength residual topography has been calculated by removing the mean continental elevation from continents (following Gurnis et al., 2000; Lithgow-Bertelloni and Silver, 1998), because the thermal and chemical structures of the continental lithosphere are not well known. In the oceans, the plate cooling model of Crosby and McKenzie (2009) was removed from the observed bathymetry after accounting for isostatic loading by sediments using the method of Sykes (1996) and the global sediment thickness map of Divins (2003). We consider that continents extend to the edge of the continental shelf (magenta contours in Fig. 1A), which roughly corresponds to the −200 m elevation contour (Harrison et al., 1983), and we calculate the equal-area arithmetic mean continental elevation as 529 m. The areas between the edge of the continental shelf and ocean floor per se as defined by (Müller et al., 2008a) were arbitrarily set to a uniform elevation of −200 m. The resulting residual topography was low-pass filtered using a filter tapered for wavelength between 2000 and 5000 km, in order to simplify comparison with published estimates of residual topography expanded in spherical harmonics to degree 12 (Kaban et al., 2003; Panasyuk and Hager, 2000; Steinberger, 2007; Figs. 1B, 1C, and 1D). The residual topography contains an isostatic component evident in high values in the Himalayas and Andes (with present crustal thicknesses greater than 50 km [Mooney et al., 1998], shown as yellow contours in Fig. 1A). In the oceans, the residual topography highs are closely associated with hotspot tracks (Cazenave et al., 1989; Schroeder, 1984) and Phanerozoic large igneous provinces (denoted by green contours in Fig. 1A). Therefore, part of the oceanic residual topography is due to the thicker crust of oceanic plateaus and at least partially isostatically compensated. Ideally, a correction for the thickness of the oceanic crust should be carried out to estimate the topography associated with mantle support. Such studies have been carried out regionally (e.g., Crosby and McKenzie, 2009; Winterbourne et al., 2009) but not yet globally. Finally, the residual topography high over the Canada Basin in the Arctic Ocean (Fig. 1A) can be attributed to the absence of sediment thickness data.

Despite being based on different isostatic models of the lithosphere, the four residual topography fields (Fig. 1) agree overall. For half-wavelengths greater than ∼3300 km, the peak-to-peak amplitude is ±2 km, and the root-mean-square amplitude is ∼0.5 km (Table 1). The peak-to-peak amplitude is smallest for the model of Panasyuk and Hager (2000) and largest for the model presented in this study (Table 1), due to the absence of correcting for isostasy in elevated areas (Fig. 1A). The root-mean-square amplitude is smallest for the model of Panasyuk and Hager (2000) and largest for the model of Kaban et al. (2003) (Table 1). In terms of pattern, in the continents, features common to all four models include the residual topography lows within central Europe, eastern North America, and western South America, and the residual topography highs within southern and eastern Africa and western North America. In the oceans, the most persistent features between models are the Australian-Antarctic Discordance (AAD) that depresses the southeast Indian mid-ocean ridge, the Argentine Basin in the southern South Atlantic, the residual topography low under the western part of the central Atlantic Ocean, the large-scale, low-amplitude residual topography low under the northeast Pacific Ocean, and the highs under the central and western Pacific Ocean (Darwin Rise), the areas offshore southern Africa, and around Iceland and the northern Atlantic mid-ocean ridge.

In principle, dynamic topography can occur over all wavelengths from 10 km to greater than 1000 km, but since crustal and lithospheric density and thickness variations and plate flexure exert a strong control on topography for wavelengths less than 100 km, most attempts to deconvolve it from the geological record have focused on transient long-wavelength changes that cannot be well explained without invoking dynamic topography. In this section, we first focus on the geological expression of long-wavelength dynamic topography, which is the main focus of this paper, before giving a brief overview of short-wavelength dynamic topography.

Long-wavelength dynamic topography is recorded by lithospheric plates moving over the constantly changing mantle flow and is thus ephemeral over millions of years to tens of millions of years. As a result, subsidence events are generally followed by uplift and vice versa. Indeed, the fact that dynamic topography is transient over tens of millions of years means that it must be recorded over long time periods, over which the preservation potential of the rock record is poor. The very nature of the long-wavelength, low-amplitude dynamic topography signal makes it difficult to isolate in the geological record, which is dominated by shorter-wavelength, larger-amplitude signals such as lithospheric extension or shortening. One further restriction on the preservation of dynamic topography comes from the geological record itself, which is biased toward areas that have been subsiding for long periods of time and have a good preservation potential, as opposed to uplifting areas that are subject to erosion and have a poor preservation potential. However, strata deposited in a dynamic topography depression also have low long-term preservation potential because dynamic topography is transient. Thus, ancient subduction-related dynamic topography is most likely to be represented by unconformities (Burgess et al., 1997). As a consequence, the stratigraphic record is the most accessible archive of the dynamic topography signal, and continental interiors that have not undergone deformation for hundreds of millions of years are the most promising places to search for the dynamic topography signal (Gurnis, 1992).

The Western Interior Seaway of North America, which covered ∼40% of the continent and resulted in >1 km of sediments deposited over ∼1000 km (Bond, 1976; Liu and Nummedal, 2004), presents a robust long-wavelength signal. Analysis of Cretaceous sediments in North America suggests that the observed flooding would require a sea-level rise of 310 m, resulting in accumulation of 700 m of sediments (Bond, 1976), but the inferred Cretaceous isopachs (Cook and Bally, 1975) are significantly thicker. Late Cretaceous subsidence was followed by Tertiary uplift with a maximum tilt amplitude of 3 km over ∼1400 km (Mitrovica et al., 1989). During the Late Cretaceous, when the sedimentary wedge was thickest, the depositional system was characterized by an eastward migration of a depocenter from Utah to Wyoming (Weimer, 1970; Fig. 2A). The dominant characteristic was a westward downward tilt with an internal depocenter that migrated eastward at 5 cm/yr within North America (Liu et al., 2011; Figs. 2B–2F). Besides this dominating Late Cretaceous to Tertiary signal, Heller et al. (2003) identified three thin but widespread alluvial conglomeratic units that recorded post-Paleozoic episodes of widespread tilting across the U.S. Cordillera that have been interpreted in terms of dynamic topography. Finally, a combination of stratigraphic, thermochronologic, and provenance constraints showed that the southwestern corner of the United States tilted to the east at ca. 80 Ma and then to the west at ca. 55 Ma over scales of ∼500 km with an amplitude of ∼1 km (Wernicke, 2011).

Episodes of anomalous long-wavelength motions have also been inferred for Australia, which has experienced little orogenic activity since the mid-Mesozoic and is relatively flat. During the Cretaceous, the marine inundation of Australia and eustasy were out of phase. Maximum flooding of Australia occurred in the late Aptian to early Albian (120–110 Ma) (Struckmeyer and Brown, 1990), when large regions of the continent experienced marine inundation. Australia became progressively exposed during the Late Cretaceous, reaching a flooding minimum in the Campanian (80–70 Ma), when global sea level was near a maximum. This record of anomalous relative sea level suggests that Australia has experienced several episodes of vertical motion, including a large-scale vertical translation of the entire continent and a localized subsidence event (Bond, 1978; Russell and Gurnis, 1994; Veevers and Conaghan, 1984). Late Cretaceous vertical motions appear to have been mostly confined to the Eromanga and Surat Basins (Gallagher and Lambeck, 1989). The Eromanga Basin recorded ∼500 m of tectonic subsidence in 10 m.y. at around 100 Ma (Gallagher and Lambeck, 1989). Apatite fission-track analysis suggests that the Cretaceous strata of the Surat Basin and basins on the eastern margin were subsequently eroded (Gallagher et al., 1994), which is consistent with the tilting up toward the east of Cretaceous marine strata (Gallagher et al., 1994; Russell and Gurnis, 1994). Near the end of the Cretaceous, Australia was ∼250 m higher than it is today (Bond, 1978; Russell and Gurnis, 1994; Veevers and Conaghan, 1984). During the Cenozoic, Australia moved northward after it separated from Antarctica and tilted down to the northeast (DiCaprio et al., 2009; Sandiford, 2007). This tilt must have had an amplitude of 300 m since the Eocene in order to reconcile interpreted patterns of marine incursion with a predicted topography accounting for global sea-level variations (DiCaprio et al., 2009). On the southern margin, long-wavelength dynamic topography was enhanced by at least 250 m of shorter-wavelength anomalous subsidence, consistent with the passage of the margin over a north-south–elongated, 500-km-wide, anomaly of dynamic topography approximately fixed with respect to the mantle (DiCaprio et al., 2009). The present-day position of this depth anomaly is aligned with the Australian-Antarctic Discordance (cf. Whittaker et al., 2010).

In addition to continental interiors, some methods have been developed to estimate dynamic topography in areas that have undergone recent tectonic deformation. For instance, postrift anomalous subsidence, the difference between observed postrift thermal subsidence and that predicted by an analytic model of lithospheric cooling, has been attributed to dynamic topography on the Queensland Plateau on the northeastern margin of Australia (Müller et al., 2000), in the South China Sea (Xie et al., 2006), between New Zealand and Antarctica (Sutherland et al., 2010), and at the New Jersey margin (see Discussion).

The previous section shows that the stratigraphic record is a powerful constraint on dynamic topography in areas that have experienced long-term subsidence, especially in areas where long-term subsidence has been followed by uplift. Different methods are required to estimate dynamic topography in areas experiencing long-term dynamic uplift, with southern Africa being a well-known example. Both the timing and amplitude of uplift are generally difficult to constrain (Gurnis et al., 2000), due to the poor preservation potential in elevated areas and to the large uncertainty of paleo-altimetry proxies, i.e., typically a few hundred meters. As a result, dynamic uplift may be constrained indirectly by estimating the timing and amplitude of rock uplift. In this respect, thermochronologic methods allow the estimation of the amount of exhumation by constraining both the amplitude and timing of the cooling of uplifted rocks. For instance, indirect estimates of vertical motions via thermochronology have been made for Eromanga Basin in Australia (Gallagher et al., 1994), the Slave Craton (Ault et al., 2009), and Kaapvaal Craton (Flowers and Schoene, 2010) using apatite (U-Th)/He (AHe) or apatite fission-track (AFT) thermochronology. One limit of thermochronologic methods, however, is the nonuniqueness of thermochronologic modeling (Redfield, 2010), which underlines the need to complement it with other, independent, observations (Green et al., 2011).

A recent approach to constraining dynamic uplift consists in “inverting” topographic river profiles by parameterizing erosion in both time and space (Pritchard et al., 2009), which is somewhat different from conventional geomorphology, in which uplift is parameterized to estimate erosion. This method has been used to estimate the uplift of the Bié dome, Angola, over the past 40 m.y. (Pritchard et al., 2009). In this case, the method constrains the uplift history of single watersheds over ≤1000 km and ≤30 m.y., and is thus complementary to the stratigraphic and thermochronological records, which provide information over longer time scales and larger spatial extents. An alternative method is that of Guillaume et al. (2009), who estimated the tilt of Miocene terraces to estimate the Neogene uplift of central and eastern Patagonia, possibly due to the opening of a slab window. The main difficulty of this method is the dating of terraces.

Mantle-driven dynamic processes may also affect topography at shorter wavelengths (100 km to a few hundred kilometers). A major challenge when studying dynamic topography at such wavelengths is to isolate the signal from the contributions of isostasy and elastic flexure, which have similar, kilometer-scale, amplitudes. The most vivid expression of dynamic topography is found in oceanic trenches that are typically <200 km wide yet several kilometers deep (Husson et al., 2012; Zhong and Gurnis, 1994).

Delamination and convective instability of the continental lithosphere are short-wavelength dynamic processes with the potential for a strong topographic expression. A compelling empirical case is that of the renewed uplift of the southern Sierra Nevada Mountains since the mid-Pliocene (Clark et al., 2005) and localized subsidence within the adjacent Great Valley (Saleeby et al., 2012), which have been related to instability of the continental lithosphere evident in a prominent high-seismic-velocity anomaly that extends from the base of the crust to ∼200 km depth (Boyd et al., 2004; Zandt et al., 2004). Other examples of vertical motions associated with continental lithospheric instability include Quaternary uplift in Baja, California (Mueller et al., 2009), late Cenozoic uplift of the Wallowa Mountains in eastern Oregon, U.S.A. (Hales et al., 2005), and uplift of Miocene marine sediments in the Central Anatolian Plateau, Turkey (Cosentino et al., 2012).

Smaller-scale modes of convection may influence surface topography. A prominent mode evident in models of mantle convection is small, low-amplitude instabilities that travel along the base of the cold, top thermal boundary layer of mantle convection (e.g., Huang et al., 2003). Such traveling convective instabilities might be responsible for vertical motions on a 2–20 m.y. time scale, as evident in sequence stratigraphic analysis of passive continental margins (Petersen et al., 2010). This raises the possibility that such localized small-scale convection could be responsible for high-frequency, second- to third-order sea-level fluctuations. Whether the second- and third-order sea-level changes evident on passive margins are global or local is controversial (e.g., Hallam, 1992; Lovell, 2010). The upwarped edges of the Colorado Plateau in the western United States, as well as the localization of late Neogene–Quaternary magmatism and steep upper-mantle seismic velocity gradients, are potentially indicative of this mode of sublithospheric, small-scale convection (van Wijk et al., 2010).

Although the uplift and high elevation of the Colorado Plateau (e.g., Flowers, 2010) have been interpreted as part of the overall long-wavelength Cenozoic uplift of the western United States (see previous), they might be viewed as a shorter-wavelength component of dynamic topography. The vertical motions evident from uplift of Mesozoic marine strata, unroofing from low-temperature thermochronology (Flowers et al., 2008), and reorganizations of drainage patterns (Wernicke, 2011) might be associated with a long-wavelength signal associated with the removal of the negatively buoyant Farallon slab (Liu and Gurnis, 2010) or with the upwelling or warm mantle that extends from the top of the lower mantle into the upper mantle (Moucha et al., 2008b, 2009), or with a combination of both processes.

Locations of active hotspot volcanism correspond to swells of residual topography that have a wavelength of several hundred kilometers and an elevation <2 km (Crosby and McKenzie, 2009; Crough, 1983). The surface expression of mantle plumes likely extends beyond this main swell. For instance, lateral magmatic underplating and pulses of hot material in the asthenosphere associated with the Iceland plume have been proposed to explain transgression and regression patterns over <10 m.y. observed in the North Atlantic (Lovell, 2010). A Paleocene–Eocene landscape off the northwest coast of Europe was recently reconstructed using a three-dimensional seismic data set on the continental shelf north of Scotland (Hartley et al., 2011; Shaw-Champion et al., 2008). Data suggest that this landscape was lifted above sea level in a series of three discrete steps of 200–400 m each, and that it was reburied after ∼1 m.y. of subaerial exposure, which points to transient uplift and subsidence. These transient vertical motions have been interpreted as a pulse of warm asthenosphere that moves in a thin low-viscosity channel radially away from the Iceland plume (Rudge et al., 2008).

Several approaches have been developed to model regional and global dynamic topography, including models that only consider the present day (instantaneous models), if they consider the time domain at all, time-dependent models, and either a forward or inverse matching of present-day mantle structure. In addition, there has been a variety of studies that have looked at generic issues associated with the physics of dynamic topography and the connection with the rock record (Burgess and Gurnis, 1995; Mitrovica and Jarvis, 1985). Models either use or predict present-day quantities (such as the geoid, mantle seismic structure, surface and/or core-mantle boundary dynamic topography, the rate of change of surface topography, and plate motions) to infer mantle properties. Inverse models and models that use data assimilation can be compared to fewer observations since assimilated data cannot be used to validate the models. This section constitutes an overview of a number of regional- and global-scale models of dynamic topography, including their strengths and limitations (Table 2) while also looking at some generic issues associated with the ways in which physical properties influence dynamic topography, time dependence of dynamic topography, and its connection to the stratigraphic record.

Hager (1984) used the approach given by Equation 2 with a thermal model of slabs consistent with Benioff zone seismicity and the observed geoid to constrain a minimum 30-fold increase of viscosity from the upper to lower mantle. The approach was then further elaborated using a mantle density structure scaled from global seismic tomography inversions (Hager et al., 1985), with results that implied that chemically stratified mantle convection could be ruled out because it produced a pattern of dynamic topography with broad subsidence over Africa. To fit the observed geoid, Hager et al. (1985) revised the minimum viscosity increase from upper to lower mantle to a factor of 10. In contrast, Lithgow-Bertelloni and Silver (1998) fit the residual topography of southern Africa with the dynamic topography predicted by an instantaneous mantle flow model with a lower mantle 50 times more viscous than the upper mantle. These differences in best-fit mantle viscosity structures are not surprising considering that the two flow models are based on different tomography models and that Hager et al. (1985) included all seismic velocity anomalies as density anomalies in their calculation of surface dynamic topography, whereas Lithgow-Bertelloni and Silver (1998) calculated the surface dynamic topography resulting from density heterogeneities below 325 km depth.

Steinberger (2007) and Conrad and Husson (2009) derived the global dynamic topography predicted by instantaneous numerical mantle flow models (Figs. 4B and 4D). These two models, hereafter referred to as S07 and CH09, are broadly similar in that they are based on a mantle density structure derived from S-wave tomography with the assumption that all lateral variations in seismic velocities are thermal, while ignoring density contrasts in the uppermost mantle (220 km in S07; 300 km in CH09). Important differences between these two models include the tomographic inversion, the scaling factor to convert seismic velocity to density anomaly, the surface velocity boundary conditions, and the viscosity structure. In addition, S07 is a semi-analytical model using a spherical harmonics approach (Hager and O’Connell, 1979) with a horizontal resolution of ∼650 km, whereas CH09 is a numerical model using a finite-element approach, CitcomS (Zhong et al., 2000), with a horizontal resolution of ∼100 km. These differences illustrate the large numbers of free parameters and boundary conditions inherent to mantle flow models.

Despite the major differences between these two models, the present-day dynamic topography fields they predict are generally similar (Figs. 4B and 4D) and compatible with proposed residual topography fields (Fig. 1). Both predicted dynamic topographies show highs within the Pacific, and southern and eastern Africa. The main dynamic topography lows extend under Central and South America, and from Australasia to eastern Europe, with a minimum under Southeast Asia. The main differences in terms of patterns are the dynamic topography lows in the central South Atlantic, in northwest Africa, and in the eastern United States predicted by S07 but not by CH09. All three of these dynamic topography lows are observed in residual topography fields (Fig. 1), although the amplitudes are substantially larger in S07 compared to CH09 and residual fields (Fig. 1). In addition, neither S07 nor CH09 reproduces the residual topography low in the northeast Pacific Ocean.

An important difference between S07 and CH09 is the amplitude of the calculated dynamic topography (Table 1), which has been the subject of a long-standing debate (e.g., Gurnis, 1990). The dynamic topography predicted by model S07 (CH09) presents a peak-to-peak amplitude of ±3 km (±1.5 km) and a root-mean-square amplitude of ∼0.9 km (∼0.5 km; Table 1). For comparison, the peak-to-peak amplitude of the residual topography fields shown in Figure 1 is approximately ±2 km, and their root-mean-square amplitude is ∼0.5 km (Table 1). Therefore, the dynamic topography predicted by model CH09 (S07) is 25% smaller (50% larger) than that of residual topography. The main parameters explaining the amplitude difference between the models are the scaling factor from seismic to density anomaly, which is greater in model S07, the depth above which density contrasts are ignored, which is 80 km shallower in S07, and the low-viscosity asthenosphere present in CH09 and not in S07. Steinberger (2007) showed that including the effect of phase transitions in such models only reduces the root-mean-square amplitude of dynamic topography by ∼5%. Other important effects to investigate include accounting for chemical heterogeneity when converting seismic velocity anomalies to density anomalies (Simmons et al., 2007) and the effect of lateral viscosity variations (Zhong and Davies, 1999). At long wavelengths, lateral variations in viscosity might be small compared to the uncertainty of tomography models, but they are expected to be significant in time-dependent mantle flow models that are largely driven by the evolution of thermal heterogeneity (Moucha et al., 2007).

The previous section shows that present-day dynamic topography can be estimated using the present-day density structure of the mantle as constrained by tomography models. However, a major challenge in estimating the evolution of dynamic topography is that the past density structure of the mantle is not directly known. Several strategies using data assimilation have been developed over the past two decades to model the evolution of mantle flow and of dynamic topography. These strategies include forward, backward, and adjoint mantle flow models based on the history of subduction and plate motions given by global tectonic reconstructions, or on mantle tomography, or on a combination of both.

Analytical models.Present-day dynamic topography. Ricard et al. (1993) developed a time-dependent analytical model of mantle flow in which subduction history (Lithgow-Bertelloni and Richards, 1998; Lithgow-Bertelloni et al., 1993) was imposed over the past 200 m.y. In this model, slabs are treated as “stokeslets,” sinking vertically at a constant rate and slowing down as they enter the lower mantle, while neglecting thermal diffusion. The best-fit model of Ricard et al. (1993) to the observed geoid was obtained for a lower mantle 40 times more viscous than the upper mantle. This estimate, based on a forward model, is closer to the analytical result of Hager (1984; lower mantle 30 times more viscous than upper mantle) than the best-fit model to the geoid of Hager et al. (1985), which showed a lower mantle only 10 times more viscous than the upper mantle. Again, such variations in results are not surprising given the differences among the three approaches. The present-day dynamic topography predicted by the model of Ricard et al. (1993) is shown in Figure 4C. It has a minimum of −1.6 km, a maximum of ∼600 m, and a root-mean square amplitude of ∼430 m (Table 1). The skew of this dynamic topography field toward negative values is due to the model being driven by slabs, while the return mantle flow is entirely passive: Negative dynamic topography is of larger amplitude and covers a smaller area than positive dynamic topography (Fig. 4C). As a consequence, the amplitude of negative dynamic topography predicted by this model is ∼20% smaller than negative residual topography (Fig. 1), whereas positive dynamic topography is only 30% of positive residual topography. In terms of spatial distribution, the dynamic topography lows predicted by the forward model of Ricard et al. (1993) under Central/South America and Southeast Asia match those predicted by the inverse models of Steinberger (2007) and Conrad and Husson (2009). However, the dynamic topography low under the Gulf of Alaska and Bering Sea predicted by Ricard et al. (1993) is not predicted by inverse models, and it is only observed in one of the four residual topography fields (Fig. 1), which suggests an overestimation of the amount of subducted material in that area in the plate model used by Ricard et al. (1993).

Evolution of dynamic topography.Gurnis (1993a) developed time-dependent models of mantle flow based on the analytical formulation proposed by Richards and Hager (1984) using the Phanerozoic subduction history of Scotese and Golonka (1992) to investigate the influence of dynamic topography and its time derivative on continental flooding. This study revealed that a significant part of the observed Phanerozoic continental flooding could be attributed to dynamic subsidence above active subduction zones, suggesting that dynamic topography affects global sea-level change at long wavelengths, and therefore that Phanerozoic stratigraphy can be used to constrain mantle flow. Building on this approach, Lithgow-Bertelloni and Gurnis (1997) used a series of five instantaneous global flow models throughout the Cenozoic to reproduce the overall observed trend of flooding and exposure of North America, Australia, and Indonesia, although the amplitudes of model uplift and subsidence were greater than that inferred from stratigraphic constraints.

Forward numerical models.Gurnis et al. (1998) developed numerical regional mantle convection models with imposed plate kinematics to investigate the evolution of the dynamic topography of the Australian plate since the Cretaceous. Gurnis et al. (1998) used the location of subduction taken from plate reconstructions (Müller et al., 1997) to build a synthetic initial temperature field for their forward models. They reproduced the regional Cretaceous flooding of Australia deduced from paleogeographic reconstructions in which it had long been known that eastern Australia achieved peak flooding with shallow marine seas ∼20 m.y. earlier than inferred peak global sea level (see previous; Veevers and Conaghan, 1984). Essentially, in the Early Cretaceous Australia drifted eastward over the long-lived subduction zone bounding Gondwanaland, which caused widespread marine inundation in the central and eastern part of the continent. Sediments filled the accommodation space as the continent slowly rose following termination of subduction such that when global sea level peaked in the Late Cretaceous, there was no remaining accommodation space, resulting in an absence of any major inundation. Gurnis et al. (1998) also showed that the position and history of the equally anomalous Australian-Antarctic Discordance (AAD; Fig. 1), a prominent present-day north-south residual topography low, could also be fit with the same model. The study is notable for two reasons. First, this rather straightforward computational approach linking plate tectonics and thermal convection showed that an entire cycle of dynamic topography could be followed from subsidence to uplift. As we shall see, most other regional applications can only follow part of the entire history (either the subsidence or the uplift phase). This is not a criticism of modeling methods but reflects fundamental limitations of the geological record. In addition, the study reinforced the need to follow the history of vertical motions in the plate frame of reference, that is, tracking the history of a point on a plate that moves with respect to the mantle. This is essential for the detailed comparison of models to geological proxies of vertical motions.

DiCaprio et al. (2010) tracked the evolution of dynamic topography for the Australian plate in a kinematically driven forward model with assimilation of plate velocities in 1 m.y. increments in which a high-resolution regional model was embedded in a global model. This model, starting at 50 Ma with subduction zones inferred from plate reconstructions, reproduced the observed fast Neogene subsidence of the Australian northeast shelf that would have contributed to the demise of carbonate reefs (DiCaprio et al., 2010). DiCaprio et al. (2010) were also able to exploit new paleogeographic tools that allowed the plate boundaries to smoothly evolve (Gurnis et al., 2012), essential for time-dependent convection models. Zhang et al. (2012) proposed a kinematically driven forward model extending back to 450 Ma and having an average global surface resolution of ∼100 km. The plate reconstruction used in this model consists of 34 stages over the past 450 Ma, and the initial condition was obtained by running the first reconstruction stage for 150 m.y. The vertical motions predicted by this model for the Kaapvaal and Slave Cratons were in agreement to the first-order with the burial histories predicted by thermochronology models.

A new model. We introduce a forward model in which we assimilate the thermal and kinematic field associated with plate reconstructions (Seton et al., 2012) with continuously closing plates (Gurnis et al., 2012) prepared using the software GPlates (Boyden et al., 2011). This model is essentially the extension of the approach of Ricard et al. (1993), Lithgow-Bertelloni et al. (1993), and Lithgow-Bertelloni and Richards (1998) to “a full three-dimensional convection problem constrained by plate motions” (Ricard et al., 1993, p. 21,896). The velocity field of plate reconstructions is continuously imposed as a surface boundary condition by linear interpolation between reconstruction stages defined every million years. The temperature field of the oceanic lithosphere is assimilated using paleo-oceanic ages (Seton et al., 2012) and the half-space cooling model (e.g., Davis and Lister, 1974). One limitation of global mantle convection models with imposed plate velocities is that they result in substantial advective thickening of slabs instead of well-defined one-sided subduction (Table 2). To overcome this limitation, we assimilate the temperature field in the proximity of the subduction zones given by the plate reconstructions. The thermal structure of slabs is derived from the age of the subducting lithosphere using a half-space cooling formulation and assuming a slab dip angle of 45°. This analytical solution is assimilated in the dynamic model within 350 km of subduction zones using a blending function to smoothly merge the analytical and numerical temperature fields. In addition, we ensure that subduction zones appearing during the modeled period are progressively inserted into the upper mantle. The method ensures that the slab buoyancy flux in the upper mantle is consistent with the plate reconstructions, and it avoids unrealistic slab advective thickening. The model, solved using CitcomS (Zhong et al., 2000), consists of 12 × 65 × 128² elements with an average surface resolution of ∼50 km. The initial condition (at 200 Ma) consists of a global temperature field (nondimensional background of 0.5) with cold slabs inserted to a depth of 1750 km and an assimilated cold surface thermal boundary layer as described previously (for a cross section of the initial condition under the Tethys, see Zahirovic et al., 2012). The temperature dependence of viscosity is implemented using an Arrhenius law with activation energy of 100 kJ/mol (33 kJ/mol) in the upper (lower) mantle, which limits viscosity variations to three orders of magnitude, and in our reference case, the lower mantle reference viscosity is 100 times larger than that of the upper mantle (for the resulting viscosity profile, see Zahirovic et al., 2012). The surface dynamic topography is computed by restarting the calculation every 5 m.y. with no-slip surface boundary conditions to suppress the effect of the surface traction imposed by plate velocities (for a discussion of the effect of model boundary conditions on predicted dynamic topography, see Thoraval and Richards, 1997) and after removing the uppermost 350 km of the mantle, i.e., the depth to which subduction zones are assimilated.

The resulting present-day dynamic topography field (Fig. 4A), broadly similar to that of Ricard et al. (1993; Fig. 4D), has a degree two pattern with a belt of low dynamic topography extending from Antarctica through the Americas to the Arctic and broadening beneath the whole Eurasian continent. In this subduction-driven model, the dynamic topography highs correspond to passive upwelling under the Pacific Ocean, southern Africa, and the Atlantic Ocean. This dynamic topography field is in overall agreement with the residual topography in terms of pattern (Fig. 1) and amplitude (between −2.2 km and +1.6 km; Table 1). Although the mean dynamic topography is zero, the modeled highs have lower amplitude than lows because upwelling is passive. The model features dynamic topography lows under central Europe, eastern North America, and western South America, and a dynamic topography high under southern Africa, which are all evident in the residual topography (Fig. 1). The dynamic topography under the Argentine Basin and the western United States has the expected sign, but its amplitude is lower than expected. The AAD is also of lower amplitude than residual topography, and it is shifted to the west. This discrepancy is tempered by the fact that only about half the excess depth of the AAD is likely dynamic, the rest being due to thin oceanic crust within the AAD (Gurnis et al., 1998; Whittaker et al., 2010). In addition, the positive dynamic topography predicted under the northeast Pacific Ocean is at odds with the large-scale, low-amplitude residual topography low (Fig. 1). Finally, the residual topography highs associated with plumes in eastern Africa, the western Pacific (Darwin Rise), and Iceland are not reproduced by the model, which only captures large-scale, passive mantle upwelling.

Backward advection. In order to investigate the time dependence of mantle flow and dynamic topography, a natural extension of instantaneous flow models consists of backward advecting the mantle density structure derived from tomography. This process essentially consists in changing the sign of gravity in the momentum equation and setting the nonreversible thermal diffusivity to zero in the energy equation (Steinberger and O’Connell, 1997). Some caution is needed, since the method accumulates light material at the core-mantle boundary and dense material toward the surface, which results in a decrease of information on the density structure of the mantle back in time. A further limit of backward advection is that thermal diffusion cannot be simply reversed (Ismail-Zadeh et al., 2004). Simple backward integration is adequate for long wavelengths where advection dominates over diffusion. The first study to model dynamic topography using this approach was that of Conrad and Gurnis (2003), who proposed that the dynamic topography high currently located under southern Africa migrated from eastern Africa over the past 36 Ma. They also pointed out the importance of thermal boundary layers in the backward advection process and used a half-space cooling model to assimilate oceanic ages in the top thermal boundary layer. This backward advection scheme, in which plate velocities are also assimilated, was reversible for ages up to between 50 and 75 Ma (Conrad and Gurnis, 2003).

Moucha et al. (2008a) estimated dynamic topography both at the surface, following the definition of Forte et al. (1993), and beneath the lithosphere, similar to the definition used here, from a backward advection model in which data were assimilated in several ways. First, the tomography model used by Moucha et al. (2008a) is a hybrid tomography/instantaneous mantle flow model that assimilates the free-air gravity anomaly, the “residual” topography field corrected for the isostasy of the crust, present-day plate velocities from DeMets et al. (1990), and the excess ellipticity of the core-mantle boundary (Simmons et al., 2007). Second, the viscosity profile is constrained using the same convection-related observables as done earlier, as well as data related to glacial isostatic adjustment (Mitrovica and Forte, 2004). Third, although plate motions are predicted by the flow rather than driving the flow in the formulation of Moucha et al. (2008a), mantle flow is backward advected using boundary conditions consistent with a no-net-rotation tectonic plate model. Aside from the data assimilation, the model of Moucha et al. (2008a) does not seem to include any specific treatment of thermal boundary layers, which is important in backward mantle advection schemes (Conrad and Gurnis, 2003), nor of lateral viscosity variations or temperature dependence of viscosity, which are important in time-dependent flow models (Moucha et al., 2007). Based on this model, Moucha et al. (2008a) proposed that the New Jersey margin has been dynamically uplifted by 100–200 m over the past 30 m.y., and its conjugate margin in northwest Africa has been uplifted by ∼100 m during the same period. The important consequence of this work is that sea-level curves derived from the well-documented stratigraphy of the New Jersey margin (Miller et al., 2005; Watts and Thorne, 1984) likely contain a dynamic component and cannot be used alone to derive eustatic sea-level curves. We will return to the way in which different models fit the New Jersey constraint in the Discussion.

Adjoint models. To overcome the shortcoming of backward advection in dealing with the conductive thermal boundary layer, adjoints of the governing equations can be introduced (Bunge et al., 2003; Ismail-Zadeh et al., 2004; Liu and Gurnis, 2008). The method consists in refining the assumed initial condition by iteratively solving the forward and backward adjoint and minimizing a cost function (in this case, between seismic tomography and forward-predicted temperature). The main advantage of this method is that it uses a forward solution of the coupled equations of motion and energy while converging upon the known tomographic structure. A significant limitation is the computational cost, as a single inversion requires multiple forward and adjoint solutions of the equations and is substantially more expensive than a backward run (Bunge et al., 2003; Table 2).

Liu et al. (2008) and Spasojevic et al. (2009) used a variant of this method to constrain mantle properties finding adjoint models (starting from present-day seismic tomography) that best fit the flooding history of North America (Fig. 5A). In these models, the mantle structure is scaled from shear wave seismic tomography (Grand, 2002) between 250 km and 2400 km depth, and the plate velocities of continuously closing plates (Gurnis et al., 2012). In addition, the flat-lying, continuous Cretaceous Farallon slab could only be backward advected (Fig. 5C) by imposing a viscous stress guide (Liu et al., 2008; Spasojevic et al., 2009). The dynamic topography (Fig. 5B) and its rate of change were used to predict the spatial pattern of marine inundation (Fig. 5A). This was then compared to the observed flooding in the Western Interior Seaway and to subsidence rates deduced from borehole data, respectively. The scaling factor between mantle tomography and density structure was constrained by minimizing the misfit between the amplitude of observed and predicted subsidence. The viscosity of the upper mantle was constrained by minimizing the misfit between observed and predicted rates of subsidence. The best fit between model output and subsidence-related observables was obtained for an upper-mantle reference viscosity of 1 × 10²¹ Pa s and a lower-mantle reference viscosity of 1.5 × 10²² Pa s (Liu et al., 2008; Spasojevic et al., 2009). Significantly, the inferred upper-mantle viscosity is within the range inferred from postglacial rebound studies (Cathles, 1975; Mitrovica and Forte, 2004). Finally, the best-fitting model from Liu et al. (2008) successfully reproduced the observed craton-ward migration of the prominent depocenter within the Cretaceous section of the interior of North America (Figs. 2C–2F).

Hybrid models.Spasojevic and Gurnis (2012) recently proposed a hybrid formulation in which the density structure of the lower mantle is entirely derived from backward advection of a seismic tomography model, whereas in the upper mantle, the positively buoyant anomalies are backward advected from seismic tomography and the negatively buoyant subduction zones are imposed from tectonic reconstructions and oceanic ages (Müller et al., 2008a). Spasojevic and Gurnis (2012) carried out a series of forward instantaneous flow models, hereafter referred to as SG12, from 90 Ma in 10 m.y. increment based on these density structures. They compared the predicted change in global dynamic topography for several models with different viscosity structures to the uplift and subsidence maps deduced from paleogeographic reconstructions (Blakey, 2008; Smith et al., 1994) and noted that different models compare best to data in different areas and suggested that the radial viscosity structure of the mantle is not laterally uniform.

This dynamic topography field of Spasojevic and Gurnis (2012; Fig. 4E) is similar to that proposed by Conrad and Husson (2009), which is expected, since model CH09 is instantaneous and model SG12 carries no information from the past to the present. Both models are also based on the same tomography model, S20RTS. The main difference between the two models is the amplitude of dynamic topography lows, which is ∼1.5 times larger in the model of Spasojevic and Gurnis (2012), which reflects (1) that the uppermost 300 km section of the mantle has been ignored in CH09, whereas the uppermost 250 km section has been removed in SG12, and (2) that the scaling factor to convert seismic velocity anomaly to density anomaly is depth independent and equal to 0.15 g cm^–3 km^–1 s in CH09, whereas it varies with depth between 0.05 and 0.3 g cm^–3 km^–1 s in SG12.

The good first-order agreement among four global residual topography fields (Fig. 1), derived independently and based on different assumptions, is encouraging. The most prominent and robust residual topography features are (1) the lows extending from central Europe to Southeast Asia, from eastern North America to western South America, the Australian-Antarctic Discordance (AAD), the Argentine Basin, and (2) the highs under southern and eastern Africa, Iceland, and the Darwin Rise.

The global dynamic topography fields predicted by five numerical mantle flow models (Fig. 4), based on two distinct methods and a hybrid of those, are also in overall agreement with each other, and with the residual topography (Fig. 1). All of these dynamic topography fields, with the exception of Steinberger (2007), are to some extent skewed toward negative values (Table 1), reflecting that negative dynamic topography anomalies are larger in amplitude but cover smaller areas than positive dynamic topography anomalies (since the mean of each field is zero). The wavelength of model dynamic topography is overall larger than that of proposed residual topographies (Fig. 1). The model of Steinberger (2007), based on a compilation of mantle tomography models, presents shorter wavelengths than models by Conrad and Husson (2009) and Spasojevic and Gurnis (2012) based on S20RTS (Ritsema et al., 2004) (Fig. 4). For instance, the dynamic topography high under Iceland and the Azores (Conrad et al., 2004) is distinct in the model of Steinberger (2007), subdued in the model of Conrad and Husson (2009), and absent in Spasojevic and Gurnis (2012) (Fig. 4). Dynamic topography lows in the equatorial South Atlantic Ocean are only captured by the model of Steinberger (2007), albeit with a larger amplitude than suggested by residual topography. With more robust processing techniques, especially with full waveform seismic tomography, future mantle tomography models will improve our understanding of global dynamic topography down to ∼1000 km scale, whereas present global studies are limited to a resolution of ∼3000–5000 km (Fig. 4).

Both synthetic time-dependent forward mantle flow models (Figs. 4A and 4C) and instantaneous models based on mantle tomography (Figs. 4B, 4D, and 4E) predict dynamic topography lows under the continents, with the notable exceptions of Africa south of the equator and of western North America (out of the five models, only that of Ricard et al. [1993] does not predict positive dynamic topography for western North America). Interestingly, all mantle flow models predict a dynamic topography low under Southeast Asia that is more pronounced than suggested by residual topography studies, and negative dynamic topography along the east coast of Asia, where the residual topography is positive (Figs. 1 and 4). One limit to constraining residual topography in Southeast Asia and along the east coast of Asia is the small extent of oceanic crust, limited to the South China, Sulu, Celebes, and Banda Seas in Southeast Asia and to the Sea of Japan and southernmost Sea of Okhotsk along the east coast of Asia (see COB contours in Fig. 1A). Because we have assigned a uniform elevation of –200 m to areas between oceans and continents, our proposed residual topography features a regional residual topography low in Southeast Asia and along the east coast of Asia (Fig. 1A), whereas other studies (Figs. 1B, 1C, and 1D) in which these areas are treated as continental suggest positive residual topography in the Gulf of Thailand, Borneo, and the east coast of Asia. The consistency of the dynamic topography low under Southeast Asia and the east coast of Asia in both instantaneous flow models (Figs. 4B, 4D, and 4E) and forward mantle flow models (Figs. 4A and 4C) suggests that the seismically fast mantle is reproduced by forward models. Given the poor control on residual dynamic topography in these regions, dynamic topography models for south and east Asia should be compared to time-dependent observables (e.g., Xie et al., 2006) and to mantle tomography in the case of synthetic forward models. However, Southeast Asia is the region with the largest extent of shallow marine seas, and such shallow seas are strongly correlated with predicted present-day dynamic topography lows (Gurnis, 1993b).

All mantle flow models predict positive dynamic topography for the entire Pacific Ocean (Fig. 4), which is at odds with residual topography studies, which consistently predict a large-scale, low-amplitude low in the northeast Pacific (Fig. 1). An anticorrelation between residual gravity and residual topography is known for this region (Crosby and McKenzie, 2009), which suggests a dynamic origin for the residual topography low, but there is no fast seismic anomaly structure under the northeast Pacific (Crosby and McKenzie, 2009). Indeed, Spasojevic et al. (2010) attributed the gravity lows under the northeast Pacific to mid- to upper-mantle upwelling above slab graveyards, thereby explaining the absence of large fast seismic velocity anomalies directly under this region and suggesting that the northeast Pacific could be presently uplifting. Perhaps this is a region with negative residual topography that is experiencing dynamic uplift (on its way up after being drawn down). Therefore, the mismatch between model dynamic topography and present-day residual topography in the northeast Pacific illustrates the limitation in using a snapshot of residual topography to constrain dynamic topography that occurs over time scales of several tens of millions of years. This underlines the concept that the time evolution of dynamic topography is arguably more important than its present-day sign and amplitude, and therefore that quality time-dependent proxies should be favored as model constraints.

Finally, the position of the AAD is well reproduced by instantaneous mantle flow models based on tomography (Figs. 4B, 4D, and 4E) but only at long wavelengths by forward models (Figs. 4A and 4C).

Arguably the most important consequence of dynamic topography studies is their influence on eustatic and relative sea level. In this section, we first discuss the effect of global dynamic topography on global sea level and then turn to the regional effects of dynamic topography.

Changes in dynamic topography could have a first-order influence on global sea level (Gurnis, 1992). Using a simple model of mantle flow associated with slabs, Gurnis (1993c) computed that the shifting motion of continents over the evolving long-wavelength dynamic topography would give rise to a 100–200 m eustatic change with two highs in the mid-Paleozoic and late Mesozoic. Husson and Conrad (2006) developed a model based on boundary layer theory to propose that changes in dynamic topography associated with geologically rapid (≤20 m.y.) changes in mean global tectonic velocities would only affect eustatic sea level by ∼20 m, which is relatively small in comparison to eustatic changes of several hundred meters associated with changes in the age distribution of the ocean floor (e.g., Hays and Pitman, 1973; Müller et al., 2008b; Spasojevic and Gurnis, 2012). Husson and Conrad (2006) proposed that the dynamic effect of longer-term (10⁸ yr) change in tectonic velocities on sea level could be up to ∼80 m.

Extrapolating from the present-day rate of change of dynamic topography, Conrad and Husson (2009) considered the effect on eustasy of continental motion across dynamic topography gradients during the aggregation and dispersal of a supercontinent. The most significant effect would occur during the dispersal phase, when fragments move away from the upwelling expected beneath a mature supercontinent (Coltice et al., 2007; Gurnis, 1988) toward dynamic topography lows associated with subduction zones. Conrad and Husson (2009) estimated that dynamic topography (Fig. 4D) currently triggers a net uplift of the ocean basins thereby contributing to a positive sea-level offset of ∼90 ± 20 m. In addition, they carried out a 1-m.y.-long forward model of mantle flow based on the present-day mantle structure and plate motions to propose a current rate of sea-level rise induced by dynamic topography of <1 m/m.y.

Spasojevic and Gurnis (2012) simultaneously calculated dynamic topography, the shifting of present-day continental isostatic topography with the plates, and the changing age distribution of the ocean floor, while forcing the water load to fit the evolving geoid in a discrete series of ten instantaneous mantle flow models between 90 Ma and present. Based on these calculations, they estimated that sea level has risen by between ∼110 m and ∼270 m since 90 Ma because of global dynamic topography, including ∼110–200 m caused by seafloor uplift and <70 m due to continental subsidence. Thus, dynamic topography influences global sea level by offsetting (by about one third) the decrease in sea level imposed by changes in the age distribution of the ocean floor. Dynamic topography should therefore not be neglected in studies of long-term global sea-level change. In computational models, it is possible to create a global sea-level curve by averaging the height of all continents with respect to sea level. Things are more complex in reality because large-amplitude swings in dynamic topography could have a disproportionate influence on the stratigraphically inferred record of sea-level change. This is particularly true on a regional scale, at which both dynamic topography and tectonic topography often dominate the stratigraphic signal.

Another important consequence of dynamic topography studies is that sea-level curves defined at a single “stable” passive continental margin (or indeed at any location) may not be representative of eustasy if the margin has been experiencing differential vertical dynamic motion. The most prominent example is the sea-level curve proposed by Miller et al. (2005) based on the careful backstripping of five wells from the New Jersey coastal plains, which predicts a much smaller decrease in sea level since the Late Cretaceous than other global sea-level curves (for details, see Müller et al., 2008b; Spasojevic et al., 2008). Recent dynamic topography studies all suggest that the New Jersey margin, although considered tectonically stable, has undergone significant vertical motions driven by mantle flow since the Late Cretaceous (Moucha et al., 2008a; Müller et al., 2008b; Spasojevic et al., 2008). Based on their backward advection model, described earlier herein, Moucha et al. (2008a) proposed that the New Jersey margin has been dynamically uplifting by ∼100–200 m over the past 30 Ma. In contrast, Spasojevic et al. (2008) argued for dynamic subsidence of the New Jersey margin by >300 m since 55 Ma, including ∼100–200 m over the past 30 Ma, based on the adjoint model of Liu et al. (2008) described earlier. Finally, Müller et al. (2008b) estimated that the New Jersey margin has subsided by 105–385 m since 70 Ma using backward advection models based on three different tomography models. Müller et al. (2008b) considered two types of models: “pure” backward advection models and “modified” backward advection models in which upwellings are continued upward to 220 km depth and downwellings are removed from the uppermost 220 km. For the past 30 Ma, their “pure” backward advection models predict between ∼50 m and ∼150 m of margin subsidence, while their modified backward advection models predict between 0 and ∼70 m of margin subsidence, and ∼90 m of subsidence between 30 Ma and 8 Ma followed by 40 m of uplift since 8 Ma for the tomography model of Grand (2002).

Comparison between observed tectonic subsidence and predicted dynamic topography for well COST-B2. To shed some light on this controversy, we focus on the tectonic subsidence history of the COST B-2 well offshore New Jersey (Scholle, 1977), which we compare to the evolution of dynamic topography predicted by published models and by the model presented in this study.

Subsidence of well COST-B2.Figure 6 shows the tectonic subsidence of well COST B-2 with respect to present day obtained by backstripping the stratigraphy (courtesy M. Kominz) for two end-member cases of sea-level change. Well COST-B2 shows a marked acceleration in postrift subsidence at ca. 16 Ma (Fig. 6A; Appendix 1) that is also evident in strain rate inversion (Fig. 6B; Appendix 1) and anomalous subsidence (Fig. 6C; Appendix 1) analysis. In addition, both strain rate inversion and anomalous subsidence analysis reveal a phase of anomalous subsidence between 85 and 15 Ma (Figs. 6B and 6C). The sign of this anomalous subsidence depends on the eustatic sea-level curve used in backstripping, with changing (constant) sea level resulting in subsidence (uplift) (Figs. 6C). The hypothesis of constant sea level is difficult to reconcile with the growth of ice sheets that would have triggered a sea-level drop by 54 m over the past 33 m.y. (Miller et al., 2005), and with the increase in the volume of ocean basins that would have caused a sea-level drop of at least 120 m since 85 Ma (Müller et al., 2008a; Spasojevic and Gurnis, 2012). This suggests that well COST-B2 has been subsiding rather than uplifting over the past 85 m.y., which is in agreement with the subsidence of New Jersey paleoshorelines by 50–200 m since the Eocene (Spasojevic et al., 2008). Given the time scale (tens of millions of years), this anomalous subsidence could result from dynamic subsidence induced by long-wavelength mantle flow.

Comparison with the dynamic topography predicted by the model from this study.Figure 7 shows the evolution of dynamic topography predicted by our forward model at 151 Ma, 100 Ma, 49 Ma, and present day, along with the coastlines and subduction zones reconstructed through time with respect to the mantle. The predicted degree two pattern of the dynamic topography field is stable through time, and in this subduction-driven model, the main changes in dynamic topography are the location, width, and amplitude of dynamic topography lows. Well COST-B2 (Figs. 1, 4, and 7) is initially located over a dynamic topography high and drifts toward the broadening dynamic topography low associated with the old Farallon slab as North America drifts to the west (Figs. 7B and 7C). Our model predicts that well COST-B2 is currently located on the northeast edge of a mild dynamic topography low (Fig. 7D), which is consistent with other dynamic topography models (Fig. 4) and in broad agreement with estimates of residual topography (Fig. 1), in which the topography low is more pronounced. In contrast to well COST-B2, our model predicts that the synthetic well located on the conjugate margin in Western Sahara would have remained on a dynamic topography high since 200 Ma (Fig. 7). For well COST-B2, the predicted evolution of dynamic topography shows an initial uplift by ∼150 m in two stages (200–170 Ma and 140–120 Ma), followed by 850 m of dynamic subsidence since 120 Ma (Fig. 8). This predicted dynamic subsidence of well COST-B2 since 120 Ma is consistent with the anomalous strain and with the anomalous subsidence discussed previously herein (Figs. 6 and 9). The model predicts 580 m of dynamic subsidence since 85 Ma, which is ∼30% greater than the anomalous subsidence deduced from backstripping using the sea-level curve of Haq et al. (1987; Fig. 9; Appendix 1,) and only 30 m of subsidence over the past 15 Ma, which is four to six times smaller than the anomalous subsidence for that period (Fig. 9; Appendix 1). As for the conjugate margin, our model predicts that it was uplifted by 250 m from 200 Ma to 110 Ma and that is has subsided by 170 m since 90 Ma (Fig. 8). This predicted subsidence is five times less than that of well COST-B2 (Fig. 8), which confirms that the margin of northwest Africa has not been subject to as much postrift vertical motion as the New Jersey margin (Moucha et al., 2008a). Unfortunately, little data are available to estimate the vertical motions of northwest Africa, particularly for the conjugate margin of New Jersey. Nevertheless, the subsidence of northwest Africa predicted by our models is to first order compatible with the marine transgression since 40 Ma deduced from paleogeographic analysis for that area (Spasojevic and Gurnis, 2012) and with the moderate anomalous subsidence of the central Atlantic margin in the Senegal Basin (Brun and Lucazeau, 1988). In addition, the relative vertical stability of northwest Africa above a dynamic topography high predicted by our models is in agreement with the thin onshore Cenozoic stratigraphic thickness (∼70 m) of the Senegal Basin (e.g., Swezey, 2009).

The dynamic subsidence inferred over tens of millions of years for well COST-B2 and from paleoshoreline analysis in New Jersey appears to be in contradiction with recent uplift reported farther south on the east coast of North America (Rowley et al., 2011). Dowsett and Cronin (1990) reported recent uplift along an early Pliocene wave-cut scarp in South Carolina and North Carolina (see Fig. 7 for location) that Rowley et al. (2011) extended from Georgia to Virginia. Together, these studies suggest up to ∼50–65 m of uplift since ca. 3.5–3.0 Ma over ∼500 km in an area located several hundred kilometers south of well COST-B2. The evolution of dynamic topography predicted by our model at one of the sections reported by Dowsett and Cronin (1990) shown on Figure 8 is overall similar to that predicted for well COST-B2. An important difference, however, is that the North Carolina location is predicted to have been dynamically uplifted by ∼20 m over the past 25 Ma, a period during which well COST-B2 is predicted to have subsided by ∼50 m. This result underlines the difficulty of identifying the evolution of dynamic topography in the geological record, since a post–mid-Pliocene uplift of relatively small amplitude could be more prominent in the present-day landscape than the dynamic subsidence of several hundred meters predicted from the early Aptian to the late Eocene. Determining spatial and temporal trends of dynamic topography in details along the east coast of North America requires further work.

Comparison with the dynamic topography predicted by previous models. To put our results for well COST-B2 in perspective, we compare them in Figure 9 with the results of three previous studies of dynamic topography for the same area, all using the same definition of mantle-flow–driven dynamic topography, excluding the dynamics of the thermal boundary layer (Moucha et al., 2008a; Müller et al., 2008b; Spasojevic et al., 2008). We use the results of Moucha et al. (2008a) for surface dynamic topography calculated using mantle density heterogeneities deeper than 200 km. These extend back to 30 Ma and were integrated over a ∼450 km × 550 km area offshore New Jersey that includes the location of well COST-B2. The results of Spasojevic et al. (2008) extend back to 50 Ma and were sampled at a point location ∼150 km west of well COST-B2, and the results of Müller et al. (2008b) extend back to 70 Ma and were sampled at the location of well COST-B2. Given the long-wavelength nature of dynamic topography, the results of each of these studies can be considered representative of the New Jersey margin and compared to our results. Two cases are presented for each model considered, including a second case for our model in which the lower-mantle reference viscosity is 50 times larger than that of the upper mantle (this ratio is 100 in the reference case). Modeling results are also compared to the anomalous subsidence deduced from Figure 6C. Note that this direct comparison implies that postrift subsidence is only affected by cooling of the lithosphere and dynamic topography, which may be an oversimplification. The Neogene phase of anomalous subsidence (Fig. 9) is difficult to reconcile with the results of Moucha et al. (2008a) and with model “ngrand m.b.a” of Müller et al. (2008b), which predict uplift since 30 Ma and since 10 Ma, respectively. All other models predict subsidence throughout the Neogene, although the amplitude is too small except for case 2 of Spasojevic et al. (2008). The models of Moucha et al. (2008a) predict uplift between 30 Ma and 10 Ma, which is consistent with backstripping for constant sea level, whereas all other models predict subsidence prior to 15 Ma, which is consistent with backstripping using the sea-level curve of Haq et al. (1987). In addition, the uplift predicted by the model of Moucha et al. (2008a) is consistent with the geologically inferred uplift from Virginia to Georgia (Dowsett and Cronin, 1990; Rowley et al., 2011), although the inferred uplift rate is three to six times larger than the modeled uplift rate. As discussed earlier herein, sea level has fallen by at least 170 m over the past 85 Ma because of a combination of decreasing volume of ocean basins and ice-sheet growth in the Cenozoic. As a consequence, we favor models that predict Cenozoic dynamic subsidence of the New Jersey margin. However, the large uncertainty associated with the dynamic models that predict between 160 m and 500 m of subsidence over the past 50 Ma (Fig. 9) makes it difficult to estimate long-term eustatic sea level from the high-quality record of the New Jersey margin by correcting for dynamic topography (Kominz et al., 2008; Müller et al., 2008b). Long-term eustatic sea-level change may be estimated from the global stratigraphic record (Haq et al., 1987), a synthesis of the flooding record from continental platforms around the world (Bond, 1978), or by modeling the change of volume of ocean basins through time (Kominz, 1984; Müller et al., 2008b; Spasojevic and Gurnis, 2012), and should be corrected for the global effect of dynamic topography (Conrad and Husson, 2009; Spasojevic and Gurnis, 2012).

The contrasting evolution of dynamic topography predicted by different models for the New Jersey margin (Fig. 9) calls for explanation. The adjoint model of Spasojevic et al. (2008) and the forward model presented here are based on numerical approaches different from the backward advection models of Müller et al. (2008b) and of Moucha et al. (2008a). In addition, lateral viscosity variations that influence the evolution of thermal heterogeneity in the mantle (Moucha et al., 2007) are taken into account via the temperature dependence of viscosity in the model presented here and in that of Spasojevic et al. (2008), whereas they are ignored in the models of Moucha et al. (2008a) and of Müller et al. (2008b), which only have radial variations in viscosity. Furthermore, the present model and that of Spasojevic et al. (2008) include a simple two-layer viscosity structure for the mantle (Paulson et al., 2007), whereas the models of Moucha et al. (2008a) and of Müller et al. (2008b) each consider ∼20 viscosity layers throughout the mantle (Mitrovica and Forte, 2004; Steinberger and Calderwood, 2006). These different viscosity structures affect the depth dependence of surface dynamic topography via the viscosity ratio between the different layers (Hager and Clayton, 1989; Fig. 3; Eq. 2). For instance, decreasing the viscosity of the lower mantle by a factor of two in our model increases the predicted subsidence since 85 Ma by ∼15% (Fig. 9).

The models of Moucha et al. (2008a) and of Müller et al. (2008b) are overall similar in design, yet they give different results for the New Jersey margin. The main differences between these two backward advection models include the mantle tomography employed, the conversion from seismic velocity to density, the viscosity structure of the mantle and the treatment of thermal boundary layers—essentially the four main steps of estimating dynamic topography (Braun, 2010). The scaling factor from seismic velocity to density clearly has an effect on the amplitude of the predicted absolute dynamic topography. A constant scaling was applied to the whole mantle by Müller et al. (2008b), but a radially and laterally variable scaling in the tomographic inversion of Simmons et al. (2007) was used by Moucha et al. (2008a). The effect of this spatially variable scaling factor on the spatio-temporal evolution of dynamic topography predicted by Moucha et al. (2008a) is likely to be significant. As for mantle tomography itself, Müller et al. (2008b) predicted Tertiary subsidence of the New Jersey margin using three different S-wave tomography models. The results of Müller et al. (2008b) also confirm that the treatment of thermal boundary layers is important in backward advection models (Conrad and Gurnis, 2003), as shown by the difference between their “pure” and “modified” models in Figure 9. The “modified” model predicts 150 m less dynamic subsidence during the Tertiary than the “pure” model, presumably due to the removal of shallow downwellings in the “modified” model. The “modified” model also predicts uplift for the past 10 m.y., likely associated with the upward continuation of upwellings, at odds with the subsidence of the New Jersey margin. The structure of the lithosphere is thus key, not only in estimating present-day residual topography, as discussed previously, but also in modeling time-dependent dynamic topography. Finally, the main difference between the viscosity structures used by Moucha et al. (2008a) and Müller et al. (2008b) is the 30-km-thick viscosity notch just above the 670 km phase transition (Mitrovica and Forte, 2004) for which the viscosity is 300 times (5 times) lower than the reference viscosity in model V1 (model V2). This narrow low-viscosity zone, potentially related to transformational superplasticity or to an internal thermal boundary layer (Pari and Peltier, 1995), should at least partially decouple the upper mantle from the lower mantle. Such a decoupling should act to suppress the effect of the subducting Farallon slab, which is imaged close to the top of the lower mantle under the present-day east coast of North America (Liu et al., 2008; Müller et al., 2008b). Thus, Cenozoic subsidence of the New Jersey margin, as opposed to the uplift predicted by Moucha et al. (2008a), appears mainly dependent on the absence or presence of a low-viscosity notch just above the lower mantle under North America.

The previous section illustrates the need for detailed comparisons between the evolution of dynamic topography predicted by models and that deduced from the geological record, as well as the complexity and uncertainty of dynamic topography models. The comparison of model predictions to the geological record is made complex by the large-wavelength nature of the dynamic topography signal and the relative scarcity of high-quality constraints. The New Jersey example discussed here illustrates the importance of the assumptions made when interpreting the geological record. Nevertheless, such comparisons are valuable since they constitute constraints for the viscosity structure of the mantle, and for the response of the lithosphere to mantle flow.

From an observational point of view, increasing amounts of geophysical and geochemical data and their compilation in databases will over time improve the knowledge of the thermal and chemical structure of the crust and lithosphere, which are critical to constraining present-day dynamic topography. Global databases of stratigraphic data will also be important to re-appraise global sea-level change and constrain the evolution of dynamic topography. From a modeling point of view, the comparison between the evolution of topography deduced from the geological record and the dynamic topography predicted by mantle flow models can be made more direct. The usual, indirect, approach consists of comparing the model topography, after removing the isostatic component (thermal boundary layer), to the estimated dynamic component of the topography deduced from the geological record. Increasing computing power will soon make it possible to explicitly model the deformation of the continental lithosphere in global mantle flow models. This would permit the simultaneous estimation of isostatic and dynamic topographies, and would facilitate the comparison between modeling results and data. Explicitly modeling both the continental and oceanic lithospheres in global mantle flow models would also open up the opportunity to study the response of the viscoelastic lithosphere to mantle stresses. Because of the computational cost of the finite-element viscoelastic treatment of global mantle flow, viscoelastic studies of sea level have largely been limited to the postglacial rebound problem (e.g., Paulson et al., 2007). Increasing computational power and new numerical methods will make calculations possible that include the full range of rheological and structural variations from hundreds of meters to global scales (Stadler et al., 2010). This would improve our understanding of topography at intermediate wavelengths (400–800 km) at which both small-scale convection and elastic flexure occur.

Considerable progress has been made over the past 30 yr on the observational constraints on and the modeling of dynamic topography. In addition to constraining the nature of mantle convection and mantle rheology, studies of dynamic topography have led to a new interpretation of the stratigraphic record and of long-term sea-level change. Current global models and observations of dynamic topography are in broad agreement at very large scale (∼10,000 km) but differ in some areas and at smaller scales (<5000 km). Increasing computing resources and modeling techniques provide the opportunity to develop higher-resolution models that will be more readily comparable to geological data. In parallel, acquiring a better understanding of the global structure of the lithosphere through data acquisition and compilation would improve the constraint on present-day global dynamic topography and therefore on mantle properties.

The updated stratigraphy of well COST-B2 was backstripped using the software of White (1993) for constant sea level and for the eustatic curve proposed by Haq et al. (1987) that features sea-level variations of up to 250 m since the Paleozoic. Both tectonic subsidence curves show ∼2.5 km of tectonic subsidence during the main rifting phase between 210 Ma and ca. 175 Ma, followed by ∼1.5 km of thermal subsidence. Considering the large error bars on paleo–water depth associated with the two anomalies in thermal subsidence at ca. 85 Ma and ca. 35 Ma, the thermal subsidence phase is overall smooth until ca. 16 Ma, when there is a marked acceleration in subsidence.

In addition to backstripping, we carried out a strain rate inversion analysis of the backstripped data following the method of (White, 1993). The best-fit subsidence curves obtained from the inversion are shown in Figure 6A, and the strain rate history is shown in Figure 6B for both sea-level curves. In both cases, the strain rate history describes a main first rifting phase between 210 Ma and ca. 165 Ma, followed by a relatively minor rifting phase between 155 Ma and 135 Ma. The inversion predicts an increase in strain rate from ∼15 Ma, of amplitude larger than the apparent second rift phase at ca. 145 Ma (Fig. 6B). In addition, the inversion using the eustatic curve of Haq et al. (1987) predicts protracted strain between ca. 85 Ma and 15 Ma.

We carried out forward modeling of pure shear lithospheric stretching following the method of Jarvis and McKenzie (1980) in order to estimate the anomalous subsidence at well COST-B2. These forward models do not take into account flexure or depth-dependent extension. We assumed a single rifting phase between 210 Ma and 175 Ma, lithosphere initially 160 km thick, continental crust initially 31 km thick, and varying stretching factors between 1.8 and 2.2. Although the solution is nonunique, these parameters give a satisfactory fit to the observed tectonic subsidence (Fig. 6A).

We estimated the amount of anomalous subsidence for well COST-B2 by subtracting the thermal subsidence predicted by the forward models from the observed data (Fig. 6C). For each sea-level curve, we used the model that best fits the early stages of postrift subsidence (130–110 Ma), namely, a stretching factor of 1.9 for the curve of Haq et al. (1987) and of 2.0 for constant sea level (Fig. 6A). On Figure 6C, positive (negative) values reflect anomalous subsidence (uplift) for which the model subsidence is shallower (deeper) than the backstripped data. Note the large error bars associated with paleo–water depth, particularly for the period between 85 Ma and 35 Ma, for which the uncertainty is ±230 m (Fig. 6). We limited our analysis to the relative changes in anomalous subsidence that exceeded this uncertainty. The case with eustatic sea-level change displays 175 ± 70 m of relative uplift between 125 Ma and 85 Ma, and then 400 ± 45 m of relative subsidence since 85 Ma, including 230 ± 110 m since 15 Ma. The case for constant sea level shows 115 ± 70 m of relative uplift between 125 Ma and 85 Ma and 130 ± 65 m of relative subsidence since 10 Ma. Note that in this latter case, the relative subsidence of 60 ± 45 m since 85 Ma is of smaller amplitude than the Neogene phase of subsidence, which implies anomalous uplift between 85 Ma and 10 Ma. The large differences in the evolution of anomalous subsidence for the case with constant sea level and for the case accounting for sea-level change illustrate the importance of the sea-level curve used when backstripping stratigraphic data.

We thank C. Conrad, Y. Ricard, S. Spasojevic, and B. Steinberger for sharing their results, M. Kominz for the updated COST-B2 well data, and L. Husson for discussions. This paper benefited from comments by two anonymous reviewers. Flament was supported by Statoil, Gurnis by the National Science Foundation (through EAR-0810303) and the Caltech Tectonics Observatory (by the Gordon and Betty Moore Foundation), and Müller by Australian Research Council grant FL0992245. Computer simulations were carried out on the Sun Constellation VAYU Cluster of the Australian National Computational Infrastructure. Figures 1, 4, and 7 were prepared using the Generic Mapping Tools (Wessel and Smith, 1998). The original CitcomS software was obtained from CIG, Computational Infrastructure for Geodynamics (http://geodynamics.org). This is contribution 216 of the Caltech Tectonics Observatory.